Data processing method and system for micro-motion detection array station distribution

Abstract

The invention discloses a data processing method and system for micro-motion detection array station distribution, and relates to the technical field of micro-motion detection, and the method comprises the steps: according to the actually measured starting point coordinates and the azimuth angle of a measuring line, setting the measuring point distance, and the radius and number of nested triangle circumcircles, generating field station coordinates in batches; collecting and preprocessing station waveform data; calculating a spatial autocorrelation coefficient of any center station and a non-center station; dBSCAN clustering is carried out based on the station distance; performing spatial averaging on the spatial autocorrelation coefficient according to a clustering result; fitting the averaged spatial autocorrelation coefficient and a zero-order Bessel function to obtain a frequency dispersion spectrum, picking up a surface wave phase velocity in the frequency dispersion spectrum, and calculating an apparent shear wave velocity; inverting the shear wave velocity of the waveform data by adopting a self-adaptive gradient inversion strategy; and generating a velocity comprehensive profile of a depth domain and a frequency domain according to the surface wave phase velocity, the apparent shear wave velocity and the inverted shear wave velocity. The invention provides a new method for station distribution to data processing for micro-motion detection.

Description

Technical Field

[0001] The present invention relates to the technical field of micro - motion detection, and particularly to a data processing method and system for array station layout in micro - motion detection. Background Technique

[0002] Micro - motion detection is developed on the basis of background noise imaging research. It uses the weak vibrations existing on the earth's surface as the signal source. This method is not affected by strong electromagnetic interference, industrial activity interference, and does not require an active seismic source. Currently, micro - motion exploration is mostly used in engineering geological exploration and urban geological surveys, and some exploratory work in mineral exploration can also be seen. Although the micro - motion exploration technology has shown many advantages in geological detection in urban and complex environments, there are still several major technical or algorithm problems to be solved urgently: (1) Currently, for the deployment of micro - motion detection stations, the common method is to manually calculate coordinate points on - site using RTK and place acquisition instruments according to the calculated coordinate points. This operation method is inefficient and prone to misplacement and omission of acquisition instruments. There is still a lack of a set of pre - observation system design methods; (2) There is a phenomenon of segmentation between the scientific research and engineering fields in micro - motion detection, and the data formats supported by different calculation methods are single. For example, the common calculation methods in the scientific research field support the SAC data format, while the engineering field methods mostly support the sg2 data format. The lack of compatibility between the two limits the flexibility and wide application of data processing; (3) The vulnerability of the SPAC algorithm. Currently, the technology defaults to performing spatial autocorrelation (SPAC) calculations centered on the central station. Once the equipment of the central station fails, is lost, or is damaged, it may lead to the failure of the entire acquisition work and the inability to carry out standard - type SPAC method calculations; (4) The existing technology does not introduce advanced clustering analysis ideas, and it is impossible to quality - control the classification of stations. Engineering and technical personnel cannot intelligently screen stations at different distances to participate in the calculation, resulting in the data of stations with large interference contaminating the dispersion energy spectrum; (5) The existing inversion algorithm generates a series of temporary model parameters by traversing different regularization parameters and search steps, calculates the corresponding mismatch function, and selects the model parameters corresponding to its minimum value as the update result. The inversion strategy has a slow convergence speed, is prone to falling into local optima, and the step - size selection depends on preset parameters, and it is impossible to balance the convergence speed and the stability of the algorithm well. The initial model setting, inversion parameter setting, etc. are more dependent on experience. If the parameter setting is improper, the iteration is prone to falling into local minima, affecting the inversion effect. (6) The drawing of the comprehensive profile is complex and cumbersome, and there is no convenient and fast method for generating and drawing the comprehensive profile of phase velocity, apparent shear - wave velocity, and inverted shear - wave velocity. Summary of the Invention

[0003] The purpose of the present invention is to propose a data processing method for array station layout in micro - motion detection to solve the problems of insufficient adaptability and stability of the micro - motion detection data processing algorithm in the existing geological exploration field, including the following steps: S1. Set the measuring point distance, the radius of the circumcircle of the nested triangle, and the number of nested triangles according to the starting point coordinates measured by RTK and the azimuth angle of the survey line, and generate the coordinates of all observation stations; S2. Collect the waveform data of the observation stations corresponding to the coordinates, and preprocess the waveform data; S3. Calculate the spatial autocorrelation coefficient of the preprocessed waveform data between any specified central station and non-central stations; S4. Perform DBSCAN clustering on the pairs of observation stations according to the distance between the pairs of stations to obtain the clustering result; S5. Perform spatial averaging on the spatial autocorrelation coefficients of the pairs of observation stations according to the clustering result, fit the averaged spatial autocorrelation coefficient with the zero-order Bessel function to obtain the dispersion energy spectrum, pick up the phase velocity dispersion data in the energy spectrum, and calculate the apparent shear wave velocity dispersion curve; S6. Invert the picked-up phase velocity dispersion data based on the adaptive gradient inversion strategy to obtain the shear wave velocity; S7. Generate a comprehensive profile in the depth domain and a comprehensive profile in the frequency domain according to the picked-up surface wave phase velocity, the calculated apparent shear wave velocity, and the inverted shear wave velocity, and load the station names or topography on the profiles to guide the field exploration work.

