A method for three-dimensional positioning of a seismic source based on multi-channel amplitude ratio

Abstract

The application discloses a three-dimensional positioning method of a seismic source based on a multi-channel amplitude ratio. The method obtains waveform data of the same event on at least four sensor channels and spatial coordinates of the sensors, pre-processes the waveforms, determines an event window corresponding to the event, extracts amplitude indicators of each channel in the event window and performs channel consistency correction, constructs an amplitude logarithm ratio observation on an arbitrary channel based on an exponential decay model to eliminate a source intensity factor, and further forms a constraint related to a distance difference of the seismic source and the sensor. A redundant channel constraint set is further constructed, and an abnormal channel is down-weighted or removed. A weighted nonlinear least square is used in combination with an iterative reweighted least square or a robust loss function to solve the problem, a residual threshold is determined and is updated with iteration. A plurality of candidate solutions are obtained by solving a multi-initial value or a coarse grid initial value in a preset feasible region, and an optimal solution is output. The method reduces the dependence on the arrival time picking, and can be used for microseismic monitoring, acoustic emission positioning and underground engineering safety monitoring.

Description

Technical Field

[0001] This invention relates to a three-dimensional source localization method, and more particularly to a three-dimensional source localization method based on multi-channel amplitude ratio. Background Technology

[0002] Existing source location methods are mostly based on arrival time information, such as P-wave first arrival picking, travel time inversion, and cross-correlation estimation. These methods typically require clear first arrivals and a high signal-to-noise ratio. In complex media, multipath propagation, near-field saturation, poor coupling, or strong noise backgrounds, arrival time picking is prone to significant errors, leading to decreased location accuracy or even location failure. Wave propagation in a medium involves absorption attenuation, scattering energy dissipation, and geometric diffusion effects, causing the observed amplitude to vary with propagation distance. However, in engineering monitoring scenarios, directly utilizing multi-channel amplitude information for location typically faces the following technical challenges:

[0003] (1) Differences in sensor gain, sensitivity and coupling state will introduce systematic amplitude deviation, resulting in incomparable amplitudes of different channels, thus causing bias in amplitude constraints;

[0004] (2) Amplitude attenuation includes not only exponential absorption terms, but also geometric diffusion terms and other distance-related terms. If only a single distance difference mapping is used, it is easy to cause model mismatch and reduce positioning accuracy.

[0005] (3) Under conditions of low signal-to-noise ratio, near-field saturation or channel anomalies, the amplitude extraction of individual channels or individual channel pairs is prone to outlier constraints. If there is a lack of redundant constraint quality control and robust solution mechanism, the localization results are prone to non-convergence or instability.

[0006] (4) When absolute value constraints or absolute value residuals are used to adapt to phase / polarity instability, multiple solutions or mirror solutions will be introduced, requiring corresponding feasible region constraints and disambiguation strategies.

[0007] Therefore, it is necessary to propose a three-dimensional source localization method that can achieve low computational cost and stable positioning while ensuring amplitude comparability, being compatible with geometric diffusion models, and suppressing the influence of outlier constraints through redundancy constraint consistency control and robust optimization. Summary of the Invention

[0008] The purpose of this invention is to provide a three-dimensional source localization method based on multi-channel amplitude ratios. This addresses the shortcomings of existing source localization methods, which often rely heavily on arrival time information (such as first arrival picking and travel time inversion). In complex media, multipath propagation, strong noise backgrounds, near-field saturation, or poor coupling conditions, arrival time picking is prone to significant errors, leading to decreased localization accuracy or even failure. Furthermore, localization schemes based on amplitude attenuation face challenges in engineering applications, including amplitude incomparability due to differences in channel gain and coupling, attenuation model mismatch (e.g., neglecting geometric diffusion terms), and solution instability caused by anomalous channels / outlier constraints. In contrast, the method of this invention offers a three-dimensional source localization method with lower computational complexity, weaker dependence on arrival time picking, and stable convergence even in the presence of anomalous channels and complex noise.

[0009] The technical solution of this invention: a three-dimensional source localization method based on multi-channel amplitude ratio, comprising the following steps:

[0010] a. Acquire waveform data of the same event on at least four sensor channels. and the spatial coordinates of at least four sensors. ;

[0011] b. Perform preprocessing on each channel waveform and determine the same event window Ω; within the event window Ω, perform multi-channel alignment based on the reference channel k to ensure that the amplitude extraction of different channels corresponds to the same event wave packet;

[0012] c. Extract the amplitude index A of each channel within the event window Ω. i And for amplitude index A i Perform channel consistency correction to obtain the corrected amplitude :

[0013] ,

[0014] ,

[0015] Where c i This is the consistency correction coefficient;

[0016] d. Based on the exponential decay model Constructing logarithmic ratio observations and establishing residuals, firstly, for any two channels, the amplitude index A... i A j Constructing logarithmic ratio observations :

[0017] ,

[0018] Based on the exponential decay model, the amplitude indices of channel i and channel j are compared and logarithmically transformed to eliminate the unknown source strength factor A0 of the event, thus obtaining the distance difference observable. :

[0019] ,

[0020] Where r i Let the propagation distance from the earthquake source to the i-th sensor be: P is the location of the earthquake source to be determined, and α is the attenuation coefficient;

[0021] When the exponential decay model further includes a geometric diffusion term When γ is the preset or calibrated geometric diffusion index; the residual term consists of a term related to the distance difference and a term related to the distance ratio, so that the logarithmic ratio observation is consistent with the attenuation model, at which point the maximum amplitude is... Determine by the following formula:

[0022] ,

