The power spectrum method is widely used for its high efficiency in reconstructing high-frequency components of turbulence. However， the power spectrum method suffers from a significant limitation： it fails to adequately sample low-frequency components due to the finite size of the simulation aperture. This inherent constraint leads to a severe truncation of the energy concentrated at large scales， resulting in a discrepancy between the simulated phase screens and the physical reality of atmospheric turbulence. Consequently， there exists a contradiction between the physical modeling requirements， which necessitate the inclusion of large-scale eddies， and the numerical implementation constraints. To address this long-standing challenge of low-frequency deficiency in the power spectrum method， this paper proposes a hybrid simulation method that leverages the frequency-domain complementarity between the Zernike polynomial method and the power spectrum method. The primary objective of this research is to develop an algorithm that can accurately reconstruct the low-frequency components missing in the power spectrum method while maintaining high computational efficiency， thereby resolving the conflict between physical fidelity and numerical feasibility.The proposed method is grounded in a rigorous theoretical analysis of the frequency-domain response characteristics of Zernike polynomials and the sampling mechanism of the power spectrum method. By deriving the analytical relationship between the Zernike mode response functions and the Modified Von Kármán （MVK） turbulence spectrum， we establish a quantitative criterion for seamless frequency-domain stitching. The core methodology decomposes the total phase screen into two distinct components： a low-frequency component and a high-frequency component. The low-frequency component， which corresponds to the energy peak of the MVK spectrum typically located below the minimum sampling frequency of the power spectrum method， is precisely reconstructed using a finite number of low-order Zernike polynomials. Through numerical analysis of the Zernike response function peaks， we determine that the first six Zernike polynomials （excluding the piston term） are sufficient to cover the low-frequency blind spot. This Zernike-reconstructed low-frequency phase is then linearly added to the high-frequency phase component generated by the standard power spectrum method. This hybrid approach ensures that the transition between the modal expansion （Zernike） and the spectral sampling （the power spectrum method） is smooth and physically consistent， effectively bridging the gap between the two methods.To validate the performance of the proposed method， extensive numerical experiments were conducted. Statistical analyses were performed on 5，000 phase screens. The results demonstrate that the generated phase screens exhibit rich high-frequency details as well as distinct large-scale low-frequency undulations， visually aligning with the expected physical behavior. Quantitatively， the calculated phase structure functions show excellent agreement with the theoretical MVK model across the entire frequency range， effectively correcting the statistical deviations observed in traditional PSRM at low frequencies. Furthermore， the variance of the Zernike coefficients derived from the simulated screens strictly adheres to the theoretical statistical laws predicted by Noll's theory. Comparative studies regarding computational efficiency reveal that the proposed method is significantly better than traditional subharmonic compensation techniques. By eliminating the need for multiple Fast Fourier Transforms （FFTs） required by subharmonic methods， the computational complexity is highly reduced. The tests also indicate that the proposed method is approximately 10 to 20 times faster than the power spectrum method with subharmonic compensation， particularly for high-resolution simulations.In conclusion， this study successfully resolves the low-frequency component deficiency inherent in the power spectrum method through a hybrid strategy based on frequency-domain complementarity. The proposed algorithm provides a robust solution that reconciles the physical requirements of atmospheric turbulence modeling with the practical constraints of numerical simulation. By utilizing the first six Zernike polynomials to compensate for the missing low frequencies， the method achieves high-fidelity simulation results without introducing redundant calculations or complex iterations. This approach not only ensures the statistical accuracy of the phase screens， as evidenced by the structure functions and coefficient variances， but also provides a practical and optimal solution for applications that demand massive datasets or real-time processing capabilities in optical engineering and atmospheric physics.