## 3 Methodology
### Extraction Method of First-Arrival Information  
In VSP, first-arrival waves, which are the earliest wave motions that reach receivers, are an important component in FWIs. First arrivals typically represent the initial wave propagations along the shortest paths through subsurface media, and their arrival times and propagation characteristics provide crucial information about velocity structures. An effective extraction and utilization of first-arrival information can significantly enhance the stability and accuracy of full-waveform inversion. In the present study, we extracted first-arrival time differences by computing cross-correlation functions between the observed and simulated data and used these differences to select appropriate seismic traces for inversion.  
Consider a single seismic trace, where \(\mathbf{d}^{s}_{\text{obs}}(t)\) and \(\mathbf{d}^{s}_{\text{cal}}(t)\) represent the observed and simulated data, s denotes the trace index, and t denotes time. Each trace consists of a time series. To extract the temporal difference in first arrivals between observed and simulated data, we computed their cross-correlation function \(\mathcal{C}^{s}(\tau)\) as follows:  
\[\mathcal{C}^{s}(\tau)=\int\mathbf{d}^{s}_{\text{obs}}(t)\mathbf{d}^{s}_{\text {cal}}(t+\tau)dt, \tag{12}\]  
where \(\tau\) denotes the time lag. We designated the peak position \(\tau^{s}_{max}\) the cross-correlation function as the optimal temporal alignment between the observed and synthetic data for this trace, expressed as follows:  
\[\tau^{s}_{\text{max}}=\text{argmax}_{\tau}\mathcal{C}^{s}(\tau), \tag{13}\]  
The peak position \(\tau^{s}_{max}\) indicates the temporal discrepancy in first-arrival times between simulated and observed data. The magnitude of \(|\tau^{s}_{\text{max}}|\) increases proportionally with the disparity between the initial velocity model and the true model. Including traces with substantial discrepancies in loss computation and gradient, backpropagation can potentially lead to erroneous parameter update directions. To quantify this temporal discrepancy, we defined a metric by subtracting the total finite difference simulation time T from \(\tau^{s}_{\text{max}}\), yielding the following:  
\[\Delta t^{s}=\tau^{s}_{\text{max}}-T, \tag{14}\]where \(\Delta t^{\rm s}\) denotes the temporal phase differential between the observed and simulated data for the s-th seismic trace. This first-arrival phase differential serves as the criterion for dynamic trace selection in the loss computation. We established a temporal threshold \(\Delta t_{\rm threshold}\), where traces satisfying \(|\Delta t^{\rm s}|<\Delta t_{\rm threshold}\) exhibited sufficient coherence between their embedded subsurface velocity information and the current initial model; conversely, traces exceeding this threshold were temporarily omitted from the loss computation. The mathematical formulation is as follows:  
\[\mathcal{S}_{\rm selected}=\{s\in\mathcal{S}|\ |\Delta t^{\rm s}|<\Delta t_{\rm threshold }\}, \tag{15}\]  
where \(\mathcal{S}\) denotes the complete ensemble of seismic traces and \(\mathcal{S}_{\rm selected}\) represents the subset of selected traces. This dynamic selection mechanism ensures that during the initial stages of inversion, when the disparity between the initial and true models is substantial, the inversion process exclusively incorporates traces that exhibit minimal first-arrival phase differentials between the observed and synthetic data. This methodology effectively mitigates the propagation of errors stemming from initial velocity model inaccuracies and enhances the inversion stability. As the inversion evolves, additional seismic traces are progressively incorporated into the inversion process. This methodological approach initiates inversion with high-fidelity data and systematically expands the dataset, thereby circumventing local minima that can arise from data insufficiency and optimizing the comprehensiveness of the inversion results.