[0023] Where A0 is the source strength factor, r i γ is the propagation distance from the earthquake source to the i-th sensor, γ is the geometric diffusion index, and α is the attenuation coefficient;

[0024] e. Construct a redundant channel constraint set And establish the distance difference residual:

[0025] ,

[0026] Channel pairs are screened, downweighted, or eliminated based on the consistency of residuals between them as a quality control criterion.

[0027] f. After screening, weighting, or elimination, establish geometric distance relationships for the epicenter location P and construct an objective function. Solve for the epicenter location using (weighted) nonlinear least squares to minimize the objective function.

[0028] ,

[0029] in Let ρ be the channel pair weights, and ρ() be the robust loss function;

[0030] The solution employs iterative reweighted least squares or introduces a robust loss function to reduce the impact of outlier channels on the solution results;

[0031] g. The residual threshold is adaptively determined based on the noise statistics or residual robust statistics estimated from the background noise segment, and updated it iteratively to achieve automatic weight reduction or elimination of outlier constraints;

[0032] h. Output the estimated seismic source location. And positioning quality indicators, which should include at least the minimum value of the objective function or the root mean square of the residuals.

[0033] In the aforementioned three-dimensional source localization method based on multi-channel amplitude ratio, the amplitude index A of the i-th channel... i Determine as follows:

[0034] Take the absolute peak value of the channel waveform within the event window Ω:

[0035] ,

[0036] Alternatively, perform a Hilbert transform on the waveform to obtain the envelope, and then take the envelope peak value.

[0037] ,

[0038] Where H{} is the Hilbert transform, This represents the waveform of the i-th channel.

[0039] In the aforementioned three-dimensional source localization method based on multi-channel amplitude ratio, when the distance difference observation in step d is expressed in absolute value form:

[0040] ,

[0041] Furthermore, step f uses the corresponding absolute value residual:

[0042] .

[0043] In the aforementioned three-dimensional source localization method based on multi-channel amplitude ratio, the determination of the attenuation coefficient α includes at least one of the following methods:

[0044] Based on the dominant frequency or center frequency f, the propagation speed v, and the quality factor Q, the following is calculated:

[0045] ,

[0046] Alternatively, the location can be determined by calibration events at known positions. Select several calibration events whose locations are known (or whose locations can be determined with high precision by other methods), and record their distance to the i-th sensor as r. i The amplitude is A i In the model without a geometric diffusion term, we have:

[0047] ,

[0048] For lnA i With r iLinear regression can be used to fit α (and simultaneously obtain lnA0). To reduce the influence of outliers, robust regression (such as Huber regression) is used. If a model with a geometric diffusion term is used, the fitted equation is:

[0049] ,

[0050] Thus, α and γ are jointly fitted in the same regression.

[0051] In the aforementioned three-dimensional source localization method based on multi-channel amplitude ratio, the attenuation coefficient α is a partitioning parameter or a layering parameter and is determined using a rolling update method. When the medium in the monitoring area is non-homogeneous, the area is divided into several sub-regions R. m α is determined for each sub-region. m When a new event arrives, if its localization residual is consistently large, then the corresponding sub-region's α is updated using this event and the M most recent high-quality events as samples. m Updated to:

[0052] ,

[0053] in The estimate is obtained based on the regression of the new sample, and λ∈(0,1] is the update coefficient.

[0054] In the aforementioned three-dimensional source localization method based on multi-channel amplitude ratio, the attenuation coefficient α is jointly estimated with the source location P. When α is uncertain and sensitive to location, α will be used as a parameter to be estimated and inverted together with P:

[0055] .

[0056] In the aforementioned three-dimensional source localization method based on multi-channel amplitude ratio, the redundant constraint set includes at least N-1 channel pairs constructed based on the reference channel, and the channel pairs are screened, weighted down, or eliminated through residual consistency to ensure that the number of effective constraints is not less than a preset threshold.

[0057] In the aforementioned three-dimensional source localization method based on multi-channel amplitude ratio, the event window is determined by any one or more of STA / LTA, energy mutation detection, threshold triggering, or cross-correlation coarse alignment.

[0058] In the aforementioned three-dimensional source localization method based on multi-channel amplitude ratio, the channel consistency correction is performed because differences in the sensitivity, coupling conditions, and gain settings of different channel sensors can lead to incomparable amplitudes of different channels under the same event. To ensure that the amplitude ratio reflects the propagation attenuation law, consistency correction is performed on the channels after amplitude extraction.

[0059] Let the observed amplitude of channel i be A.i The consistency correction factor is c. i The corrected amplitude is defined as:

[0060] ,

[0061] Correction factor c i This is used to compensate for differences in channel sensitivity, so that the amplitudes of different channels are in the same dimension and scale under the same propagation conditions.

[0062] (1) Coefficient determination based on calibration events

[0063] Using a set of calibrated events with known locations For each event Calculate its geometric distance r to each channel. i,e Under the exponential decay model, the expected outcome is:

[0064] ,

[0065] Where r i,e Let be the geometric distance from the source location of event e to the i-th sensor. Let e ​​be the source strength factor. Let be the corrected amplitude index for event e on the i-th channel.

[0066] With a fixed reference channel k, let c k =1. For each channel i ≠ k, estimate c by minimizing the following equation. i :

[0067] ,

[0068] in Let r be the amplitude index extracted from event e on the i-th channel. k,e Let be the geometric distance from the source location of event e to the sensor in reference channel k.

[0069] To avoid the impact of abnormal events, median statistics or robust regression can be used for the event set.

[0070] (2) Reference normalization under calibration conditions