[0004] Further, S1 is specifically as follows: Generate a series of measuring point coordinates according to the given starting point coordinates, azimuth angle, measuring point distance, and quantity, and draw nested triangles with a specified data quantity at each measuring point. The generation of the nested triangles is based on the geometric characteristics of the circumcircle radius and equilateral triangles, and multiple layers of nesting are realized by recursive calculation of the midpoints. The center coordinates of all circumcircles and the vertex coordinates of the nested triangles are the coordinates of the observation stations; When the number of nested triangles is set to 0, the observation stations are linearly arranged; when the number of nested triangles is set to 1 - 4, the observation stations are in a two-dimensional array.

[0005] Further, the preprocessing methods include downsampling, de-meaning, detrending, maximum normalization or moving absolute average normalization, spectral whitening, and band-pass filtering.

[0006] Further, S4 is specifically as follows: Calculate the distances between pairs of stations, form a distance list, and convert the distance list into an array. According to the allowable deviation percentage dev parameter, calculate the neighborhood radius eps of DBSCAN: eps = (dev / 100) × mean_distance, where mean_distance is the mean of the distance list; extract the indices and distance values of the pairs of stations in each cluster. The clustering result includes: station names, indices of pairs of stations, and distance list.

[0007] Further, the spatially autocorrelated coefficients after fitting and averaging are used to obtain the dispersion energy spectrum with the zero-order Bessel function. The phase velocity dispersion data in the energy spectrum are picked up, specifically as follows: Through velocity scanning, the spatially autocorrelated coefficients at different station spacings after clustering at the same frequency are fitted with the first-kind zero-order Bessel function to generate a dispersion energy spectrum with frequency varying with velocity. Under the interface of the dispersion energy spectrum, press the left mouse button to trigger the picking function, obtain the frequency and velocity at the current mouse position to analyze the dispersion energy spectrum, find the velocity point with the lowest energy at each frequency point. During the process of dragging the mouse, trigger the picking function. For each new frequency point, check whether it is necessary to continue picking the velocity point with the lowest energy at the frequency point. Release the left mouse button to complete the range selection and finish picking the dispersion curve.

[0008] Further, according to the phase velocity dispersion data in the picked-up dispersion energy spectrum, the apparent shear wave velocity calculation formula is as follows: Apparent shear wave velocity =

[0009] where, represents the frequency of the first dispersion point, represents the frequency of the second dispersion point, represents the phase velocity of the first dispersion point, represents the phase velocity of the second dispersion point.

[0010] Further, an adaptive gradient inversion strategy is adopted for inversion, introducing a momentum term, an adaptive parameter update strategy, and a line search mechanism. The momentum term is updated by the momentum decay factor:

[0011] where, V represents the momentum term, is the search direction vector, represents the decay factor; The decay factor decays according to the following method based on the number of iterations:

[0012] The step size that satisfies the Armijo condition is selected through the line search mechanism to update the model;

[0013] In the formula, is the search step size, is the updated formation model, represents the current formation model, is the parameter in Armijo, is the search gradient direction; If the Armijo criterion is not satisfied, the current step size is adjusted by the bisection method until the condition is met; Perform adaptive adjustment on: where:

[0014]

[0015] Among them, is the adjustment rate, which is used for smooth adjustment, represents the standard deviation of the function value sequence, k represents the number of loops, and maxIter represents the maximum number of iterations.

[0016] Furthermore, according to the picked surface wave phase velocity, the calculated apparent shear wave velocity, and the inverted shear wave velocity, a comprehensive profile in the depth domain and a comprehensive profile in the frequency domain are generated. Specifically: Using the Kriging interpolation method, based on the coordinate data and velocity data of the stations, a two-dimensional velocity profile is generated. The y-axis is set as the frequency to obtain the profile in the frequency domain, and the y-axis is set as the depth to obtain the profile in the depth domain. Elevation static correction is performed on the profile in the depth domain according to the elevation data.

[0017] The present invention also proposes a data processing system for the layout of a microtremor detection array, including: An observation station generation unit, which is used to set the measurement point distance, the radius of the circumcircle of the nested triangle, and the number of nested triangles according to the starting point coordinates measured by RTK and the azimuth of the survey line, and generate the coordinates of all observation stations; A data preprocessing unit, which is used to collect the waveform data of the observation stations corresponding to the coordinates and preprocess the waveform data; A spatial autocorrelation coefficient calculation unit, which is used to pair the observation stations in pairs and calculate the spatial autocorrelation coefficient of the pairs of observation stations according to the waveform data of the observation stations; A clustering analysis unit, which is used to perform DBSCAN clustering on the pairs of observation stations according to the distance between the pairs of stations to obtain the clustering result; A dispersion point generation unit, which is used to perform spatial averaging on the spatial autocorrelation coefficients of the pairs of observation stations according to the clustering result, fit the averaged spatial autocorrelation coefficient with the zero-order Bessel function to obtain the dispersion energy spectrum, pick up the phase velocity dispersion data in the energy spectrum, and calculate the apparent shear wave velocity dispersion curve; An inversion unit, which is used to invert the picked phase velocity dispersion data based on the adaptive gradient inversion strategy to obtain the shear wave velocity; A profile generation unit, which is used to generate a comprehensive profile in the depth domain and a comprehensive profile in the frequency domain according to the picked surface wave phase velocity, the calculated apparent shear wave velocity, and the inverted shear wave velocity, and load the station name or terrain on the profile to guide the field exploration work.

