A data processing method and system for the layout of a micro-motion detection array

Abstract

The present invention discloses a data processing method and system for the layout of a microseismic detection array, which relates to the technical field of microseismic detection. The method includes: setting the measuring point distance, the radius and quantity of the circumscribed circles of nested triangles according to the measured starting point coordinates and the azimuth angle of the survey line, and batch generating the coordinates of field stations; collecting the waveform data of the stations and preprocessing it; calculating the spatial autocorrelation coefficient between any central station and non-central stations; performing DBSCAN clustering based on the station distances; performing spatial averaging on the spatial autocorrelation coefficient according to the clustering results; fitting the averaged spatial autocorrelation coefficient with the zero-order Bessel function to obtain the dispersion spectrum, picking up the surface wave phase velocity in the dispersion spectrum, and calculating the apparent shear wave velocity; adopting an adaptive gradient inversion strategy to invert the shear wave velocity of the waveform data; generating a velocity comprehensive profile in the depth domain and frequency domain according to the surface wave phase velocity, the apparent shear wave velocity, and the inverted shear wave velocity. The present invention provides a new method for microseismic detection from station layout to data processing.

Description

Technical Field

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

[0002] Microseismic detection is developed based on the research of background noise imaging. 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, microseismic exploration is mostly used in engineering geological exploration and urban geological surveys, and some exploratory work in mineral exploration can also be seen. Although microseismic exploration technology has shown many advantages in geological detection in urban and complex environments, there are still several major technical or algorithmic problems that need to be solved urgently: (1) Currently, for the deployment of microseismic 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 microseismic detection, and the data formats supported by different calculation methods are single. For example, the commonly used calculation methods in the scientific research field support the SAC data format, while the engineering field methods mostly support the sg2 data format. The two lack compatibility, which 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 technicians cannot intelligently screen stations at different distances to participate in the calculation, resulting in the contamination of the dispersion energy spectrum by the data of stations with large interference; (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 comprehensive profile drawing 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 microseismic detection to solve the problems of insufficient adaptability and stability of the microseismic detection data processing algorithm in the existing geological exploration field, including the following steps:

[0004] 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 of the survey line, and generate the coordinates of all observation stations;

[0005] S2. Collect the waveform data of the observation stations corresponding to the coordinates, and preprocess the waveform data;

[0006] S3. Calculate the spatial autocorrelation coefficient of the preprocessed waveform data between any specified central station and non - central stations;

[0007] S4. Perform DBSCAN clustering on the pairs of observation stations according to the distance between each pair of stations to obtain the clustering result;

[0008] 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;

[0009] S6. Invert the picked - up phase velocity dispersion data based on the adaptive gradient inversion strategy to obtain the shear wave velocity;

[0010] 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 profile to guide the field exploration work.

[0011] Further, S1 is specifically as follows:

[0012] Generate a series of measuring point coordinates according to the given starting point coordinates, azimuth, measuring point distance and quantity, and draw nested triangles with a specified data volume at each measuring point. The generation of the nested triangles is based on the geometric properties of the outer circle radius and equilateral triangles, and multiple - layer nesting is realized by recursively calculating the mid - points. The center coordinates of all circumcircles and the vertex coordinates of the nested triangles are the coordinates of the observation stations;

[0013] 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.

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

[0015] Further, S4 is specifically as follows: Calculate the distances between every two station pairs to form a distance list, and convert the distance list into an array. According to the allowed 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 station pairs in each cluster. The clustering result includes: station names, indices of station pairs, and distance lists.

[0016] Further, fit the spatially averaged autocorrelation coefficient with the zero-order Bessel function to obtain the dispersion energy spectrum, and pick up the phase velocity dispersion data in the energy spectrum, specifically: Through velocity scanning, fit the spatially 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 with frequency varying with velocity. Under the dispersion energy spectrum interface, 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 for 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 velocity point with the lowest energy of the frequency point. Release the left mouse button to complete the range selection and finish picking the dispersion curve.

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

[0018] Apparent shear wave velocity =

[0019] 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.

[0020] Further, adopt an adaptive gradient inversion strategy for inversion, introduce a momentum term, an adaptive parameter update strategy, and a line search mechanism, and update the momentum term through the momentum decay factor:

[0021]

[0022] where, V represents the momentum term, is the search direction vector, represents the decay factor;

[0023] Decay factor Decays according to the following method based on the number of iterations:

[0024]

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

[0026]

[0027] 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;

[0028] If the Armijo criterion is not satisfied, the current step size is adjusted by the bisection method until the condition is met;

[0029] For perform adaptive adjustment:

[0030]

[0031]

[0032] Among them, is the adjustment rate, 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.

[0033] 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:

[0034] Using the Kriging interpolation method, based on the coordinate data and velocity data of the stations, a two-dimensional velocity profile is generated. Setting the y-axis as the frequency gives the profile in the frequency domain, and setting the y-axis as the depth gives the profile in the depth domain. Elevation static correction is performed on the profile in the depth domain according to the elevation data.

[0035] The present invention also proposes a data processing system for the microtremor detection array layout, including:

[0036] An observation station generation unit, configured to set the measuring 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;

[0037] A data preprocessing unit, configured to collect the waveform data of the observation stations corresponding to the coordinates and perform preprocessing on the waveform data;

[0038] A spatial autocorrelation coefficient calculation unit, configured to pair observation stations in pairs and calculate the spatial autocorrelation coefficient of pairs of observation stations according to the waveform data of the observation stations;

[0039] A clustering analysis unit, configured to perform DBSCAN clustering on pairs of observation stations according to the distances between pairs of stations to obtain a clustering result;

[0040] A dispersion point generation unit, configured to perform spatial averaging on the spatial autocorrelation coefficients of pairs of observation stations according to the clustering result, fit the averaged spatial autocorrelation coefficient and 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;

[0041] 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;

[0042] 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 for guiding field exploration work.

[0043] The beneficial effects brought by the technical solution provided by the present invention are:

[0044] According to the measured coordinates of the measuring points and the azimuth of the survey line, the present invention generates the coordinates of all observation stations according to a nested triangle or a linear observation system; performs DBSCAN clustering on the observation stations, and according to the clustering result, averages the spatial autocorrelation coefficients of the same clustering clusters to improve the calculation effect of the dispersion energy spectrum; generates a comprehensive profile in the depth domain and a comprehensive profile in the frequency domain by batch interpolation according to the picked-up 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 for guiding 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

[0045] Figure 1 is a flowchart of the data processing method for the microtremor detection array layout in the embodiment of the present invention;

[0046] Figure 2 is a spatial distribution map of the acquisition of stations generated in the embodiment of the present invention;

[0047] Figure 3 is a 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 (b) is the seismic waveform data preprocessed according to the preprocessing method of the present invention;

[0048] Figure 4 is the spatial autocorrelation curve diagram generated in the embodiment of the present invention; Figure 4 In (a) is the spatial autocorrelation curve diagram after spatial averaging when the station spacing is 15m, Figure 4 In (b) is the spatial autocorrelation curve diagram after spatial averaging when the station spacing is 30m;

[0049] 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. Among them, Figure 5 In (a) is the dispersion energy spectrum calculated with the 73rd station as the central station, Figure 5 In (b) is the phase velocity corresponding to the 73rd station as the central station and the calculated apparent shear wave velocity curve;

[0050] 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. Among them, Figure 6 In (a) is the dispersion energy spectrum calculated with the 72nd non - central station, Figure 6 In (b) is the phase velocity corresponding to the 72nd non - central station and the calculated apparent shear wave velocity curve;

[0051] Figure 7 is the dispersion fitting and model inversion result in the embodiment of the present invention. Among them, Figure 7 In (a) is the comparison between the observed and predicted values of frequency and phase velocity, Figure 7 In (b) is the initial value and inversion result of shear wave velocity and depth;

[0052] Figure 8 are all the dispersion data points extracted in the embodiment of the present invention;

[0053] Figure 9 is the phase velocity profile in the frequency domain in the embodiment of the present invention;

[0054] Figure 10 is the apparent shear wave velocity profile in the frequency domain in the embodiment of the present invention;

[0055] Figure 11 is the inverted shear wave velocity profile in the frequency domain in the embodiment of the present invention;

[0056] Figure 12 is the phase velocity profile in the depth domain in the embodiment of the present invention;

[0057] Figure 13 is the apparent shear wave velocity profile in the depth domain in the embodiment of the present invention;

[0058] Figure 14 It is the shear wave velocity profile of the deep domain inversion in the embodiment of the present invention. Specific Embodiments

[0059] 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.

[0060] The flowchart of the data processing method for the microseismic detection array layout in the embodiment of the present invention is as Figure 1 , and specifically includes the following steps:

[0061] S1. According to the starting point coordinates measured by RTK (Real-Time Kinematic) 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.

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

[0063] Initialize 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, counterclockwise is 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.

[0064] Measuring point generation: Generate a series of uniformly distributed measuring point coordinates along the specified direction according to the starting point coordinates, azimuth angle, measuring point distance, and number of measuring points: Convert the azimuth angle angle from degrees to radians, and the formula is:

[0065]

[0066] For the i-th measuring point, its east coordinate x and north coordinate y are calculated respectively through the following formulas:

[0067]

[0068]

[0069] Triangle generation: Generate the coordinates of the three vertices of the initial equilateral triangle according to the outer circle radius, and generate the nested triangle levels as needed. The calculation idea is that the included angle between the three vertices of the equilateral triangle is 120 degrees, and by adjusting the angle_offset, the overall direction of the triangle can be rotated:

[0070]

[0071]

[0072] where i is 0, 1, or 2.

[0073] 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.

[0074] Two points 、 The midpoint calculation formula is: which is used to generate the three vertices of the next-level nested triangle.

[0075] Drawing nested triangles: According to the coordinates generated by the measuring points, draw nested equilateral triangles and number them. The calculation idea is: Based on the coordinates of each measuring 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, and recursively generate num_nested layers of nested triangles. Number and label each measuring point and the vertices of its nested triangles.

[0076] Exporting coordinate point data: Write the east and north coordinates of all measuring points and triangle vertices line by line to a CSV file for subsequent field measurements.

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

[0078] By exporting the calculated coordinates, the rationality of the calculation can be detected and corrected on the Ovi map. 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 arranged in nested triangles, the SPAC or ESPAC algorithm can be calculated according to the characteristics of the nested arrangement. At the same time, for the nested triangle stations, a certain station can be selected as the center point, which greatly improves the rule that only the station with the smallest number is defaulted as the center point in the previous processing technology, and also avoids the artificial limitation that the stations are deployed clockwise or counterclockwise in some processing technologies.

[0079] S2. Collect the waveform data of the observation stations corresponding to the coordinates, read the microtremor data in multiple formats (sac, sg2), and preprocess the waveform data.

[0080] The present invention supports waveform data in sac and sg2 formats, writes the station number into the sac file header, solves the problem of station information loss in station data, and can perform correct subsequent data processing by matching the station number in the coordinate file with the data with station numbers in the waveform data.

[0081] When collecting waveform data of observation stations, the observation time is long and the data volume is large. For many systems, directly importing the original data often has limited computing memory and is difficult to process. Downsampling the original waveform to a reasonable sampling rate range makes the inner layer of a single waveform file meet the requirements of computing memory and computing speed, greatly improving the work efficiency.

[0082] The present invention designs a variety of optional preprocessing methods, including downsampling (optional), de-meaning, detrending, maximum normalization or sliding absolute average trace normalization, spectral whitening, and band-pass filtering. Any one of the processing methods can be manually selected. By comparing the changes in the waveform data after preprocessing with the original waveform data, the processing method that conforms to the geological background and acquisition background can be selected. The automated and anti-aliasing batch downsampling 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 downsampling is ensured.

[0083] The implementation process of the preprocessing method is as follows:

[0084] Downsampling is to reduce the sampling rate of a digital signal and reduce the number of sampling points per unit time. For example, if the original signal sampling rate is , and the downsampling factor is M , then the new sampling rate after downsampling is which is:

[0085]

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

[0087]

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

[0089] The purpose of de-meaning is to eliminate the DC component in the signal and make the average value of the signal zero. This 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 signal after de-meaning is The calculation is as follows:

[0090] ,

[0091] ,

[0092] The purpose of detrending is to eliminate the linear trend in the signal, making the signal more stable and facilitating subsequent processing. This is achieved by fitting the trend line of the signal and subtracting it from the signal. For example, if the signal after mean removal is , the fitted linear trend , where and are determined as the slope and intercept through the least squares method. The detrended signal is calculated as follows:

[0093] ,

[0094] Moving absolute average normalization realizes local normalization of the signal by calculating the moving average of the absolute values 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 , .

[0095] 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 maximum value normalization can be expressed as:

[0096] ,

[0097] Spectral whitening aims to flatten the spectrum of the 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, perform amplitude normalization on the frequency-domain signal, and then perform an inverse Fourier transform to obtain the whitened signal in the time domain. For example, if the input signal is the normalized signal , its Fourier transform is , and the frequency-domain signal after spectral whitening is calculated as: , and performing an inverse Fourier transform on it gives the time-domain data after spectral whitening: .

[0098] 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 unexpected situations such as the central station device malfunctioning, being lost, or damaged, it will not be possible to meet the principle of the SPAC algorithm to carry out work, which may lead to the invalidation of the entire acquisition work. In the present invention, the central station can be arbitrarily specified, and calculating SPAC is more in line with the actual production scenario.

[0099] 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 traditional clustering analysis methods 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.

[0100] 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 at least one sample point in each cluster.

[0101] 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 cluster index list and distance value list in idx_list and all_list respectively. Through the DBSCAN clustering method, the clustering results not only include the distance list and the indices of the station pairs, but also can be associated with the specific station names, facilitating subsequent calculation of the spatial autocorrelation coefficient and the dispersion energy spectrum.

[0102] S5. Perform spatial averaging on 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.

[0103] The process of performing a dispersion energy spectrum scan and picking up the dispersion curve using manual and semi-automatic methods includes:

[0104] Through velocity scanning, fit the spatial autocorrelation coefficients of different station spacings after clustering at the same frequency with the first kind of zero-order Bessel function to generate a dispersion energy spectrum with frequency varying with velocity. The calculation formula for the dispersion energy spectrum is:

[0105]

[0106] Among them, fitting error, is the spatial autocorrelation coefficient, is the first kind of zero-order Bessel function, , is the spatial distance of the jth station. Under the dispersion energy spectrum interface, the user presses the left mouse button to trigger the picking function, obtains the frequency and velocity at the current mouse position to analyze the dispersion energy spectrum, automatically searches for the velocity point with the lowest energy near each frequency point. If a valid point is found, add it to the list and update the data and graph. During the mouse dragging process, trigger the picking function 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, call the "find the lowest energy value" function to automatically find the lowest energy point and add it to the list, and update the data and graph. When the user releases the left mouse button, automatically save all the picked dispersion points to a file. This method reduces manual operations, improves efficiency and accuracy, and solves the problem of the previous need for frequent manual and automatic switching operations. The present invention still retains the function of manually adding and deleting points, providing flexibility and controllability. At the same time, introduce the sliding window technology during the picking process to perform mean smoothing processing on the picked dispersion points, reduce the influence of noise, enhance the data continuity and reliability; based on the automatically picked dispersion points, use the relationship between velocity and frequency to automatically calculate the apparent shear wave velocity to ensure the physical rationality of the calculation results.

[0107] Apparent shear wave velocity =

[0108] Among them, 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.

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

[0110] According to the adaptive gradient algorithm strategy, introduce the momentum term, the adaptive parameter update strategy, and the line search mechanism. The search direction is mainly determined by the momentum decay factor Update the momentum term. The calculation formula is:

[0111]

[0112] 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 in the following way:

[0113]

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

[0115]

[0116] In the formula, 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 will be adjusted by the bisection method until the condition is met.

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

[0118] If the volatility of the objective function is large, it can be suppressed by reducing . That is:

[0119]

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

[0121] This method will Set in the following form with the aim of optimizing in the initial stage, where the adjustment speed of the step size is fast to avoid premature convergence, and in the later stage, the step size adjustment becomes more meticulous.

[0122]

[0123] Among them, k represents the number of cycles, and maxIter represents the maximum number of iterations.

[0124] 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, calculated apparent shear wave velocity, and inverted shear wave velocity. Optionally, the profile can load the station name or topography to guide the field exploration work.

[0125] Read the station coordinate file, including the station name, east coordinate, north coordinate, and elevation coordinate; read each velocity file, where the velocity file includes: picked surface wave phase velocity data, calculated apparent shear wave velocity data, and inverted shear wave velocity data;

[0126] Using the Kriging interpolation method, generate a two-dimensional velocity profile based on 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.

[0127] In an exemplary embodiment, it includes a data processing system for the layout of a microtremor detection array, including:

[0128] An observation station generation unit for setting 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 generating the coordinates of all observation stations;

[0129] A data preprocessing unit for collecting the waveform data of the observation stations corresponding to the coordinates and preprocessing the waveform data;

[0130] A spatial autocorrelation coefficient calculation unit for pairing the observation stations in pairs and calculating the spatial autocorrelation coefficient of the pairs of observation stations according to the waveform data of the observation stations;

[0131] A clustering analysis unit for performing DBSCAN clustering on the pairs of observation stations according to the distance between the pairs of stations to obtain the clustering result;

[0132] A dispersion point generation unit for spatially averaging the spatial autocorrelation coefficients of the pairs of observation stations according to the clustering result, fitting the averaged spatial autocorrelation coefficient with the zero-order Bessel function to obtain the dispersion energy spectrum, picking up the phase velocity dispersion data in the energy spectrum, and calculating the apparent shear wave velocity dispersion curve;

[0133] An inversion unit, which is used to invert the picked phase velocity dispersion data based on an adaptive gradient inversion strategy to obtain the shear wave velocity;

[0134] 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.

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

[0136] 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 3 in which (b) is the seismic waveform data after mean removal, trend removal, sliding 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.

[0137] Figure 4 is based on Figure 3 the spatial autocorrelation curve calculated from the microtremor preprocessed waveform data in, and is also compared and shown with the theoretical zero-order Bessel function. Figure 4 contains 2 subgraphs, Figure 4 in which (a) is the spatially averaged spatial autocorrelation curve graph when the station spacing is 15m, Figure 4Figure (b) is the spatial autocorrelation curve after spatial averaging when the station spacing is 30 m. They are all autocorrelation curves averaged by distance after DBSCAN clustering, and are divided into 3 pairs of station pairs with a distance of 15 m each, namely station pairs 73 - 65, 73 - 70, and 73 - 72, and 3 pairs of station pairs with a distance of 30 m each, namely station pairs 73 - 37, 73 - 68, and 73 - 71.

[0138] The method of arbitrarily specifying the central channel developed in the present invention provides a good solution when there are problems with the central channel station data and it affects dispersion imaging. Refer to Figure 5 and Figure 6 , Figure 5 Figures (a) and (b) are the dispersion energy spectrum and velocity curve obtained through the integration of manual and automatic operations with the 73rd station as the central station in the embodiment of the present invention. Among them, Figure 5 Figure (a) is the dispersion energy spectrum calculated with the 73rd station as the central station, Figure 5 Figure (b) is the phase velocity corresponding to the 73rd station as the central station and the calculated apparent shear wave velocity curve. Figure 6 Figures (a) and (b) are the dispersion energy spectrum and velocity curve obtained through the integration of manual and automatic operations with the 72nd non - central station in the embodiment of the present invention. Among them, Figure 6 Figure (a) is the dispersion energy spectrum calculated with the 72nd non - central station, Figure 6 Figure (b) is the phase velocity corresponding to the 72nd non - central station and the calculated apparent shear wave velocity curve. It can be seen from the analysis that the dispersion imaging effect of the non - central station is very close to the effect calculated based on the theoretical method, that is, 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.

[0139] Figure 7 Figures (a) and (b) are the dispersion fitting and model inversion results of the embodiment of the present invention using the adaptive gradient algorithm strategy. Among them, Figure 7 Figure (a) is the comparison between the observed and predicted values of frequency and phase velocity, Figure 7 Figure (b) 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 the inversion of dispersion curves, the convergence effect and stability performance are better, and it is very suitable for the inversion work of micro - motion data.

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

[0141]

[0142] The present invention can display various dispersion spectra and velocity profiles.Figure 8 All the station dispersion data points extracted in the embodiments of the present invention Figure 8 In it, 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 its depth. Figures 9 - 11 is the frequency domain velocity (phase velocity, apparent shear wave velocity, inverted shear wave velocity) profile of the embodiments of the present invention. Figures 12 - 14 is the depth domain velocity (phase velocity, apparent shear wave velocity, inverted shear wave velocity) profile of the embodiments of the present invention. Figures 12 - 14 The positions of the acquisition stations are marked on the terrain line of the depth domain comprehensive profile. The present invention can directly guide the field exploration drilling deployment work, improving the basis and accuracy of the drilling deployment.

[0143] 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 apparent 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 rather to the widest scope consistent with the principles and novel features disclosed herein.

Claims

1. A data processing method for the layout of a micro-motion detection array, characterized in that It includes 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 measuring 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 each pair 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 the comprehensive profile in the depth domain and the 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.

2. The data processing method for the layout of a micro-motion detection array according to claim 1, wherein, Specifically, S1 is 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 outer circle radius and equilateral triangles, and multiple layers of nesting are realized through 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.

3. The data processing method for the layout of a micro-motion detection array according to claim 1, characterized in that The preprocessing methods include downsampling, de-meaning, detrending, maximum normalization or moving absolute average normalization, spectral whitening, and band-pass filtering.

4. The data processing method for the layout of a micro-motion detection array according to claim 1, characterized in that Specifically, S4 is as follows: Calculate the distances between each pair of stations to form a distance list, convert the distance list into an array, and calculate the neighborhood radius eps of DBSCAN according to the allowable deviation percentage dev parameter: 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 lists.

5. The data processing method for the layout of a micro-motion detection array according to claim 1, characterized in that Fit the averaged spatial autocorrelation coefficient with the zero-order Bessel function to obtain the dispersion energy spectrum, and pick up the phase velocity dispersion data in the energy spectrum. Specifically: 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 the dispersion energy spectrum of frequency varying with velocity; Under the interface of the dispersion energy spectrum, click the left mouse button to trigger the picking function, obtain the frequency and velocity of the current mouse position to analyze the dispersion energy spectrum, find the velocity 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 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.

6. A data processing method for a micro-motion detection array station layout according to claim 1, characterized in that According to the phase velocity dispersion data in the picked dispersion energy spectrum, the apparent shear wave velocity calculation formula is as follows: Apparent shear wave velocity = Among them, 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. A data processing method for the layout of a micro-motion detection array according to claim 1, characterized in that An adaptive gradient inversion strategy is adopted for inversion. A momentum term, an adaptive parameter update strategy, and a line search mechanism are introduced. The momentum term is updated by the momentum decay factor: where V represents the momentum term, is the search direction vector, represents the decay factor; Attenuation factor Attenuate according to the number of iterations in the following manner: The step size that satisfies the Armijo condition is selected through the line search mechanism to update the model; wherein, 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; For perform adaptive adjustment: wherein, is the adjustment rate, 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.

8. The data processing method for the layout of a micro-motion detection array according to claim 1, characterized in that, Based on 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, according to 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.

9. A data processing system for the layout of a micro - motion detection array, characterized in that, Including: An observation station generation unit, which is used to set the measuring 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 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 pair 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 pair of observation stations according to the distance between the pairs of stations to obtain the clustering result; A dispersion curve generation unit, which is used to perform spatial averaging on the spatial autocorrelation coefficient of the pair 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.

Images