[0071] When calibration events are lacking, the relative gain is obtained by statistically analyzing multiple events within a time window. With a fixed reference channel k, for each channel i ≠ k, the following definition is used:

[0072] ,

[0073] Where exp() is an exponential function with base e. To calculate the median for event index e.

[0074] Even in multiple events A i With A k This method achieves the same scale in a statistical sense. It is independent of the epicenter location but requires the number of events to meet a preset lower limit E. min (Minimum number of events included in the statistics, for example, no less than 20).

[0075] In the aforementioned three-dimensional source localization method based on multi-channel amplitude ratio, the weights in the weighted nonlinear least squares method... Determined according to the following formula:

[0076] ,

[0077] SNR i Let be the signal-to-noise ratio of the i-th channel. A positive constant is preset to prevent the denominator from being zero.

[0078] In the aforementioned three-dimensional source localization method based on multi-channel amplitude ratio, outlier channels or outlier constraints are removed during the solution process, or a robust loss function is used to reduce the impact of outliers.

[0079] ,

[0080] Where ρ() is the robust loss function. Furthermore, iterative reweighted least squares (IRLS) can be used: in the l-th iteration, with...

[0081] ,

[0082] Update the equivalent weights, where ω(⋅) is the weight function derived from the selected robust loss.

[0083] In the aforementioned three-dimensional source localization method based on multi-channel amplitude ratio, the root mean square residual used for anomaly detection is determined as follows:

[0084] .

[0085] In the aforementioned three-dimensional source localization method based on multi-channel amplitude ratio, when the number of sensor channels is four, the spatial arrangement of the four sensors satisfies the geometric non-degradation condition. The geometric non-degradation condition includes: the four points are not coplanar or the volume of the geometric body formed by the sensor coordinates is greater than a preset threshold. If the geometric non-degradation condition is not met, the number of sensor channels is increased, or a feasible domain constraint of the source space is introduced and an uncertainty prompt is output.

[0086] In the aforementioned three-dimensional source localization method based on multi-channel amplitude ratio, when the sensor channels are four channels 1, 2, 3, and 4, the constructed distance difference set includes at least one of the following:

[0087] ,

[0088] ,

[0089] And respectively, the amplitude was calculated as follows:

[0090] ,

[0091] .

[0092] In the aforementioned three-dimensional source localization method based on multi-channel amplitude ratio, step f uses weighted nonlinear least squares or weighted robust minimization to solve for the source location P, and superimposes the source space feasible domain constraint P∈V.

[0093] In the aforementioned three-dimensional source localization method based on multi-channel amplitude ratio, when using absolute value form, multiple candidate solutions are obtained by using multiple initial values ​​or coarse grid initial value strategies within the preset feasible domain of the source space. The final solution is selected based on the criterion of minimizing the objective function value and satisfying the feasible domain constraint. When the objective function values ​​of multiple candidate solutions are close, the confidence level or uncertainty index is output.

[0094] In the aforementioned three-dimensional source localization method based on multi-channel amplitude ratio, the residual threshold τ in step g is determined based on the statistics of the background noise segment outside the event window, satisfying:

[0095] ,

[0096] Where σ is the standard deviation of the background noise segment, and β is the preset multiplier coefficient;

[0097] When the residual threshold τ is determined based on robust statistics of the residuals, the following conditions are met:

[0098] ,

[0099] in β is a preset coefficient.

[0100] In the aforementioned three-dimensional source localization method based on multi-channel amplitude ratio, the residual threshold τ in step g is updated during the iterative solution process, satisfying:

[0101] ,

[0102] in λ is the threshold calculated based on the residual statistics in the l-th iteration, and λ∈(0,1] is the update coefficient.

[0103] The beneficial effects of the present invention: Compared with the prior art, the present invention has at least the following beneficial effects:

[0104] (1) Reduce reliance on first arrival picking and improve localization under weak signal conditions. This invention does not require first arrival picking of waveforms in each channel, and can construct localization constraints based on the maximum amplitude extracted within the event window; therefore, in scenarios where the first arrival is unclear, the signal-to-noise ratio is low, the reflection / scattering component is dominant, or the consistency of manual picking is poor, it can still form constraint information that can be used for three-dimensional inversion, thereby improving the availability and engineering adaptability of localization.

[0105] (2) Channel consistency correction improves amplitude comparability and reduces systematic amplitude deviation. By performing channel consistency correction on the amplitude index within the event window, the systematic deviation caused by differences in sensor sensitivity, gain settings and coupling conditions is compensated, making the amplitudes of different channels comparable. This reduces the impact of amplitude offset on distance difference constraints from the source and improves positioning accuracy and cross-event consistency.

[0106] (3) By using redundant channels for consistency screening and robust optimization, the influence of abnormal channels and outlier constraints is suppressed, thereby improving the convergence stability in scenarios with low signal-to-noise ratio, saturation, or poor coupling. This invention constructs a set of redundant channel pairs as constraints and uses the consistency of channel pair residuals as a quality control criterion to reduce the weight of or eliminate abnormal channel pairs. This mechanism can effectively suppress outlier constraints caused by individual channel saturation, shearing, transient interference, coupling failure, or amplitude extraction errors, avoiding the problem of "a few abnormal constraints dominating the solution," thereby improving the localization convergence stability and the reliability of the results.

[0107] (4) Low computational load, suitable for real-time or near-real-time positioning. The present invention uses amplitude extraction, distance difference constraint construction and (weighted) least squares optimization as the main computational processes. The overall process structure is clear, the number of parameters is small, and efficient iterative solution can be achieved. Therefore, it can achieve rapid positioning on conventional computing platforms and is suitable for real-time or near-real-time source monitoring, early warning and online analysis scenarios. Attached Figure Description