[0018] The beneficial effects brought by the technical solution provided by the present invention are as follows: Based on the measured coordinates of the measuring points and the azimuth of the measuring line, the present invention generates the coordinates of all observation stations according to the nested triangle or linear observation system; performs DBSCAN clustering on the observation stations, and according to the clustering results, averages the spatial autocorrelation coefficients of the same clustering clusters to improve the calculation effect of the dispersion energy spectrum; batch interpolates and generates a comprehensive profile in the depth domain and a comprehensive profile in the frequency domain based on the picked surface wave phase velocity, the calculated apparent shear wave velocity, and the inverted shear wave velocity, and loads the station names or topography on the profile to guide the field exploration work. The present invention integrates the entire process from station deployment to data processing, greatly improving the field construction efficiency of microtremor exploration in fields such as engineering geological exploration and mineral survey, improving the accuracy and precision of station deployment, and improving the overall process matching degree and operability from data acquisition, data processing, data inversion to result output. Description of the Drawings

[0019] Figure 1 is a flowchart of the data processing method for the microtremor detection array station layout in the embodiment of the present invention; Figure 2 is the acquisition space distribution map of the stations generated in the embodiment of the present invention; Figure 3 is the waveform data quality control chart generated in the embodiment of the present invention, Figure 3 in which (a) is the 30-minute microtremor raw waveform data collected according to the observation system, Figure 3 in which (b) is the seismic waveform data preprocessed according to the preprocessing method of the present invention; Figure 4 is the spatial autocorrelation curve graph generated in the embodiment of the present invention; Figure 4 in which (a) is the spatially averaged spatial autocorrelation curve graph when the station spacing is 15 m, Figure 4 in which (b) is the spatially averaged spatial autocorrelation curve graph when the station spacing is 30 m; Figure 5 is the dispersion energy spectrum and velocity curve obtained by the integration of manual and automatic methods with the 73rd station as the central station in the embodiment of the present invention, where, Figure 5 in which (a) is the dispersion energy spectrum calculated with the 73rd station as the central station, Figure 5 in which (b) is the phase velocity corresponding to the 73rd station as the central station and the calculated apparent shear wave velocity curve; Figure 6 is the dispersion energy spectrum and velocity curve obtained by the integration of manual and automatic methods with the 72nd non-central station in the embodiment of the present invention, where, Figure 6 in which (a) is the dispersion energy spectrum calculated with the 72nd non-central station, Figure 6 in which (b) is the phase velocity corresponding to the 72nd non-central station and the calculated apparent shear wave velocity curve; Figure 7 These are the dispersion fitting and model inversion results of the embodiments of the present invention. Among them, Figure 7 in (a), it is the comparison between the observed and predicted values of frequency and phase velocity; Figure 7 in (b), it is the initial value and inversion result of shear wave velocity versus depth; Figure 8 These are all the dispersion data points extracted from the stations in the embodiments of the present invention; Figure 9 This is the phase velocity profile in the frequency domain of the embodiments of the present invention; Figure 10 This is the apparent shear wave velocity profile in the frequency domain of the embodiments of the present invention; Figure 11 This is the shear wave velocity profile inversed in the frequency domain of the embodiments of the present invention; Figure 12 This is the phase velocity profile in the depth domain of the embodiments of the present invention; Figure 13 This is the apparent shear wave velocity profile in the depth domain of the embodiments of the present invention; Figure 14 This is the shear wave velocity profile inversed in the depth domain of the embodiments of the present invention. Specific embodiments

[0020] To make the objectives, technical solutions, and advantages of the present invention clearer, the embodiments of the present invention will be further described below in conjunction with the accompanying drawings.

[0021] The flowchart of the data processing method for the microtremor detection array station layout in the embodiments of the present invention is as shown in Figure 1 , and specifically includes the following steps: S1. According to the starting point coordinates measured by RTK (real-time kinematic positioning) and the azimuth angle of the survey line, set the measuring point distance, the circumradius of the nested triangle, and the number of nested triangles to generate the coordinates of all observation stations.

[0022] The setting of the circumradius of the nested triangle array is related to the exploration depth. Generally, the empirical relationship is that the exploration depth is 3 to 5 times the circumradius.

[0023] Initialize the parameters: the starting point coordinates (start_x, start_y), that is, the east coordinate and north coordinate of the measurement starting point, the azimuth angle of the survey line (angle), that is, the angle of the survey line relative to the positive direction of the X-axis, with the counterclockwise direction being positive, the measuring point distance (distance), that is, the distance between two adjacent measuring points, the number of measuring points (count), that is, the total number of generated measuring points, the outer circle radius (initial_radius), that is, the circumradius of the outer equilateral triangle, and the number of nested triangles (num_nested), that is, the number of layers of the drawn nested equilateral triangles.

[0024] Measurement point generation: Based on the starting point coordinates, azimuth angle, measurement point distance, and the number of measurement points, a series of measurement point coordinates evenly distributed along the specified direction are generated. Convert the azimuth angle angle from degrees to radians using the formula:

[0025] For the i-th measurement point, its east coordinate x and north coordinate y are calculated respectively using the following formulas:

[0026]

