Herein, a method is developed, within the framework of the multivariate linear prediction (MLP), for conditionally simulating spatially correlated nonstationary ground motions at unobserved points measured in this local region. Given the limitation of recognizing the characteristics of the original signal, this scheme is based on the extension of the MLP to the wavelet packet transform (WPT). In this process, the WPT algorithm is used to estimate the corresponding wavelet packet coefficients from the decomposed orthogonal component time-histories corresponding to nonoverlapping frequency bands, which can accurately describe the time-varying nonstationary characteristics of the original signals. The original signal used in the MLP is replaced; instead, the wavelet packet coefficients are regarded as signals and are decomposed into different nonoverlapping subfrequency bands. Subsequently, this simulating process is transformed into a solution of the wavelet packet coefficients at the target points. Finally, an inverse WPT is added to obtain the corresponding ground motions. In addition to matching the specified coherency values, the earthquake time-histories generated by the proposed scheme inherit the physical properties of the predefined recordings well, including the waveform, cumulative energy, response spectrum, and nonstationary characteristics. Two numerical examples and blind tests are employed to simulate an ensemble of spatially correlated ground motions using data from the SMART-1 array. The proposed method will have a significant influence on engineering applications and regional hazard analyses.