[0108] Figure 1 This is a schematic diagram of the three-dimensional layout of the sensors;

[0109] Figure 2 Define the schematic diagram for the waveform and event window;

[0110] Figure 3 The flowchart shows the three-dimensional source localization method based on multi-channel amplitude ratio redundancy constraints and robust optimization.

[0111] Figure reference numerals: 1—Sensor; 2—Distance r from the seismic source to the sensor; 3—Seismic source; 4—Event window Ω; 5—Amplitude index A i t0—Reference time; t p —Reserved lead time; t s —Post-time. Detailed Implementation

[0112] The present invention will be further described below with reference to the accompanying drawings and embodiments, but this should not be construed as limiting the present invention.

[0113] This invention provides a three-dimensional source localization method based on amplitude attenuation, which includes at least the following steps: data acquisition, event window determination, maximum amplitude extraction, attenuation coefficient determination, distance difference constraint construction, simultaneous solution, and output. For ease of understanding, the following text first provides notation conventions, followed by an example.

[0114] Symbol and terminology conventions:

[0115] i=1, …., N: Sensor (or detector) channel number, N≥4.

[0116] : Spatial coordinates of the i-th sensor.

[0117] Euclidean norm.

[0118] : Distance from the earthquake source to the i-th sensor.

[0119] t: Time variable.

[0120] : Waveform data of the i-th channel.

[0121] Ω: Same event window.

[0122] : The amplitude index of the i-th channel within the event window Ω.

[0123] c i : Consistency correction coefficient for the i-th channel.

[0124] A0: An unknown scale factor strongly correlated with the event source.

[0125] α: Attenuation coefficient.

[0126] γ: Geometric diffusion index.

[0127] C: Set of channel pairs .

[0128] Channel pair weights.

[0129] : Channel pair residual function.

[0130] ρ( ): Robust loss function.

[0131] V: Feasible spatial domain of the seismic source.

[0132] RMS: Root Mean Square of Residuals.

[0133] SNR i : Signal-to-noise ratio of the i-th channel.

[0134] τ: Outlier detection threshold.

[0135] MAD(): Median absolute deviation.

[0136] β,λ: Threshold coefficient / Update coefficient.

[0137] M: Rolling update of sample count.

[0138] A three-dimensional source localization method based on amplitude attenuation consists of the following steps:

[0139] S1: Acquire waveform data of the same event on at least four sensor channels. And obtain the spatial coordinates of each sensor. When the number of sensor channels is 4, the spatial layout of the sensors should meet the geometric non-degenerate condition (e.g., the four points are not coplanar or the volume / rank formed by the sensor coordinates is greater than a preset threshold); if the geometric non-degenerate condition is not met, the number of sensor channels should be increased, or a feasible domain constraint of the source space should be introduced and an uncertainty prompt should be output.

[0140] S2: For waveform data of each channel Perform preprocessing operations. The preprocessing operations include one or more of the following: DC removal, trend removal, bandpass filtering, and noise reduction.

[0141] S3: Determine the event window Ω corresponding to the event. The event window Ω can be manually specified or automatically determined through threshold triggering, short-to-long window energy ratio, energy mutation, etc. To ensure the comparability of amplitude indicators extracted from different channels, time alignment can be performed on the event windows of each channel after the event window is determined. This time alignment includes energy peak alignment, cross-correlation alignment, or a combination of both. Preferably, envelope peaks or energy indicators are extracted within the same passband, so that the amplitude indicators correspond to the propagation response of the same wave packet / same frequency band.

[0142] S4: Extract the amplitude index A of each channel within the event window Ω. i The amplitude index is then subjected to channel consistency correction to obtain the corrected amplitude:

[0143] ,

[0144] .

[0145] Differences in sensor sensitivity, coupling conditions, and gain settings across different channels can lead to incomparable amplitudes across different channels for the same event. To ensure that the amplitude ratio reflects the propagation attenuation pattern, consistency correction is performed on the channels after amplitude extraction. Correction may include:

[0146] (1) Coefficient determination based on calibration events

[0147] Using a set of calibrated events with known locations For each event Calculate its geometric distance r to each channel. i,e Under the exponential decay model, the expected outcome is:

[0148] ,

[0149] With a fixed reference channel k, let c k =1. For each channel i ≠ k, estimate c by minimizing the following equation. i :

[0150] ,

[0151] To avoid the impact of abnormal events, median statistics or robust regression can be used for the event set.

[0152] (2) Reference normalization under calibration conditions

[0153] When calibration events are lacking, relative gain can be obtained by statistically analyzing multiple events within a time window. With a fixed reference channel k, for each channel i ≠ k, the following definition is provided:

[0154] ,

[0155] Even in multiple events A i With A k This method achieves the same scale in a statistical sense. It is independent of the epicenter location but requires the number of events to meet a preset lower limit E. min (For example, no fewer than 20).

[0156] S5: Determine the medium attenuation coefficient α. α is determined by calibration experiments, historical data inversion, or regional medium parameters.

[0157] ,

[0158] Where f is the main frequency or center frequency, v is the propagation speed, and Q is the quality factor.

[0159] S6: When the decay model contains only an exponential decay term, the exponential decay model is used.

[0160] ,

[0161] Where A0 is the source strength factor and α is the attenuation coefficient. Construct the logarithmic ratio observation for any channel pair (i,j):

[0162] ,

[0163] Substituting into the attenuation model, we get

[0164] ,

