## 1 Introduction  
Elastic wave numerical simulation is an important means by which to study the law of seismic wave propagation within the ground, and it plays a vital role in the whole process of seismic exploration, including data acquisition, processing, and interpretation. Many numerical simulation methods, including the finite element [1, 2], finite difference [3, 4, 5], spectral element [6], and pseudo spectral [7] methods have been successfully applied to elastic wave forward modeling. The staggered grid finite difference method (SGFD) is widely used in the forward modeling of elastic wave equations owing to its high calculation efficiency, high accuracy, and convenient implementation process. However, there are two problems with SGFD. One is that the fixed grid step may discretize the interface to the wrong depth, resulting in an inaccurate travel time for the reflected waves. The second is that a stepped interface will appear when the undulating interface is discretized with a regular grid, which will generate false diffracted waves. To eliminate grid diffraction, many scholars have carried out extensive research, including using variable grid algorithms [8], vertical grid mapping based on coordinate transformation [9], body fitted grids [10], and other methods. These methods effectively suppress stepped diffraction to a certain extent, but they also have problems. For example, the variable mesh algorithm essentially uses a smaller mesh step size to discretize the undulating interface, so it cannot fundamentally eliminate the stepped diffraction; the difference method based on coordinate transformation requires a corresponding transformation of the wave equation; and the implementation process is relatively complex. In recent years, the meshless method, which avoids mesh generation, discretizes the solution area into a series of independent nodes and constructs approximate functions on these discrete nodes to obtain linear equations [11]. Since the nodes change flexibly with the velocity model, the meshless method can effectively avoid grid diffraction [12]. The generalized finite difference method is a meshless method with a simple principle, flexible node arrangement, and high calculation accuracy, and it has been widely used to solve a variety of mathematical and engineering problems [13; 14; 15]. [PERSON] (1972) first used the Taylor series expansion of several adjacent nodes around a center point, based on a distance function, to solve differential equations [16]. Later, [PERSON] et al. (2001) developed the explicit formula of the generalized finite difference method [17]. Based on the Taylor series expansion and weighted least square fitting, the partial derivative of the unknown quantity in the differential equation is expressed as a linear combination of the function values of several adjacent nodes. Thus, the basic theory of the generalized finite difference method was formed. Since then, many scholars have improved the theory of this method and its implementation technology. For example, [PERSON] and [PERSON] systematically analyzed various factors that affect the accuracy of the GFDM calculations, such as node generation [18; 19], star of nodes shape [20], and weight function [17], and they pointed out that as the discretization (cloud of nodes) becomes more regular, the results obtained become more stable [21]. [PERSON] et al. (2011) first applied the generalized finite difference method to the forward modeling of elastic wave equations, and proposed an adaptive method to minimize the effect of the irregularity of node distributions on dispersion [22]; on this basis, [PERSON] et al. (2012) discussed the dispersion and stability of elastic wave forward modeling using the generalized finite difference method under regular and irregular node conditions [23]. [PERSON] et al. (2013) further studied the generalized finite difference method for solving the problem of determining seismic wave propagation in homogeneous isotropic media [24]. [PERSON] et al.