Abstract: A mode-superposition-based equivalent linearization method was proposed based on the mode perturbation method and mode superposition method for fast seismic response analysis of two-dimensional sites with non-proportional hysteresis damping. In this method, based on the modal analysis results of the original system, the eigenvalue of the new system can be obtained by the modal perturbation method with modified iterative stiffness. Based on the modal superposition method, the dimension of matrices in the motion equation of a non-proportional complex damping system is reduced, and the modal-superposition-based equivalent linearization method of generalized coordinate coupling and decoupling are formed. Then, the influence of coupling coefficients in a damping matrix is discussed. The variation law of the horizontal acceleration amplification factor under different input ground motion amplitudes is studied by an engineering site in Hangzhou. And the results are compared with the two-dimensional frequency domain equivalent linearization to verify the accuracy and effectiveness of the proposed method. The numerical analysis results show that: the calculation error is less than 5% under different input amplitudes with the generalized coordinate coupled method. But by the generalized coordinate decoupling method, the calculation errors increase with the increase of the coupling factor. Hence, when the ground motion amplitude is less than 0.1 g, the generalized coordinate decoupling method can be used to calculate the site seismic response. For the computational efficiency, the iterative times of mode-superposition-based equivalent linearization method is basically the same with that of the two-dimensional frequency-domain equivalent linearization method, but the computational time is about half that of the frequency-domain equivalent linearization method.