Abstract: Since ancient times, earthquakes have been an unavoidable natural phenomenon. As one of the regions with frequent seismic activities, China urgently needs advanced underground detection technology to characterize the fine structure of the lithosphere. This will provide a better understanding of the generation and development laws of earthquakes and minimize their impact on human society. Full waveform inversion （FWI） is an inversion method that utilizes the dynamic characteristics of seismic wave propagation to obtain subsurface physical parameters, providing important insights into the fine structure of the subsurface. This method further enriches the theoretical framework of conventional velocity modeling methods and has gradually become a research hotspot due to its high accuracy, multi-parameter, and multi-dimensional modeling advantages.　　As a high-resolution seismic inversion method, FWI minimizes the objective function that measures the discrepancy between the synthetic waveform and seismic observation data through wavefield simulation and optimization iteration, obtaining optimal subsurface medium parameters and constructing imaging. In the early stages of research, FWI was primarily performed in the time domain. Subsequently, researchers extended the inversion theory to the frequency domain. In comparison, frequency-domain inversion methods are conducive to parallel computing, reducing computational burden, and easy characterization of attenuation effects. By selecting a small number of frequencies that sufficiently cover the wavenumber, the fine structural details of the inversion model can be characterized without compromising inversion accuracy, thereby further reducing the computational cost of frequency-domain FWI. Additionally, frequency-domain wavefield simulation exhibits higher efficiency, especially in the case of multiple sources and increased attenuation factors. Therefore, conducting FWI in the frequency domain holds significant importance.　　With the continuous development of seismic research and rapid improvement in computer hardware, traditional acoustic wave equations have difficulty meeting the requirements of fine subsurface detection. FWI has gradually evolved towards the direction of elastic parameter inversion based on isotropic, anisotropic, viscoelastic wave equations, among others. Compared to acoustic waves, elastic waves include the propagation laws of P-waves and S-waves, providing richer wavefield information. Furthermore, in the near-surface region where the propagation time is short, it is challenging to separate body waves and surface waves, and attenuating or filtering surface waves is not easy. Both types of waves need to be considered for imaging. The elastic wave equation accurately simulates seismic wave propagation in subsurface media. Therefore, this study adopts the elastic wave equation as a mathematical model to simulate the propagation laws of seismic waves and conducts corresponding research on inversion methods.　　 In FWI, the main computational burden lies in the forward modeling process, and the computational efficiency and accuracy of the inversion heavily depend on the quality of the forward simulation method. Considering the complex structure of the elastic wave equation, traditional numerical solving methods require more computational costs. Therefore, employing efficient forward simulation methods is crucial for improving computational efficiency and the accuracy of inversion results. The nearly analytie discrete （NAD） operator has emerged as an important class of differencing operators in recent years. Its construction principle involves approximating high-order partial derivatives using the cross-combination of nodal displacements and their gradient values. The NAD operator offers high accuracy and significantly enhances computational efficiency. Therefore, in order to improve computational efficiency and the accuracy of inversion results, this study introduces the NAD operator for frequency-domain elastic wave equation forward simulation.　　This study establishes a frequency-domain elastic wave FWI framework based on the NAD operator. It derives the gradient calculation formula for the inversion objective function using the sparse block structure of the impedance matrix. The idea of separating compressional and shear waves is implemented throughout the entire solution process, greatly reducing the scale of gradient calculation and suppressing interference between them, further improving inversion computational efficiency. Finally, numerical experiments using three classic media models are conducted, obtaining ideal inversion results. The effectiveness and correctness of the proposed elastic wave inversion method are validated through the trend curves of inversion errors at different frequencies.　　In summary, the frequency-domain elastic wave FWI method based on the NAD operator derived in this study exhibits significant advantages in the inversion of subsurface physical parameters. Through numerical experiments, we have demonstrated the effectiveness and accuracy of this method in capturing fine subsurface structures. We believe that this method has broad application prospects in geological exploration and resource development and provides valuable references for researchers in related fields.