[0027] Triangle generation: Generate the coordinates of the three vertices of the initial equilateral triangle based on the circumradius, and generate nested triangle levels as needed. The calculation idea is that the angle between the three vertices of an equilateral triangle is 120 degrees each. By adjusting the angle_offset, the overall direction of the triangle can be rotated:

[0028]

[0029] , where i is 0, 1, 2.

[0030] For an equilateral triangle, the relationship between the circumradius R and the side length is: , where a is the distance between any two vertices of the triangle.

[0031] Two points 、 The midpoint calculation formula is: , which is used to generate the three vertices of the next layer of nested triangles.

[0032] Drawing nested triangles: Based on the coordinates generated by the measurement points, draw nested equilateral triangles and perform numbering and annotation. The calculation idea is: Based on the coordinates of each measurement point and the circumradius, generate the vertices of the initial equilateral triangle. The vertices of each layer of nested triangles are obtained by calculating the midpoints of the sides of the previous layer of triangles. Recursively generate num_nested layers of nested triangles, and perform numbering and annotation on each measurement point and the vertices of its nested triangles.

[0033] Coordinate point data export: Write the east and north coordinates of all measurement points and triangle vertices row by row to a CSV file for subsequent field measurements.

[0034] The centers of all circumcircles and the vertex coordinates of the nested triangles are the coordinates of the observation stations. And number each station.

[0035] By exporting the calculated coordinates, the rationality of the detection calculation on the Ovi map can be projected and post - correction can be carried out. When the array is linearly arranged, the SPAC or ESPAC algorithm can be calculated according to the characteristics of the linear arrangement. When the array is nested triangularly arranged, the SPAC or ESPAC algorithm can be calculated according to the characteristics of the nested arrangement. At the same time, for the nested triangular stations, a certain station can be selected as the center point, which greatly improves the rule that the previous processing technology only defaults the station with the smallest number as the center point, and also avoids the limitation that some processing technologies artificially restrict the stations to be deployed clockwise or counter - clockwise.

[0036] S2. Collect the waveform data of the observation stations corresponding to the coordinates, read the micro - motion data in multiple formats (sac, sg2), and pre - process the waveform data.

[0037] The present invention supports waveform data in sac and sg2 formats, writes the station numbers into the sac file header, solves the problem of missing station information in station data, and by matching the station numbers in the coordinate file with the data with station numbers in the waveform data, correct data processing can be carried out later.

[0038] When collecting the waveform data of the observation stations, the observation time is long and the data volume is large. When many systems directly import the original data, the computing memory is often limited and it is difficult to process. Down - sampling the original waveform to a reasonable sampling rate range, the inner layer of a single waveform file meets the requirements of computing memory and computing speed, greatly improving the work efficiency.

[0039] The present invention designs a variety of optional pre - processing methods, including down - sampling (optional), de - meaning, de - trending, maximum normalization or moving absolute average trace normalization, spectral whitening, band - pass filtering. Any one of the processing methods can be manually selected. By comparing the changes in the waveform data after pre - processing with the original waveform data, the processing method that meets the geological background and acquisition background can be selected. The automated and anti - aliasing batch down - sampling method can efficiently process a large number of SAC - format seismic data files. By combining low - pass filtering and anti - aliasing, the signal quality after down - sampling is ensured.

[0040] The implementation process of the pre - processing method is as follows: Down - sampling is to reduce the sampling rate of the digital signal and reduce the number of sampling points per unit time. For example, if the original signal sampling rate is and the down - sampling factor is M , then the new sampling rate after down - sampling is which is:

[0041] To prevent aliasing, the cut - off frequency of the low - pass filter is set to half of the new Nyquist frequency:

[0042] The calculation process is as follows: batch read SAC files to obtain waveform data and the original sampling rate; set a new decimation factor according to the numerical characteristics of the original sampling rate, and perform decimation calculation using the anti-aliasing low-pass filter method; batch output the decimated SAC format data.

[0043] The purpose of de-meaning is to eliminate the DC component in the signal and make the average value of the signal zero, which is achieved by calculating the mean value of the signal and subtracting it from each sampling point. For example, if the input signal is , and its length is N , the de-meaned signal is The calculation is as follows: , ,

[0044] The purpose of detrending is to eliminate the linear trend in the signal, make the signal more stable and facilitate subsequent processing, which is achieved by fitting the trend line of the signal and subtracting it from the signal. For example, if the de-meaned signal is , the fitted linear trend is , where and are determined by the least squares method for the slope and intercept. The detrended signal is calculated as follows: ,

[0045] Moving absolute average normalization is to achieve local normalization of the signal by calculating the moving average of the absolute value of the signal. For example, if the input signal is the detrended signal , the moving window size is , and the moving average is calculated as follows: , and the normalized signal , .

[0046] Maximum value normalization scales the amplitude range of the signal to between [-1, 1] by processing each sampling point of the signal with the maximum value of the absolute value of the signal. The detrended signal is , and the maximum value normalization can be expressed as: ,

[0047] Spectral whitening aims to flatten the spectrum of a signal, making the energy distribution of each frequency component more uniform. The process is to perform a Fourier transform on the input signal to obtain frequency-domain data, calculate the amplitude or power spectrum of each frequency component, normalize the amplitude of the frequency-domain signal, and then perform an inverse Fourier transform to obtain the whitened signal in the time domain. If the input signal is a normalized signal , its Fourier transform is , and the frequency-domain signal after spectral whitening is calculated as: . Performing an inverse Fourier transform on it gives the time-domain data after spectral whitening: .

