Full-waveform inversion (FWI) imaging is a high-resolution seismic imaging technique that directly produces subsurface images by inverting the full recorded wavefield. However, its reliability is often limited by numerical dispersion errors arising from finite-difference (FD) forward modeling. One key approach for reducing dispersion is to optimize the FD coefficients using an optimization algorithm. However, conventional methods for optimizing FD weights focus only on reducing spatial dispersion, which can weaken numerical stability, especially when using large time steps (i.e., high Courant–Friedrichs–Lewy [CFL] numbers). To address this issue, we introduce a new optimization approach that improves both simulation accuracy and stability. The proposed method combines error functions from both the time–space domain and the spatial domain into a single adaptive objective function. A dynamic weighting factor, which depends on the CFL number, facilitates a trade-off between accuracy and stability of the optimal FD weights. We also use the seismic wavelet spectrum as prior information to constrain the optimization. The optimization problem is solved by the least-squares method. In the theoretical test, the proposed weights significantly reduce wavefield simulation errors across a wide range of wavenumbers, with a higher CFL number than conventional approaches. When applied to FWI, this method reduces phase distortion and local minima in the objective function. In a test using the Marmousi model at 40 Hz, our approach produced clear and continuous deep structures, closely matching results from dispersion-free benchmarks. In contrast, conventional methods failed due to severe dispersion. This work provides a more robust numerical foundation for high-frequency FWI imaging by improving both accuracy and stability.