Introduction: The dynamic response of structures subjected to earthquake ground motion is usually calculated by the direct integral method in time domain and response spectrum method in frequency domain. For simple structures, the direct integral method is widely used and accuracy requirement can also be achieved. However, for complex MDOF systems, the computational time of the direct integral method is remarkably long. Therefore, the response spectrum method is frequently used to estimate the peak value of structural seismic response in engineering practice. In recent years, the response spectrum method is extended to calculate the peak response of the structures with different characteristics such as long-span structure subjected to multi-support seismic excitations (Hao 1989; Kiureghian and Neummhofer 1992; Su et al. 2006; Tian and Li 2010; Tian et al. 2014), non-classically damping system (Foss 1958; Yu and Zhou 2008; Xu 2010; Guo et al. 2011, 2013; Liu et al. 2015) and the structure considering ductility and strain rate effects (Li and Li 2010, 2012; Li et al. 2013). However, for the MDOF system subjected to multi-support excitations, current studies on the dynamic analysis in frequency domain are usually focus on linear structures. Therefore, the response spectrum method is inapplicable for the nonlinear system such as the structure equipped with nonlinear base isolations (Lu et al. 2005; Xu et al. 2009). In the present paper, an effective approach is proposed to investigate PDF for the response of stochastic MDOF system with nonlinear elastic supports under multi-support excitations. It is known that the PDF solution of nonlinear stochastic dynamic system is governed by FPK equation and the exact solution of this equation is very difficult to obtain. Only under some restrictive conditions, the stationary exact solutions are obtainable for some one, two, or few-dimensional system. Therefore, some methods were introduced or extended to obtain the approximate probabilistic solutions of nonlinear stochastic dynamic systems. The most frequently used approximation method is the equivalent linearization or Gaussian closure procedure (Iyengar and Dash 1978). If the structure is highly nonlinear or the multiplicative noise is involved, the probability distribution of the system response is usually far from being Gaussian. In such a case the equivalent linearization and Gaussian closure methods are considered inapplicable. A non-Gaussian closure method (Assaf and Zirkie 1976; Grandall 1980) is proposed using Cram-Charlier series or Hermite polynomial to approximate the PDF solution. This method improves the accuracy of the approximate solution but may lead to negative probability. Stratonovich (1963) proposed the stochastic average method to calculate the PDF solution of response amplitudes of nonlinear systems. This method is merely suitable for the weakly damped systems subjected to weak excitations. Sobczyk and Trebicki (1990) extended the maximum entropy method to the general class of stochastic nonlinear systems and apply it to stochastic nonlinear oscillators. Numerical analysis indicates that the maximum entropy principle can serve as a general method for determining the probability distributions of oscillators. However, it is difficult to extend for MDOF problems.  
    In the present paper, the structural dynamic equation in matrix form is firstly decoupled by mode decomposition method and then reconstructed using state-space technique (Kalman 1959). From the reconstructed SDOF equation describing the modal response of structure, the FPK formula is deduced and solved by EPC method. According to the PDFs of the modal responses, the PDF for the total response is calculated based on probability theory. The advantage of the proposed approach is that due to the application of the traditional mode decomposition method, the MDOF structure is conveniently transformed into SDOF systems so that the problems in solving the MDOF FPK equation are avoided. Moreover, the probabilistic solution is more reasonable and scientific than the peak response solution of the conventional response spectrum method. A two-span bridge is taken as an example to calculate the PDF of structural response by Monte Carlo and the proposed methods, and numerical results show that the proposed method is accurate and efficient.

     For the seismic analysis of spatial system with underground structure, one may input the multi-point earthquake excitations at the base of the building. In this case, the underground excitation is generated by he seismic wave propagating to the base of the underground structure at a certain depth below the ground.The existing ground PSD, response spectrum and coherence models cannot be directly employed to simulate the multi-support underground motions. Therefore, the underground PSD and underground response spectrum are firstly deduced and then the underground cross-correlation function is further presented. These three theoretical models are firstly employed to establish the total power spectrum matrix and then the underground motions are generated by decomposing the aforementioned matrix using Cholesky decomposition method. Moreover, these three theoretical models are not only necessary conditions to simulate the multi-support earthquake underground motions but also the standards to verify the simulation results. A visual program for simulating multi-support earthquake underground motions is developed and its reliability is verified. Finally, a two-span bridge is taken as an example to calculate its seismic response under both multi-support earthquake ground and underground motions. Numerical results show that the simulation of multi-support earthquake underground motions is significant for both study and engineering application.