[0048] S3. Set any station as the central station, pair up the N observation stations in pairs to obtain (N - 1) pairs of station pairs, and calculate the spatial autocorrelation coefficient (SPAC) of the station pairs based on the waveform data of the observation stations. The mainstream processing technology defaults the station at the center point as the central station and then calculates using the SPAC method. If there are accidents such as equipment failure, loss, or damage at the central station, it will not be possible to meet the principle of the SPAC algorithm to carry out the work, which may lead to the invalidation of the entire acquisition work. In the present invention, the central station can be arbitrarily designated, and calculating SPAC is more in line with the actual production scenario.

[0049] S4. Perform distance clustering on all observation station pairs. In data processing, the distances between station pairs are not the same. To effectively classify station pairs with similar distances for subsequent dispersion spectrum calculation and analysis, the present invention uses the DBSCAN (Density-Based Spatial Clustering of Applications with Noise) density-based clustering analysis algorithm to group points in high-density regions into one cluster and regard points in low-density regions as noise. Different from the traditional clustering analysis method based on mean deviation, DBSCAN automatically identifies the number of clusters through the concept of density and can effectively handle noise points, making the clustering process more efficient and robust.

[0050] The implementation process includes: calculating the distances between all stations to form a distance list dist_list, converting the distance list into a Numpy array, and converting the one-dimensional data into two-dimensional data to meet the input requirements of the DBSCAN algorithm; calculating the neighborhood radius eps of DBSCAN according to the allowed deviation percentage dev parameter, that is, eps = (dev / 100) × mean_distance, where mean_distance is the mean of the distance list. Set the parameters eps and min_samples = 1, where min_samples = 1 allows each cluster to contain at least one sample point.

[0051] Using the OrderedDict class in the Python library, sort the clustering results based on the mean of each cluster to ensure that the output results are arranged in ascending order of the mean; traverse all unique cluster labels, extract the indices and actual distance values of the elements in each cluster, and calculate the mean of the cluster; store the list of cluster indices and the list of distance values in idx_list and all_list respectively. Through the DBSCAN clustering method, the clustering results not only contain the distance list and the indices of the station pairs, but also can be associated with the specific station names, facilitating the subsequent calculation of the spatial autocorrelation coefficient and the dispersion energy spectrum.

[0052] S5. Perform spatial averaging on the spatial autocorrelation coefficients of the observed station pairs according to the clustering results, fit the averaged spatial autocorrelation coefficients with the zero-order Bessel function to obtain the dispersion energy spectrum, pick up the phase velocity dispersion data in the energy spectrum, and calculate the apparent shear wave velocity dispersion curve.

[0053] The process of performing a dispersion energy spectrum scan and picking up the dispersion curve using manual and semi-automatic methods includes: Through velocity scanning, fit the spatial autocorrelation coefficients of different station spacings after clustering at the same frequency with the first-kind zero-order Bessel function to generate a dispersion energy spectrum of frequency varying with velocity. The formula for the dispersion energy spectrum is:

[0054] where is the fitting error, is the spatial autocorrelation coefficient, is the first-kind zero-order Bessel function, , is the spatial distance of the j-th station. Under the dispersion energy spectrum interface, when the user presses the left mouse button, the picking function is triggered to obtain the frequency and velocity analysis of the dispersion energy spectrum at the current mouse position. Near each frequency point, the velocity point with the lowest energy is automatically searched for. If a valid point is found, it is added to the list, and the data and graph are updated. During the mouse dragging process, the picking function is triggered to check whether it is necessary to continue picking the velocity point with the lowest energy (based on the frequency interval threshold). For each new frequency point, the "find the lowest energy value" function is called to automatically find the lowest energy point and add it to the list, and the data and graph are updated. When the user releases the left mouse button, all picked dispersion points are automatically saved to a file. This method reduces manual operations, improves efficiency and accuracy, and solves the problem of the need for frequent manual and automatic switching operations in the past. The present invention still retains the function of manually adding and deleting points, providing flexibility and controllability. At the same time, the sliding window technology is introduced during the picking process to perform mean smoothing processing on the picked dispersion points, reduce the influence of noise, and improve the continuity and reliability of the data; based on the automatically picked dispersion points, using the relationship between velocity and frequency, the apparent shear wave velocity is automatically calculated to ensure the physical rationality of the calculation results.

[0055] Apparent shear wave velocity =

[0056] where represents the frequency of the first dispersion point, represents the frequency of the second dispersion point, represents the phase velocity of the first dispersion point, represents the phase velocity of the second dispersion point.

[0057] S6. Use the adaptive gradient algorithm strategy to invert the phase velocity dispersion data picked in the previous step to obtain the shear wave velocity.

[0058] According to the adaptive gradient algorithm strategy, a momentum term, an adaptive parameter update strategy, and a line search mechanism are introduced. The search direction is mainly determined by the momentum decay factor to update the momentum term. The calculation formula is:

[0059] In the formula, is the momentum term, is the search direction vector. On the premise of improving the search process, the decay factor decays according to the number of iterations as follows:

[0060] After determining the search direction, the step size that satisfies the Armijo criterion condition is selected through the line search mechanism to update the model, so as to achieve the purpose of fast and stable inversion.

[0061]