[0165] Thus, the distance difference observation is obtained:

[0166] ,

[0167] Establish the location of the earthquake source To each sensor S i Geometric distance relationship r i Constructing a channel And establish the distance difference residual:

[0168] ,

[0169] When the attenuation model further includes a geometric diffusion term When doing so, both the distance difference-related term and the distance ratio-related term should be considered in the residual of the objective function to ensure that the logarithmic ratio observation is consistent with the attenuation model, rather than directly dividing the logarithmic ratio by α to obtain the distance difference. In this case, the amplitude index is:

[0170] ,

[0171] Where A0 is the source strength factor, r i Let be the propagation distance from the earthquake source to the i-th sensor, γ be the geometrical diffusion exponent, and α be the attenuation coefficient. The logarithmic ratio can be written as:

[0172] ,

[0173] Based on this, a hybrid residual (constraining both distance difference and distance ratio) is constructed:

[0174] ,

[0175] S7: Establish earthquake source location Geometric distance relationship with each sensor location:

[0176] ,

[0177] S8: Combine the distance difference constraint obtained in step S6 with the geometric constraint obtained in step S7, and select the channel pair set C( Construct a robust objective function and solve for the location of the earthquake source:

[0178] ,

[0179] Or adopt with The corresponding absolute value residual form, where ρ(⋅) is the weight, and ρ(⋅) is the robust loss function. The threshold can be adaptively determined based on the background noise statistic or the robust residual statistic, and can be updated with iteration.

[0180] ,

[0181] SNR i Let be the signal-to-noise ratio of the i-th channel. A positive constant is preset to prevent the denominator from being zero.

[0182] During the iterative solution process, an outlier detection threshold is adaptively determined based on the background noise statistic or the residual robust statistic, and is used to update the channel pair weights or remove abnormal channel pairs; the threshold can be updated with iteration.

[0183] (1) The residual threshold τ is determined based on the statistics of the background noise segment outside the event window, satisfying:

[0184] ,

[0185] Where σ is the standard deviation of the background noise segment, and β is the preset multiplier coefficient;

[0186] When the residual threshold τ is determined based on robust statistics of the residuals, the following conditions are met:

[0187] ,

[0188] in β is a preset coefficient.

[0189] (2) The residual threshold τ is updated during the iterative solution process, satisfying:

[0190] ,

[0191] in λ is the threshold calculated based on the residual statistics in the l-th iteration, and λ∈(0,1] is the update coefficient.

[0192] When using absolute value constraints or absolute value residuals, multiple candidate solutions are obtained by solving within the preset feasible region of the seismic source space through multiple initial values ​​or coarse grid initial values. The final solution is selected based on the criterion of minimizing the objective function value and satisfying the feasible region constraints. When the objective function values ​​of multiple candidate solutions are close, the confidence level or uncertainty index is output.

[0193] S9: Initial value P (0) (The monitoring area center) begins iteration until any of the following stopping conditions are met:

[0194] (1) ,

[0195] (2) ,

[0196] (3) The number of iterations reaches the upper limit L max .

[0197] S10: Output the final estimate And positioning quality indicators such as root mean square residuals (RMS):

[0198] ,

[0199] And reselect the channel pair set when the residual exceeds the threshold. Or adjust the weight It is updated iteratively to achieve automatic outlier constraint removal and reweighting. Specific Implementation Example 1:

[0201] The method of this invention will be explained below using four sensor channels 1, 2, 3, and 4 as an example. The maximum amplitude is constructed through four channels. The distance difference set is used to solve the three-dimensional coordinates of the earthquake source using nonlinear least squares.

[0202] S1: Preparation and Data Acquisition

[0203] (1) Obtain the spatial coordinates of the four sensors S1, S2, S3, and S4. When the number of sensors is 4, in order to ensure the solvability and stability of the three-dimensional positioning, the four sensors should preferably be arranged in a non-coplanar spatial configuration to avoid instability of the solution caused by geometric degradation. If the configuration conditions limit the coplanar configuration, the feasible domain constraint of the seismic source space can be introduced in the solution stage to improve the stability and uniqueness of the solution. The coordinates can come from construction surveys, layout maps, or spatial calibration results; the coordinate system can be the local coordinate system of the mine / tunnel or the geographic coordinate system, and the unit is meters.

[0204] (2) Set the sampling rate F s(Engineering available range of 500Hz–20kHz), and ensure four-channel acquisition time synchronization (multi-channel synchronization on the same acquisition unit or synchronization using a unified clock).

[0205] (3) Acquire the four-channel waveforms W1(t)~W4(t) corresponding to the same event. Optionally, the background noise segment can be acquired at the same time for subsequent signal-to-noise ratio estimation.

[0206] S2: Waveform Preprocessing and Quality Control

[0207] (1) Perform DC removal, trend removal and bandpass filtering on each channel; the bandpass frequency can be set around the main event frequency (10–500Hz or 100–3000Hz, depending on the application scenario).

[0208] (2) If the acquisition system has an amplitude limit, detect whether the waveform is clipped / saturated; if saturated, mark the channel as low confidence and assign it a lower weight or remove it in the future.

[0209] (3) Calculate the root mean square noise N in the background noise segment. i Calculate the signal peak value or energy S within the event window. i The signal-to-noise ratio is obtained as follows:

[0210]

[0211] When SNR i When the SNR is below a threshold (e.g., 3–10), the channel can be marked as low confidence. The processing rules for low confidence channels are as follows: when a channel is determined to be saturated, has an SNR below the threshold, or has obvious anomalies (such as missing data, abrupt changes, breakpoints, or strong interference peaks dominating), the channel can be marked as invalid, and channel pairs related to this channel will be removed during subsequent channel pair construction, or the weight of the related channel pairs will be set to 0; if it is only low confidence (not invalid), then the related channel pairs will be given a smaller weight in the subsequent objective function to reduce their impact on the localization results.

[0212] S3: Event Detection and Event Window Determination

[0213] (1) The event start time t0 is determined by threshold triggering, short-to-long window energy ratio (STA / LTA) or energy mutation detection.

[0214] (2) Set the event window based on t0. , where t p For the reserved lead time, t s For example, t is the post-time; p Take 5–50 ms, t s The time interval is set to 50–1000 ms, and can be adjusted adaptively based on the propagation distance and main frequency.

[0215] (3) If there is a deviation in the triggering time of the four channels, the start and end boundaries of Ω can be fine-tuned by energy peak alignment or cross-correlation alignment (fine-tuning occurs within the window and is used for subsequent amplitude extraction consistency) to improve the consistency of maximum amplitude extraction.

[0216] Energy peak alignment: Calculate the energy sequence or envelope sequence of each channel within the event window, and use the time corresponding to the energy peak or envelope peak as the alignment reference.

[0217] Cross-correlation alignment: Calculate the cross-correlation function of each channel and the reference channel within the event window, and use the time delay corresponding to the peak cross-correlation value as the alignment amount.

[0218] S4: Calculate the maximum amplitude within the event window Ω:

[0219] ,

[0220] To reduce the impact of instantaneous spike noise, the analytical envelope can be calculated first:

[0221] ,

[0222] In the formula, H{} is the Hilbert transform, and then take .

[0223] For cases where different amplitude definition methods are used, distance difference constraints can be constructed and the seismic source location can be completed by "logarithmic elimination of A0 by amplitude ratio".

[0224] Differences in sensor sensitivity, coupling conditions, and gain settings across different channels can lead to incomparable amplitudes across different channels for the same event. To ensure that the amplitude ratio reflects the propagation attenuation pattern, consistency correction is performed on the channels after amplitude extraction. Correction may include:

[0225] (1) Coefficient determination based on calibration events

[0226] Using a set of calibrated events with known locations For each event Calculate its geometric distance r to each channel. i,e Under the exponential decay model, the expected outcome is:

[0227] ,

[0228] With a fixed reference channel k, let c k =1. For each channel i ≠ k, estimate c by minimizing the following equation. i :

[0229] ,

[0230] To avoid the impact of abnormal events, median statistics or robust regression can be used for the event set.

[0231] (2) Reference normalization under calibration conditions

[0232] When calibration events are lacking, relative gain can be obtained by statistically analyzing multiple events within a time window. With a fixed reference channel k, for each channel i ≠ k, the following definition is provided:

[0233]

[0234] Even in multiple events A i With A k This method achieves the same scale in a statistical sense. It is independent of the epicenter location but requires the number of events to meet a preset lower limit E. min (For example, no fewer than 20).

[0235] S5: Determination of the attenuation coefficient α

[0236] (1) Calculated from the dominant frequency or center frequency f, propagation speed v, and quality factor Q:

[0237] ,

[0238] (2) Select a calibration event at a known location (such as a burst point / calibration source) P cal Given P cal Calculate the geometric distance from each sensor to the calibration source under the given conditions:

[0239] ,

[0240] Extract the maximum amplitude of each channel within the event window of the calibration event. Based on the exponential decay model:

[0241] ,

[0242] Construct the following for any channel (i,j):

[0243] ,

[0244] Thus, after eliminating the unknown source strength factor A0, α is obtained by linear fitting or least-squares fitting of the data across multiple channels. Optionally, the above process is repeated using multiple calibration events, and the fitted α is statistically averaged or robustly estimated.

[0245] (3) Divide the monitoring space into several medium blocks or depth intervals, and pre-set α candidate values ​​or candidate ranges for each block; during the solution process, select the corresponding α according to the initial value of the source location or the block to which the iteration location belongs. Furthermore, α can be reselected from the candidate set according to the statistics of the distance difference residual (e.g., residual RMS or robust cost value) to improve the consistency between the model and the actual propagation attenuation.

[0246] (4) When the positioning residual is systematically large in multiple consecutive events and no obvious abnormalities are found in the channel quality control, the correction process of α can be triggered: preferably, select the calibration event in the most recent time window, repeat the "logarithmic fitting of amplitude ratio" in step (2) to update α, or perform rolling updates of α based on the distance difference residual statistics. Optionally, α and the source position P can be used as joint parameters to be estimated and jointly optimized and updated in the same objective function.

[0247] S6: Explicit mapping and distance difference calculation

[0248] Construct a set of channel pairs for the four-channel case:

[0249] ,

[0250] And define the set of absolute values ​​of distance differences:

[0251] ,

[0252] ,

[0253] Corresponding to calculation based on amplitude:

[0254] ,

[0255] ,

[0256] like and If there is a statistically significant inconsistency (exceeding a threshold), it can be determined that there is an abnormal channel or abnormal amplitude, triggering a reweighting / removal process. The "statistically significant inconsistency" can be determined by residual exceeding a threshold, for example, when some channels continuously satisfy the residual... hour( For a single pair of residual thresholds, the corresponding weight is reduced or the channel pair is removed and the solution is recalculated based on noise level, sensor layout geometry, sampling rate and engineering experience (or determined adaptively through historical event statistics).

[0257] S7: Establish geometric constraints and construct the objective function for solving the problem.

[0258] Establish geometric distance relationships:

[0259] ,

[0260] Construct the absolute value residual objective function (and) correspond):

[0261] ,

[0262] in for Corresponding item, weight It can be given in the following manner:

[0263] ,

[0264] SNR i Let be the signal-to-noise ratio of the i-th channel. A positive constant is preset to prevent the denominator from being zero.

[0265] Or associate saturated / coupled poor channels Set to a smaller value. Weight It can be determined based on signal-to-noise ratio, coupling quality, or amplitude stability; if any channel is saturated or ineffective, then the corresponding... .

[0266] S8: Determining the location of the earthquake source

[0267] Take the geometric center of the sensor as the initial value P (0) Alternatively, a coarse-grid search can be performed within the monitoring volume to select the point that minimizes the residual as the initial value. Then, a nonlinear least squares algorithm such as Levenberg-Marquardt is used to iteratively update P. (l) →P (l +1) .when Alternatively, iteration can stop when the objective function decreases below a threshold; the maximum number of iterations can be set to 20–100. If a channel consistently shows a significantly larger residual, the corresponding... Decrease or remove the constraint, and then solve the problem again.

[0268] S9: Outputs the three-dimensional coordinates (x, y, z) of the seismic source and outputs the positioning quality indicators for subsequent automatic quality inspection or manual verification. Specific Implementation Example 2:

[0270] S1: Set up 4 sensors with the following coordinates (where S4 has a significant height difference in the Z direction to form a non-coplanar geometry. When the number of sensors is 4, in order to ensure the solvability and stability of the 3D positioning, the 4 sensors should form a non-coplanar spatial layout): S1 (0,0,0); S2 (30,0,0); S3 (0,25,0); S4 (10,10,25) m.

[0271] S2: Select a calibration source point at a known location (e.g., artificially triggered or known detonation point): P cal =(16,10,6)m.

[0272] S3: Extract the peak amplitudes A1, A2, A3, and A4 of the four channels from the calibration event. Based on the exponential decay model.

[0273] ,

[0274] Eliminating A0 from channel pair (i,j) yields:

[0275] ,

[0276] By performing least-squares fitting on the constraints using multiple sets of calibration events or multiple channels, the attenuation coefficient in this embodiment is obtained: α = 0.035m -1 .

[0277] S4: For this calibration event, obtain the peak amplitude of the four channels:

[0278] A1=1.0071, A2=1.0740, A3=0.9143, A4=1.3339.

[0279] S5: Construct distance difference constraints based on amplitude ratio:

[0280] ,

[0281] And establish residuals:

[0282] ,

[0283] The 3D coordinates P are solved using weighted least squares or LM iteration. Using the sensor's geometric center as the initial value, the positioning result is obtained:

[0284] m,

[0285] The three-dimensional error relative to the known calibration point is:

[0286] m.

Claims

1. A three-dimensional source localization method based on multi-channel amplitude ratio, characterized in that: Includes the following steps: a. Acquire waveform data of the same event on at least four sensor channels. and the spatial coordinates of at least four sensors. ; b. Perform preprocessing on each channel waveform and determine the same event window Ω; within the event window Ω, perform multi-channel alignment based on the reference channel k to ensure that the amplitude extraction of different channels corresponds to the same event wave packet; c. Extract the amplitude index A of each channel within the event window Ω. i And for amplitude index A i Perform channel consistency correction to obtain the corrected amplitude : ， ， Where c i This is the consistency correction coefficient; d. Based on the exponential decay model Constructing logarithmic ratio observations and establishing residuals, firstly, for any two channels, the amplitude index A... i A j Constructing logarithmic ratio observations : ， Based on the exponential decay model, the amplitude indices of channel i and channel j are compared and logarithmically transformed to eliminate the unknown source strength factor A0 of the event, thus obtaining the distance difference observable. : ， Where r i Let the propagation distance from the earthquake source to the i-th sensor be: P is the location of the earthquake source to be determined, and α is the attenuation coefficient; When the exponential decay model further includes a geometric diffusion term When γ is the preset or calibrated geometric diffusion index; the residual term consists of a "term related to the distance difference" and a "term related to the distance ratio" to ensure that the logarithmic ratio observation is consistent with the attenuation model, at which point the maximum amplitude is... Determine by the following formula: ， Where A0 is the source strength factor, r i γ is the propagation distance from the earthquake source to the i-th sensor, γ is the geometric diffusion index, and α is the attenuation coefficient; e. Construct a redundant channel constraint set And establish the distance difference residual: ， Channel pairs are screened, downweighted, or eliminated based on the consistency of residuals between them as a quality control criterion. f. After screening, weighting or elimination, establish the geometric distance relationship for the source location P and construct the objective function. Solve the source location through nonlinear least squares to minimize the objective function. ， in Let ρ be the channel pair weights, and ρ() be the robust loss function; The solution employs iterative reweighted least squares or introduces a robust loss function to reduce the impact of outlier channels on the solution results; g. The residual threshold is adaptively determined based on the noise statistics or residual robust statistics estimated from the background noise segment, and updated it iteratively to achieve automatic weight reduction or elimination of outlier constraints; h. Output the estimated seismic source location. And positioning quality indicators, which should include at least the minimum value of the objective function or the root mean square of the residuals.

2. The three-dimensional source localization method based on multi-channel amplitude ratio according to claim 1, characterized in that: The amplitude index A of the i-th channel i Determine as follows: Take the absolute peak value of the channel waveform within the event window Ω: ， Alternatively, perform a Hilbert transform on the waveform to obtain the envelope, and then take the envelope peak value. ， Where H{} is the Hilbert transform, This represents the waveform of the i-th channel.

3. The three-dimensional source localization method based on multi-channel amplitude ratio according to claim 1, characterized in that: When the distance difference observation in step d is expressed in absolute value form: ， Furthermore, step f uses the corresponding absolute value residual: 。 4. The three-dimensional source localization method based on multi-channel amplitude ratio according to claim 1, characterized in that: The attenuation coefficient α is determined in at least one of the following ways: Based on the dominant frequency or center frequency f, the propagation speed v, and the quality factor Q, the following is calculated: ， Alternatively, the location can be determined by calibration events at known positions. Select several calibration events whose locations are known or whose locations can be determined with high precision by other methods, and denote their distance to the i-th sensor as r. i The amplitude is A i In the model without a geometric diffusion term, we have: ， For lnA i With r i Linear regression can fit α and simultaneously obtain lnA0. To reduce the influence of outliers, robust regression is used. If a model with a geometric diffusion term is used, the fitted equation is: ， Thus, α and γ are jointly fitted in the same regression; The attenuation coefficient α is a partitioning parameter or a layering parameter and is determined using a rolling update method. When the medium in the monitoring area is non-uniform, the area is divided into several sub-regions R. m α is determined for each sub-region. m When a new event arrives, if its localization residual is consistently large, then the corresponding sub-region's α is updated using this event and the M most recent high-quality events as samples. m Updated to: ， in The estimate is obtained based on regression of the new sample, and λ∈(0,1] is the update coefficient; The attenuation coefficient α is jointly estimated with the source location P. When α is uncertain and sensitive to location, α will be used as a parameter to be estimated and inverted together with P: 。 5. The three-dimensional source localization method based on multi-channel amplitude ratio according to claim 1, characterized in that: The channel consistency correction addresses the issue that differences in sensor sensitivity, coupling conditions, and gain settings across different channels can lead to incomparable amplitudes across different channels under the same event. To ensure that the amplitude ratio reflects the propagation attenuation law, consistency correction is performed on the channels after amplitude extraction. Let the observed amplitude of channel i be A. i The consistency correction factor is c. i The corrected amplitude is defined as: ， Correction factor c i Used to compensate for differences in channel sensitivity, so that the amplitudes of different channels are in the same dimension and scale under the same propagation conditions; (1) Coefficient determination based on calibration events Using a set of calibrated events with known locations For each event Calculate its geometric distance r to each channel. i,e Under the exponential decay model, the expected outcome is: ， Where r i,e Let be the geometric distance from the source location of event e to the i-th sensor. Let e ​​be the source strength factor. Let e ​​be the corrected amplitude index for event e on the i-th channel; With a fixed reference channel k, let c k =1, for each channel i ≠ k, estimate c by minimizing the following equation. i : ， in Let r be the amplitude index extracted from event e on the i-th channel. k,e Let be the geometric distance from the epicenter location of event e to the sensor in reference channel k. To avoid the impact of abnormal events, median statistics or robust regression can be used for the event set; (2) Reference normalization under calibration conditions When calibration events are lacking, the relative gain is obtained by statistically analyzing multiple events within a time window. With a fixed reference channel k, the following definition is used for each channel i ≠ k: ， Where exp() is an exponential function with base e. To calculate the median for event index e; Even in multiple events A i With A k To achieve the same statistical scale, this method is independent of the epicenter location, but requires the number of events to meet a preset lower limit E. min .

6. The three-dimensional source localization method based on multi-channel amplitude ratio according to claim 1, characterized in that: The weights in the weighted nonlinear least squares Determined according to the following formula: ， SNR i Let be the signal-to-noise ratio of the i-th channel. A positive constant is preset to prevent the denominator from being zero.

7. The three-dimensional source localization method based on multi-channel amplitude ratio according to claim 1, characterized in that: In the solution process, outlier channels or constraints can be removed, or a robust loss function can be used to reduce the impact of outliers. ， Where ρ() is the robust loss function, and further, it can be implemented using iterative reweighted least squares (IRLS): in the l-th iteration, with ， Update the equivalent weights, where ω(⋅) is the weight function derived from the selected robust loss.

8. The three-dimensional source localization method based on multi-channel amplitude ratio according to claim 1, characterized in that: The root mean square of the residuals used for anomaly detection is determined as follows: 。 9. The three-dimensional source localization method based on multi-channel amplitude ratio according to claim 1, characterized in that: When the number of sensor channels is four, the spatial arrangement of the four sensors satisfies the geometric non-degradation condition. The geometric non-degradation condition includes: the four points are not coplanar or the volume of the geometric body formed by the coordinates of the sensors is greater than a preset threshold. If the geometric non-degradation condition is not met, the number of sensor channels is increased, or the feasible domain constraint of the seismic source space is introduced and an uncertainty prompt is output. When the sensor channels are four channels 1, 2, 3, and 4, the constructed distance difference set includes at least one of the following: ， ， And respectively, the amplitude was calculated as follows: ， 。 10. The three-dimensional source localization method based on multi-channel amplitude ratio according to claim 1, characterized in that: In step g, the residual threshold τ is determined based on the statistics of the background noise segment outside the event window, satisfying the following: ， Where σ is the standard deviation of the background noise segment, and β is the preset multiplier coefficient; When the residual threshold τ is determined based on robust statistics of the residuals, the following conditions are met: ， in β is a preset coefficient; In step g, the residual threshold τ is updated during the iterative solution process, satisfying: ， in λ is the threshold calculated based on the residual statistics in the l-th iteration, and λ∈(0,1] is the update coefficient.

Images