## 1 Introduction  
Over the past three decades, a variety of surface-wave exploration methods has become an increasingly important means for inferring the properties of subsurface structures, especially in the near-surface region [1; 2]. Compared with the commonly employed seismic body waves, the seismic surface waves (SWs), which primarily include Rayleigh and Love waves in land acquisition, usually propagate in the vicinity of the free surface and exhibit conspicuous dispersion characteristics in seismic records. This dispersion characteristic of SWs is generally delineated by the so-called phase- and/or group-velocity dispersion curves (DCs). The near-surface shear-wave velocities (Vs) can be obtained by inverting these DCs. The inverse problem of DCs, inherently nonlinear, can be solved by using global optimization methods, for instance, genetic algorithm [3], simulated annealing [4], particle swarm optimization [5], etc., or by using various linearized inversion strategies [1; 6; 7]. The success of either DCs inversion strategy is based on the accurate and efficient computation of surface-wave DCs. In addition, if the latter strategies are chosen, reliable calculations of partial derivatives (PDs) of DCs with respect to the parameters of subsurface strata are indispensable to create the Jacobian matrix (or the sensitivity matrix) in this class of strategies, especially due to the demand of their repeated calculations during the entire inversion process. The efficient computation of surface-wave DCs benefits not only from the pioneering groundwork of [PERSON] and [PERSON] [8; 9], but also from the rapid development of computer techniques after the 1960s. Inspired by [PERSON] and [PERSON], numerous approaches have been proposed since then [10; 11; 12; 13; 14; 15; 16; 17; 18; 19; 20; 21; 22; 23], particularly to remedy the issue of loss of numerical precision at high frequencies, which is an inherent shortcoming of [PERSON]'s theory (see reference [23] for more details). Among all these approaches, the one put forward by [PERSON] [13] and its improved variants [14; 23] are frequently utilized as the foundation of various surface-wave DCs inversion strategies due to their simplicity, efficiency, and effectiveness. [PERSON]'s approach is based on the so-called delta matrix theory [13], with which the shortcoming of [PERSON]'s theory are thoroughly overcome. [PERSON] et al. [10] developed the first computer program for computation of Rayleigh- and Love-wave DCs. In spite of the high efficiency in computing phase-velocity DCs, they applied the numerical differentiation of phase velocity to acquire the group velocity, which is a brute-force, low-precision, and unreliable method. The first category of analytical approaches for calculating the surface-wave group-velocity DCs and their PDs is based on the variational theory [24; 25; 26; 27; 28; 29]. By virtue of [PERSON]'s principle, [PERSON] [24] first derived expressions for computing the group velocities of surface waves analytically; nonetheless, the PDs of DCs were approximated to the first order accuracy with numerical differentiation. Thereafter, the same idea was adapted to tackle more sophisticated cases, such as [PERSON] et al. [25; 26], [PERSON] [27], [PERSON] [28], and [PERSON] et al. [29], just to name a few. On the other hand, to bypass the theoretical complexity of these approaches, [PERSON] introduced the implicit function method, another analytical approach, first for Love waves [30; 31] and then for Rayleigh waves in the case of a layer overlying a half-space [32]. Then [PERSON] et al. [33] succeeded in generalizing this approach to the case of multiple layers for Rayleigh waves. With the help of [PERSON]'s algorithm [11; 15; 16], in combination with the implicit function method, they successfully acquired the analytical expressions for computing the PDs of Rayleigh-wave phase velocities, and the group velocities are also acquired by further using [PERSON]'s method [34].