[0062] Wherein, is the updated formation model, is the parameter in the Armijo criterion, is the search gradient vector. is the search step size. If the Armijo criterion is not satisfied, the current step size is adjusted by the bisection method until the condition is met.

[0063] In order to balance the convergence speed and algorithm stability, the present invention provides an adaptive adjustment based on the change of the objective function . Mainly according to the standard deviation of the function value sequence to measure the volatility of the current optimization process.

[0064] If the volatility of the objective function is large, the growth of the step size can be suppressed by reducing . That is:

[0065] Wherein, is the adjustment rate, is used for smooth adjustment to ensure that will not decrease too low or too fast.

[0066] This method sets in the following form, aiming to have a fast adjustment speed of the step size in the initial stage of optimization to avoid premature convergence, and in the later stage, the step size adjustment becomes more delicate.

[0067]

[0068] where k represents the number of loops and maxIter represents the maximum number of iterations.

[0069] S7. Generate a comprehensive profile in the depth domain and a comprehensive profile in the frequency domain according to the picked surface wave phase velocity, the calculated apparent shear wave velocity, and the inverted shear wave velocity. The optional profile can load the station name or terrain, which is used to guide the field exploration work.

[0070] Read the station coordinate file, including the station name, east coordinate, north coordinate, and elevation coordinate; read each velocity file, and the velocity file includes: the picked surface wave phase velocity data, the calculated apparent shear wave velocity data, and the inverted shear wave velocity data; Using the Kriging interpolation method, generate a two-dimensional velocity profile according to the geographical distribution of the stations and the velocity data. Set the y-axis as the frequency to obtain the profile in the frequency domain, set the y-axis as the depth to obtain the profile in the depth domain, and perform elevation static correction on the profile in the depth domain.

[0071] In an exemplary embodiment, a data processing system including a microseismic detection array station layout includes: An observation station generation unit, configured to set the measurement point distance, the circumradius of the nested triangle, and the number of nested triangles according to the starting point coordinates measured by RTK and the azimuth angle of the survey line, and generate the coordinates of all observation stations; A data preprocessing unit, configured to collect the waveform data of the observation stations corresponding to the coordinates and preprocess the waveform data; A spatial autocorrelation coefficient calculation unit, configured to pair the observation stations in pairs and calculate the spatial autocorrelation coefficient of the pair of observation stations according to the waveform data of the observation stations; A clustering analysis unit, configured to perform DBSCAN clustering on the pairs of observation stations according to the distances between the pairs of stations to obtain a clustering result; A dispersion point generation unit, configured to perform spatial averaging on the spatial autocorrelation coefficients of the pairs of observation stations according to the clustering result, fit the averaged spatial autocorrelation coefficient with the zero-order Bessel function to obtain a dispersion energy spectrum, pick up the phase velocity dispersion data in the energy spectrum, and calculate the apparent shear wave velocity dispersion curve; An inversion unit, configured to invert the picked-up phase velocity dispersion data based on an adaptive gradient inversion strategy to obtain the shear wave velocity; A profile generation unit, configured to generate a comprehensive profile in the depth domain and a comprehensive profile in the frequency domain according to the picked-up surface wave phase velocity, the calculated apparent shear wave velocity, and the inverted shear wave velocity, and load the station names or topography on the profile to guide the field exploration work.

[0072] Taking a section of field-measured microseismic data as an example, 7 three-component ground pulsation meters were deployed during the work (the station serial numbers of each measurement point were 37, 65, 68, 70, 71, 72, 73), and the observation system was designed as a double-nested triangle, where the 73rd station was designed as the center point. By collecting the starting point coordinates of the measurement points x = 764456.01, y = 2513553.63, the survey line azimuth angle -45° (positive counterclockwise along the X-axis, negative clockwise), the outer circle radius 30m, and the measurement point spacing 10m, the acquisition spatial distribution map of the stations can be generated at one time using the method of the present invention. Refer to Figure 2 , Figure 2 is the acquisition spatial distribution map of the stations generated by the embodiment of the present invention, Figure 2 only shows 15 measurement points and all the coordinates of the observation stations ( Figure 2 shows some of the center point coordinate values). Through the station array layout method of the present invention, the time for field station layout is greatly reduced, and the work efficiency is effectively improved.

[0073] Figure 3 is the waveform data quality control chart generated by the embodiment of the present invention, Figure 3In (a), it is the original microtremor waveform data of 30 minutes collected according to the observation system. Figure 3 In (b), it is the seismic waveform data preprocessed according to the preprocessing method of the present invention. Figure 3 In (b), it is the seismic waveform data after mean removal, trend removal, moving absolute smoothing normalization (window of 2001 sampling points), spectral whitening, and band-pass filtering according to the system method of the present invention. It can be seen that the waveform data after preprocessing by the system of the present invention is cleaner and conforms to the waveform trend of stationary randomness.

[0074] Figure 4 It is according to Figure 3 The spatial autocorrelation curve calculated from the preprocessed microtremor waveform data in (a), and at the same time, it is compared and shown with the theoretical zero-order Bessel function. Figure 4 It contains 2 subgraphs. Figure 4 In (a), it is the spatial autocorrelation curve graph after spatial averaging when the station spacing is 15m. Figure 4 In (b), it is the spatial autocorrelation curve graph after spatial averaging when the station spacing is 30m. Both are the autocorrelation curves averaged by distance in space after DBSCAN clustering, divided into 3 pairs of station pairs with a distance of 15m each, namely station pairs 73 - 65, 73 - 70, 73 - 72, and 3 pairs of station pairs with a distance of 30m each, namely station pairs 73 - 37, 73 - 68, 73 - 71.

[0075] The method of arbitrarily specifying the central trace developed by the present invention provides a good solution when there are problems with the central trace station data and it affects dispersion imaging. Refer to Figure 5 and Figure 6 , Figure 5 are the dispersion energy spectrum and velocity curve obtained through the integration of manual and automatic methods with the 73rd station as the central station in the embodiment of the present invention. Among them, Figure 5 In (a), it is the dispersion energy spectrum calculated with the 73rd station as the central station. Figure 5 In (b), it is the phase velocity corresponding to the 73rd station as the central station and the calculated apparent shear wave velocity curve. Figure 6 are the dispersion energy spectrum and velocity curve obtained through the integration of manual and automatic methods with the 72nd non - central station in the embodiment of the present invention. Among them, Figure 6 In (a), it is the dispersion energy spectrum calculated with the 72nd non - central station. Figure 6 In (b), it is the phase velocity corresponding to the 72nd non - central station and the calculated apparent shear wave velocity curve. Through analysis, it can be seen that the dispersion imaging effect of the non - central station is very close to the effect calculated based on the theoretical method using the central station. It can be seen that this method is a new idea for standard SPAC calculation and very meets the actual production requirements.

[0076] Figure 7This is the dispersion fitting and model inversion result of the embodiment of the present invention adopting the adaptive gradient algorithm strategy. Among them, Figure 7 In (a), it is the comparison between the observed values and predicted values of frequency and phase velocity. Figure 7 In (b), it is the initial value and inversion result of shear wave velocity and depth. Using the method of the present invention to invert a phase velocity dispersion curve containing 36 dispersion points, the calculation takes 21.76 seconds, and the fitting error is 3.75 m / s. Compared with the traditional optimization algorithm, the adaptive gradient optimization strategy adopted by the present invention has inverted 32 dispersion curves, and the fitting error of each measurement point is smaller (Table 1). For inverting dispersion curves, the convergence effect and stability performance are better, and it is very suitable for the inversion work of microtremor data.

[0077] Table 1 Statistical table of inversion fitting errors of different inversion strategies

[0078] The present invention can display various dispersion spectra and velocity profiles. Figure 8 These are all the dispersion data points extracted from the stations in the embodiment of the present invention. Figure 8 Among them, circles represent phase velocity dispersion points, squares represent apparent shear wave velocity dispersion points, and triangles represent inverted shear wave velocity dispersion points. Different colors of each shape represent the change in their depth. Figures 9 - 11 This is the velocity (phase velocity, apparent shear wave velocity, inverted shear wave velocity) profile in the frequency domain of the embodiment of the present invention. Figures 12 - 14 This is the velocity (phase velocity, apparent shear wave velocity, inverted shear wave velocity) profile in the depth domain of the embodiment of the present invention. Figures 12 - 14 The positions of the acquisition stations are marked on the terrain line of the comprehensive profile in the depth domain. The present invention can directly guide the deployment work of field exploration boreholes, improving the basis and accuracy of borehole deployment.

[0079] The above description of the disclosed embodiments enables those skilled in the art to implement or use the present invention. Various modifications to these embodiments will be obvious to those skilled in the art, and the general principles defined herein can be implemented in other embodiments without departing from the spirit or scope of the present invention. Therefore, the present invention will not be limited to these embodiments shown herein, but will be accorded the widest scope consistent with the principles and novel features disclosed herein.

Claims

1. A data processing method for micro-motion detection array station layout, characterized in that: The following steps are involved: S1. According to the coordinates of the starting point measured by RTK and the azimuth of the survey line, the measuring point distance, the radius of the circumscribed circle of the nested triangle, and the number of nested triangles are set to generate the coordinates of all observation stations; S2, collecting the waveform data of the observation station corresponding to the coordinates, and preprocessing the waveform data; S3, calculating the spatial autocorrelation coefficient of the waveform data after preprocessing of any designated central station and non-central station; S4, perform DBSCAN clustering on the observation station pairs according to the distance between each pair of stations to obtain the clustering results; S5. According to the clustering results, the spatial autocorrelation coefficients of the observation station pairs are spatially averaged, the averaged spatial autocorrelation coefficients are fitted with the zero-order Bessel function to obtain the dispersion energy spectrum, the phase velocity dispersion data in the energy spectrum is picked up, and the apparent shear wave velocity dispersion curve is calculated; S6. Based on the adaptive gradient inversion strategy, the picked phase velocity dispersion data is inverted to obtain the shear wave velocity; S7. Generate a comprehensive profile in the depth domain and a comprehensive profile in the frequency domain based on the picked surface wave phase velocity, the calculated apparent shear wave velocity, and the inverted shear wave velocity, and load the station name or terrain on the profile to guide field exploration work.

2. The data processing method for micro-motion detection array station layout according to claim 1 is characterized in that: S1 is specifically: According to the given starting point coordinates, azimuth, distance and number of measuring points, a series of measuring point coordinates are generated, and nested triangles with a specified amount of data are drawn at each measuring point. The generation of nested triangles is based on the outer circle radius and the geometric characteristics of equilateral triangles. Multi-layer nesting is achieved by recursively calculating the midpoint. The center coordinates of all circumscribed circles and the vertex coordinates of the nested triangles are the coordinates of the observation station. When the number of nested triangles is set to 0, the observation stations are arranged linearly; when the number of nested triangles is set to 1~4, the observation stations are arranged in a two-dimensional array.

3. The data processing method for micro-motion detection array station layout according to claim 1 is characterized in that: Preprocessing methods include downsampling, de-averaging, de-trending, maximum value normalization or sliding absolute mean normalization, spectral whitening, and bandpass filtering.

4. The data processing method for micro-motion detection array station layout according to claim 1 is characterized in that: S4 is as follows: calculate the distance between any two pairs of stations to form a distance list, convert the distance list into an array, calculate the neighborhood radius eps of DBSCAN according to the allowed deviation percentage dev parameter: eps = (dev / 100) × mean_distance, where mean_distance is the mean of the distance list; extract the index and distance value of the station pair in each cluster, and the clustering results include: station name, station pair index, and distance list.

5. The data processing method for micro-motion detection array station layout according to claim 1 is characterized in that: The averaged spatial autocorrelation coefficient is fitted with the zero-order Bessel function to obtain the dispersion energy spectrum, and the phase velocity dispersion data in the energy spectrum is picked up. Specifically, through velocity scanning, the spatial autocorrelation coefficients of all clustered stations with different spacings at the same frequency are fitted with the first-kind zero-order Bessel function to generate a dispersion energy spectrum in which the frequency changes with the velocity. In the dispersion energy spectrum interface, press the left mouse button to trigger the picking function, obtain the frequency and speed analysis dispersion energy spectrum of the current mouse position, and find the speed point with the lowest energy at each frequency point. During the mouse dragging process, trigger the picking function. For each new frequency point, check whether it is necessary to continue picking the speed point with the lowest energy at the frequency point. Release the left mouse button to complete the range selection and complete the dispersion curve picking.

6. The data processing method for micro-motion detection array station layout according to claim 1 is characterized in that: According to the phase velocity dispersion data in the picked dispersion energy spectrum, the apparent shear wave velocity is calculated as follows: Apparent shear wave velocity = in, represents the frequency of the first dispersion point, represents the frequency of the second dispersion point, represents the phase velocity of the first dispersion point, represents the phase velocity of the second dispersion point.

7. The data processing method for micro-motion detection array station layout according to claim 1 is characterized in that: Adopting the adaptive gradient inversion strategy for inversion, the momentum term, adaptive parameter update strategy and line search mechanism are introduced to update the momentum term by the momentum attenuation factor: Where V represents the momentum term, is the search direction vector, represents the attenuation factor; Attenuation Factor Depending on the number of iterations, it decays as follows: The model is updated by selecting the step size that satisfies the Armijo condition through the line search mechanism; In the formula, is the search step length, For an updated stratigraphic model, Indicates the current formation model. is the parameter in Armijo, To search for gradient direction; If the Armijo criterion is not met, the current step size is adjusted using the dichotomy method until the condition is met; right To make adaptive adjustments: in, To adjust the rate, For smooth adjustment, represents the standard deviation of the function value sequence, k represents the number of loops, and maxIter represents the maximum number of iterations.

8. The data processing method for micro-motion detection array station layout according to claim 1 is characterized in that: The comprehensive profiles in the depth domain and the frequency domain are generated according to the picked surface wave phase velocity, the calculated apparent shear wave velocity, and the inverted shear wave velocity. Specifically: The Kriging interpolation method is used to generate a two-dimensional velocity profile based on the coordinate data and velocity data of the station. The y-axis is set to frequency to obtain the profile in the frequency domain, and the y-axis is set to depth to obtain the profile in the depth domain. The depth domain profile is corrected for elevation statics based on the elevation data.

9. A data processing system for micro-motion detection array station, characterized in that: include: The observation station generation unit is used to set the measurement point distance, the radius of the circumscribed circle of the nested triangle, and the number of nested triangles according to the starting point coordinates and the azimuth of the survey line measured by RTK, and generate the coordinates of all observation stations; A data preprocessing unit, used for collecting waveform data of observation stations corresponding to coordinates and preprocessing the waveform data; A spatial autocorrelation coefficient calculation unit is used to pair the observation stations in pairs and calculate the spatial autocorrelation coefficients of the observation station pairs based on the waveform data of the observation stations; The cluster analysis unit is used to perform DBSCAN clustering on observation station pairs according to the distance between each pair of stations to obtain clustering results; A dispersion curve generating unit is used to spatially average the spatial autocorrelation coefficients of the observation station pairs according to the clustering results, fit the averaged spatial autocorrelation coefficients with the zero-order Bessel function to obtain the dispersion energy spectrum, pick up the phase velocity dispersion data in the energy spectrum, and calculate the apparent shear wave velocity dispersion curve; An inversion unit, used to invert the picked phase velocity dispersion data based on an adaptive gradient inversion strategy to obtain the shear wave velocity; The profile generation unit is used to generate a comprehensive profile in the depth domain and a comprehensive profile in the frequency domain according to the picked surface wave phase velocity, the calculated apparent shear wave velocity, and the inverted shear wave velocity, and load the station name or terrain on the profile to guide field exploration work.

Images