Audio magnetotelluric data quality rapid evaluation method

Abstract

The present application relates to the field of digital processing, and discloses a method for quickly evaluating the quality of audio magnetotelluric data, comprising: reading original observation data of audio magnetotelluric measuring points stored in a hard disk or memory; obtaining a sorted frequency sequence and a sorted response parameter sequence; performing data domain transformation; first weighting, second weighting and third weighting; performing second weighting spline fitting using the third weight to obtain a final reference curve; reading the data to be evaluated stored in the hard disk or memory, comparing the data to be evaluated with the reference curve point by point, calculating the comprehensive similarity of each effective frequency point, judging the qualified frequency points according to the comprehensive similarity, counting the proportion of the qualified frequency points, calculating the comprehensive score, dividing the quality level and outputting the evaluation result. The present application only uses the curvature characteristics of the original data itself, prior band knowledge and residual statistical characteristics to construct a reference curve which objectively reflects the main trend of the data and is suitable for any electrical structure.

Description

Technical Field

[0001] This invention relates to the field of digital electrical processing, and more particularly to a method for rapid evaluation of the quality of audio magnetotelluric data. Background Technology

[0002] Audio-frequency magnetotellurics (AMT) uses natural alternating electromagnetic fields as its source and inverts the resistivity distribution of subsurface media by observing the mutually orthogonal electric and magnetic field components at the Earth's surface. Compared to traditional electrical exploration methods, AMT has significant advantages such as greater detection depth, higher resolution, lighter equipment, and immunity to high-resistivity shielding, making it an important tool for shallow and intermediate-depth geophysical exploration.

[0003] With the widespread adoption of analog-to-mass measurement (AMT) in mineral resource exploration, geothermal surveys, and engineering geology, the number of measurement points collected in a single exploration mission has increased dramatically, resulting in massive amounts of measured data. However, the natural source signals upon which AMT relies are extremely weak, making the acquisition process highly susceptible to factors such as power frequency interference, dead-band effects, static distortion, and human noise, leading to inconsistent data quality. The quality of the data directly determines the reliability of subsequent impedance estimation, inversion imaging, and geological interpretation. Therefore, how to conduct rapid, objective, and automated quality assessment of massive AMT data has become a critical issue that urgently needs to be addressed in current exploration data processing workflows.

[0004] Currently, the evaluation of AMT data quality both domestically and internationally mainly employs the following methods: (1) Visual identification method: According to industry standards, technicians manually observe the continuity and smoothness of the apparent resistivity and phase curves and the proportion of "qualified frequency points" across the entire frequency band to classify data quality into different levels. This method is essentially a semi-quantitative evaluation based on "statistics + expert experience", which is highly dependent on subjective judgment, inefficient, and cannot meet the batch processing needs of large data volumes.

[0005] (2) Comparison method based on external reference curve: Select a known high-quality reference curve (such as the curve of nearby high-quality measuring points, the forward model curve of theoretical model or regional statistical curve), and compare the data to be evaluated with the reference curve point by point. The quality is quantified by the degree of deviation of apparent resistivity and phase. The accuracy of this method is highly dependent on the reliability of the reference curve. However, the actual underground structure is complex and it is difficult to obtain a universally applicable standard curve.

[0006] (3) Full tensor fitting evaluation method: This method does not rely on external reference curves, but instead fits an optimal one-dimensional anisotropic medium model using the test data itself, and uses the theoretical response of this model as an "internal reference curve". This method works well when the data quality is good and the electrical structure is close to one-dimensional / two-dimensional anisotropy, but when the underground structure is complex (such as three-dimensional, with static effects), the fitted reference curve itself may deviate from the true response, causing evaluation bias.

[0007] When applied to audio magnetotelluric data, the above method has the following specific limitations: (1) The non-uniform distribution of interference in the AMT band is not fully considered: The AMT frequency range is usually 1Hz to 10kHz, with obvious "dead bands" (the natural field energy is extremely low in the range of 1000 to 5000Hz, and the curve jumps and is discontinuous), power frequency harmonic interference (50Hz and its harmonics), static effects caused by shallow inhomogeneities (the apparent resistivity curve is shifted as a whole, and the phase is normal), and high-frequency statistical jitter caused by short observation time. Existing methods mostly treat the entire frequency band equally, which may lead to dense jump points in the dead band incorrectly dominating the scoring, or overfitting in the fitting process.

[0008] (2) Model dependency problem in reference curve acquisition: Both external standard curves and full tensor fitting curves implicitly assume prior assumptions about underground electrical structure. For data affected by static effects, superphase, etc., the fitted reference curve may be "skewed", thus making it impossible to objectively evaluate the degree of agreement between the original data and the actual geological response.

[0009] (3) Existing smoothing methods lack targeted weight design: Conventional spline smoothing or median filtering uses equal weight or simple distance weights, which cannot distinguish between high-frequency noise and curve changes caused by the real electrical interface, and are prone to over-smoothing or retaining noise.

[0010] (4) Low evaluation efficiency and inability to provide real-time feedback: Existing methods that require iterative inversion or manual intervention are time-consuming, making it difficult to achieve real-time quality monitoring at the field collection site, and also unable to support large-scale batch rapid screening of thousands of measurement points.

[0011] Therefore, there is an urgent need for a reference curve construction method that does not rely on underground model assumptions, can adapt to the interference characteristics of each frequency band of AMT, and is computationally efficient, so as to achieve objective and rapid quality evaluation. Summary of the Invention

[0012] The purpose of this invention is to overcome the shortcomings of the prior art and provide a rapid evaluation method for audio magnetotelluric data quality.

[0013] The objective of this invention is achieved through the following technical solution: A first aspect of the present invention provides a method for rapid evaluation of the quality of audio magnetotelluric data, comprising the following steps: Read the raw observation data of the audio magnetotelluric measurement points stored on the hard disk or in memory. The raw observation data includes a frequency array and a response parameter array. The response parameter array includes an apparent resistivity array or an impedance phase array. The frequency array and the response parameter array are preprocessed to convert them into column vectors and sorted in ascending order of frequency to obtain sorted frequency sequences and sorted response parameter sequences. Perform a logarithmic transformation on the sorted frequency array to generate a logarithmic frequency sequence, which serves as the sorted frequency array after data domain transformation. ; and perform data domain transformation on the sorted response parameter sequence according to the type of the response parameter array: when the response parameter array is an apparent resistivity array, perform a logarithmic transformation on the sorted response parameter sequence to generate a logarithmic response sequence as the sorted and data domain transformed response parameter array. When the response parameter array is an impedance phase array, the units of the impedance phases in the sorted response parameter sequence are converted to radians, thus obtaining the sorted and data domain-transformed response parameter array. , The array index; The first weighting is a cubic spline fit based on geometric weights of the second derivative: calculating the frequency array after sorting and data domain transformation. With the sorted and domain-transformed response parameter array Approximate value of the second derivative Approximating the second derivative The first weight is calculated by substituting the first linear truncation weight function based on the threshold, and the approximate value of the second derivative in the first linear truncation weight function is obtained. The larger the point, the smaller the initial weight; the second derivative approximation. The smaller the point, the greater its initial weight. The threshold in the first linear cutoff weight function is calculated using an approximation of the second derivative. Multiply the standard deviation by a preset factor; substitute the first weight into the weighted cubic smooth spline function to perform the first weighted spline fitting, and obtain the first fitted curve and the first fitted value of each point; The second weighting is a bias weight based on the fitting residuals: the first fitted values ​​are compared with the sorted and domain-transformed response parameter array. residual The second weight is calculated using a threshold-based second linear truncation weight function. In this function, points with smaller residuals have larger second weights, and points with larger residuals have smaller second weights. The threshold in the second linear truncation weight function is calculated using the residuals... The average or residual the median; The third weighting is a forced weight reduction of dead frequency bands based on prior knowledge: the frequency index corresponding to the dead frequency band of audio magnetotelluric data is identified based on prior knowledge, and the frequency points in the dead frequency band are forcibly set with minimal weights, while the weights of the other frequency points remain unchanged from the second weight. The overall weight is used as the third weight. A second weighted spline fitting is performed using the third weight to obtain the final reference curve; Read the data to be evaluated stored in the hard disk or memory, compare the data to be evaluated with the reference curve point by point, calculate the comprehensive similarity of each effective frequency point, determine the qualified frequency points based on the comprehensive similarity, count the proportion of qualified frequency points, calculate the comprehensive score, divide the quality level based on the comprehensive score and output the rating result.

[0014] Furthermore, the first linear truncation weight function based on the threshold is: When the curvature threshold is zero, the first weight is one; when the absolute value of curvature is between zero and the curvature threshold, the first weight decreases linearly; when the absolute value of curvature is greater than or equal to the curvature threshold, the first weight is truncated to zero. The curvature threshold is adaptively determined based on the overall curvature distribution of the data, that is, the curvature threshold is equal to a preset multiple multiplied by the standard deviation of the absolute value of curvature.

[0015] Furthermore, the threshold-based second linear truncation weight function is: The second weight is one when the residual is zero; the second weight decreases linearly when the residual is between zero and the residual mean; and the second weight is truncated to zero when the residual is greater than or equal to the residual mean. The residual mean is the arithmetic mean of the residuals.

[0016] Furthermore, the step of identifying the frequency point index corresponding to the dead frequency band of audio magnetotelluric data based on prior knowledge specifically involves: The frequency range of the dead band is determined based on the collected parameters, and the frequency range of the dead band is from 1000 Hz to 5000 Hz; the set of frequency point indices falling within this frequency range is determined on the logarithmic frequency axis. The minimum weight is 0.001.

[0017] Furthermore, the step of comparing the data to be evaluated with the reference curve frequency by frequency includes a frequency alignment step: Take the base-10 logarithm of the frequency arrays of the data to be evaluated and the reference data respectively; For each frequency point in the data to be evaluated, the nearest frequency point is found in the logarithmic frequency sequence of the reference data. If the minimum distance is less than the preset tolerance, the match is accepted; otherwise, the frequency point to be evaluated is discarded, and the effective frequency point set and the corresponding reference data subset are obtained. The preset tolerance is 0.02 logarithmic units.

[0018] Furthermore, the overall similarity is calculated in the following manner: Apparent resistivity similarity and phase similarity are calculated separately, and then the apparent resistivity similarity and phase similarity are fused together with equal weight to form a comprehensive similarity. The apparent resistivity similarity is obtained by linearly mapping the relative error in the logarithmic domain to the interval of 0 to 1, and the phase similarity is obtained by linearly mapping the absolute value of the phase difference to the interval of 0 to 1 after considering the 360-degree periodic rotation. In the process of determining qualified frequency points based on comprehensive similarity, the qualified threshold for comprehensive similarity is 0.8.

[0019] Furthermore, the comprehensive score is calculated using a weighted scoring method: Weights are assigned based on the distribution density of valid frequency points on the logarithmic frequency axis. The weight of each frequency point is determined according to the width of the frequency range represented by that frequency point on the logarithmic frequency axis. The endpoints are spaced on one side, and the intermediate frequency points are spaced by half the distance between adjacent frequency points. After normalizing the weights of each frequency point, the normalized weights are multiplied by the qualified mark of that frequency point and summed to obtain the weighted comprehensive score.

[0020] Furthermore, the quality grades include four levels: excellent, good, medium, and unqualified. The classification rules are as follows: an overall score greater than or equal to 0.95 is rated as excellent; an overall score less than 0.95 but greater than or equal to 0.80 is rated as good; an overall score less than 0.80 but greater than or equal to 0.60 is rated as medium; and an overall score less than 0.60 is rated as unqualified.

[0021] Furthermore, the smoothing parameter for the first weighted spline fitting is 0.985, and the smoothing parameter for the second weighted spline fitting is 0.955.

[0022] The beneficial effects of this invention are: In an exemplary embodiment of the present invention, a weighted spline fitting framework consisting of curvature weight, residual weight, and dead-band forced weighting is proposed. Each weight has a clear physical meaning and works synergistically: First weighting (curvature weight): Utilizing the geometric features of the curve, local high-frequency spikes and static effect peaks are automatically suppressed, preserving smooth macroscopic electrical interface information; Second weighting (residual weight): Based on the residual after the previous fitting, outlier interference points are further eliminated, making the reference curve more faithful to the main electrical trend and avoiding overfitting; Third weighting (dead-band forced weighting): For dead-band interference unique to audio magnetotellurics, the corresponding frequency weight is reset to a minimum value (e.g., 0.001) based on prior knowledge, so that the fitted curve smoothly passes through this region without being skewed by jump points.

[0023] Furthermore, it requires no assumptions about underground models and is data-driven. This exemplary embodiment utilizes only the curvature characteristics, prior frequency band knowledge, and residual statistical characteristics of the original data itself, without introducing any prior geological models. The constructed reference curve objectively reflects the main trend of the data and is applicable to any electrical structure (one-dimensional, two-dimensional, three-dimensional, including static effects, etc.). Attached Figure Description

[0024] Figure 1 A flowchart of a rapid evaluation method for audio magnetotelluric data quality provided as an exemplary embodiment of the present invention; Figure 2 A comparison diagram of the original apparent resistivity curve of noisy synthetic data provided for an exemplary embodiment of the present invention and reference curves obtained by different methods; Figure 3 A comparison diagram of the original phase curve of noisy synthetic data and reference curves obtained by different methods, provided for an exemplary embodiment of the present invention; Figure 4 In an exemplary embodiment of the present invention, the original data curve is the apparent resistivity curve, and the A-level AMT measurement point curve of the constructed reference data is provided in the XY mode. Figure 5 In the YX mode provided as an exemplary embodiment of the present invention, the original data curve is the apparent resistivity curve, and the A-level AMT measurement point curve of the constructed reference data is plotted. Figure 6 In an exemplary embodiment of the present invention, an A-level AMT measurement point curve diagram of the original data curve as the phase curve and the constructed reference data is provided in the XY mode. Figure 7 In the YX mode provided as an exemplary embodiment of the present invention, the original data curve is a phase curve, and an A-level AMT measurement point curve diagram of the constructed reference data is shown. Figure 8 In an exemplary embodiment of the present invention, the original data curve is the apparent resistivity curve, and the B-level AMT measurement point curve of the constructed reference data is provided in the XY mode. Figure 9 In the YX mode provided as an exemplary embodiment of the present invention, the original data curve is the apparent resistivity curve, and the B-level AMT measurement point curve diagram of the constructed reference data is shown. Figure 10 In an exemplary embodiment of the present invention, a B-level AMT measurement point curve diagram of the original data curve is a phase curve and reference data is provided in the XY mode; Figure 11 In an exemplary embodiment of the present invention, a B-level AMT measurement point curve diagram of the original data curve is provided, which is a phase curve and a reference data is constructed under the YX mode. Figure 12In an exemplary embodiment of the present invention, under the XY mode, the original data curve is the apparent resistivity curve, and a C-level AMT measurement point curve diagram for constructing reference data is provided. Figure 13 In the YX mode provided as an exemplary embodiment of the present invention, the original data curve is the apparent resistivity curve, and the C-level AMT measurement point curve of the constructed reference data is plotted. Figure 14 In an exemplary embodiment of the present invention, a C-level AMT measurement point curve diagram of the original data curve as the phase curve and the constructed reference data is provided in the XY mode. Figure 15 In an exemplary embodiment of the present invention, a C-level AMT measurement point curve diagram of the original data curve is a phase curve and reference data is provided under the YX mode; Figure 16 In an exemplary embodiment of the present invention, the original data curve is the apparent resistivity curve, and the D-level AMT measurement point curve of the constructed reference data is provided in the XY mode. Figure 17 In the YX mode provided as an exemplary embodiment of the present invention, the original data curve is the apparent resistivity curve, and the D-level AMT measurement point curve of the constructed reference data is plotted. Figure 18 In an exemplary embodiment of the present invention, a D-level AMT measurement point curve diagram of the original data curve as the phase curve and the constructed reference data is provided in the XY mode. Figure 19 The original data curve is a phase curve, and the D-level AMT measurement point curve diagram of the reference data is provided in the YX mode as an exemplary embodiment of the present invention. Detailed Implementation

[0025] The technical solution of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention. Furthermore, the technical features involved in the different embodiments of the present invention described below can be combined with each other as long as they do not conflict with each other.

[0026] See Figure 1 , Figure 1 A flowchart of a rapid evaluation method for audio magnetotelluric data quality provided by an exemplary embodiment of the present invention is shown, including the following steps: S1: Read the raw observation data of the audio magnetotelluric measurement point stored in the hard disk or memory. The raw observation data includes a frequency array and a response parameter array. The response parameter array includes an apparent resistivity array or an impedance phase array. S2: Preprocess the frequency array and the response parameter array by converting both into column vectors and arranging them in ascending order of frequency to obtain a sorted frequency sequence. and the sorted response parameter sequence ; S3: Perform a logarithmic transformation on the sorted frequency array (preferably a base-10 logarithmic transformation) to generate a logarithmic frequency sequence, which serves as the sorted and data-domain-transformed frequency array. ,Right now , The array index is used; and the data domain transformation of the response parameter sequence is performed according to the type of the response parameter array: when the response parameter array is an apparent resistivity array, a logarithmic transformation (preferably a base-10 logarithmic transformation) is performed on the sorted response parameter sequence to generate a logarithmic response sequence as the sorted and data domain transformed response parameter array. ,Right now When the response parameter array is an impedance phase array, the units of the impedance phases in the sorted response parameter sequence are converted to radians, thus obtaining the sorted and data domain-transformed response parameter array. ,Right now , Indicates radian conversion; S4: The first weighting is a cubic spline fitting based on the geometric weights of the second derivative: In this exemplary embodiment, instead of using traditional median filtering or equal-weighted spline smoothing, the weights for the first fitting are constructed using the geometric characteristics of the curve itself—the second derivative (curvature). Specifically: S41: Calculate the frequency array after sorting and data domain transformation. With the sorted and domain-transformed response parameter array Approximate value of the second derivative of For data points with non-uniform horizontal axis (logarithmic frequency) The approximate values ​​of the second derivative at each point are calculated using the non-uniform second-order difference formula. For the midpoint : ; For endpoints Forward or backward difference approximation can be used.

[0027] S42: Approximate the second derivative value The first weight is calculated by substituting the first linear truncation weight function based on a threshold: In this exemplary embodiment, the larger the absolute value of the second derivative, the greater the curvature at that point, and the more drastic the local change in the curve. In AMT data, points with large second derivatives often represent sharp changes caused by noise interference or static effects, while curve transitions reflecting the true electrical interface are usually relatively gentle. Therefore, in this exemplary embodiment, it is set that: points with larger second derivatives have smaller weights (closer to 0); points with smaller second derivatives have larger weights (closer to 1).

[0028] Therefore, in this step, a first linear truncation weighting function based on a threshold is used, and the second derivative approximation value in the first linear truncation weighting function is... The larger the point, the smaller the initial weight; the second derivative approximation. The smaller the point, the greater its initial weight, and the threshold in the first linear cutoff weight function is calculated using an approximation of the second derivative. The standard deviation multiplied by a preset factor; S43: Substitute the first weights into the weighted cubic smooth spline function to perform the first weighted spline fitting, and obtain the first fitted curve and the first fitted value of each point; S5: The second weighting is based on the deviation weights of the fitting residuals: In this exemplary embodiment, the first fitting in step S4 has initially filtered out obvious spikes and suppressed dead-band interference, but there may still be large deviations in some frequency bands due to strong interference. Therefore, in this exemplary embodiment, the first fitted values ​​are used in conjunction with the sorted and data-domain-transformed response parameter array. residual Construct a second set of weights: points with smaller residuals indicate that the original data is more consistent with the main trend and should be given higher weights; points with larger residuals may be strong interferences or local distortions that have not been completely filtered out and should be given lower weights.

[0029] Therefore, in this step, the first fitted values ​​are used in conjunction with the sorted and domain-transformed response parameter array. The residuals are used to calculate the second weight, employing a threshold-based second linear truncation weight function. In this function, points with smaller residuals have a larger second weight, and points with larger residuals have a smaller second weight. The threshold in the second linear truncation weight function is calculated based on the residuals. The average or residual the median; S6: The third weighting method is a forced weighting of dead frequency bands based on prior knowledge: Audio magnetotelluric data contains dead frequency bands with extremely low natural field energy in specific frequency ranges (typically 1000–5000 Hz, corresponding to logarithmic frequencies of approximately 3.0–3.7). Data quality in this band is generally poor, with apparent resistivity and phase curves exhibiting discontinuities, jumps, or even omissions. Conventional smoothing methods, without special processing for this region, will skew the fitted curve by these jump points, leading to distortion of the reference curve. This exemplary embodiment proposes a forced weighting strategy based on prior knowledge: In this step, the frequency index corresponding to the dead frequency band of the audio magnetotelluric data is identified based on prior knowledge. The frequency points within the dead frequency band are forcibly assigned a very small weight, while the other frequency points retain the second weight unchanged. The whole is used as the third weight. S7: Use the third weight to perform a second weighted spline fitting to obtain the final reference curve; S8: Read the data to be evaluated stored in the hard disk or memory, compare the data to be evaluated with the reference curve point by point, calculate the comprehensive similarity of each effective frequency point, determine the qualified frequency points based on the comprehensive similarity, count the proportion of qualified frequency points, calculate the comprehensive score, divide the quality level based on the comprehensive score and output the rating result.

[0030] After obtaining the reference curve, the data to be evaluated is quality-rated. The similarity calculation and weighted scoring in this step also reflect the adaptability of this invention to AMT band characteristics.

[0031] In summary, in this exemplary embodiment, the first weighting (curvature weight) utilizes the geometric characteristics of the electromagnetic response curve itself—the degree of local curvature reflected by the second derivative—to identify and suppress noise spikes. In audio magnetotelluric data, the disturbed frequency points often manifest as severe local jitter in the curve, with a curvature significantly greater than the gentle transitions reflecting the true electrical interface. Therefore, points with small curvature are given large weights, and points with large curvature are given small weights, enabling spline fitting to preserve the macroscopic electrical trend while suppressing local interference.

[0032] The second weighting (residual weighting) further identifies and eliminates outliers using the fitted residuals. After the first fitting, some frequency points that are strongly disturbed may still deviate from the main trend, and their residuals are significantly larger than those of normal frequency points. These are excluded from the second fitting by the residual weighting function, making the reference curve more robust.

[0033] The third weighting (forced deweighting in dead bands) is a specific design in this exemplary embodiment to address interference unique to audio-frequency magnetotelluric data. Natural audio-frequency magnetotelluric fields have extremely low energy in the 1000–5000 Hz frequency band, and measured data in this band generally exhibits severe distortion. Without processing, spline curves will be distorted by these jump points. Based on this prior knowledge, this exemplary embodiment forcibly sets the weight of this frequency band to a minimum value, allowing the spline curve to smoothly pass through this region without being pulled by interference points.

[0034] The reference curve obtained after triple weighted fitting is then compared with the data to be evaluated on a frequency-by-frequency basis. Through the fusion evaluation of apparent resistivity similarity and phase similarity, combined with the weighted ratio of qualified frequency points and frequency distribution density, the final output reflects the overall quality of the data.

[0035] Therefore, the core innovation of this exemplary embodiment lies in proposing a weighted spline fitting framework with three weights: curvature weight, residual weight, and dead-band forced weighting. Each weight has a clear physical meaning and works synergistically: First weighting (curvature weight): Utilizing the geometric features of the curve, it automatically suppresses local high-frequency spikes and static effect peaks, preserving smooth macroscopic electrical interface information; Second weighting (residual weight): Based on the residuals after the previous fitting, it further eliminates outlier interference points, making the reference curve more faithful to the main electrical trend and avoiding overfitting; Third weighting (dead-band forced weighting): For the dead-band interference unique to audio-to-geomagnetic signals, based on prior knowledge, the weights of the corresponding frequency points are reset to a minimum value (e.g., 0.001), allowing the fitted curve to smoothly pass through this region without being skewed by jump points.

[0036] Furthermore, it requires no assumptions about underground models and is data-driven: Traditional full tensor fitting methods rely on the assumption of a one-dimensional anisotropic model, and the fitting results are unreliable when the underground structure is complex. This exemplary embodiment only utilizes the curvature characteristics, prior frequency band knowledge, and residual statistical characteristics of the original data itself, without introducing any geological model priors. The constructed reference curve objectively reflects the main trend of the data and is applicable to any electrical structure (one-dimensional, two-dimensional, three-dimensional, including static effects, etc.).

[0037] More preferably, in an exemplary embodiment, in step S4, the first linear truncation weight function based on the threshold is: When the curvature threshold is zero, the first weight is one; when the absolute value of curvature is between zero and the curvature threshold, the first weight decreases linearly; when the absolute value of curvature is greater than or equal to the curvature threshold, the first weight is truncated to zero. The curvature threshold is adaptively determined based on the overall curvature distribution of the data, that is, the curvature threshold is equal to a preset multiple multiplied by the standard deviation of the absolute value of curvature.

[0038] Specifically, in this exemplary embodiment, the technical principle of adaptively determining the curvature threshold is as follows: The apparent resistivity amplitude and curve shape variation ranges vary greatly at different measuring points. If a fixed threshold is used, it may be too strict for some measuring points (misjudging normal electrical interfaces as noise), and too lenient for others (failing to effectively suppress noise). This exemplary embodiment uses the standard deviation of the absolute curvature multiplied by a preset multiple as the threshold, ensuring that the threshold matches the curvature distribution scale of the data itself, eliminating the need for individual parameter tuning for each measuring point. The linear truncation function takes the form: when the absolute curvature is zero, the weight is 1 (completely retained); when the absolute curvature reaches the threshold, the weight is 0 (completely excluded), decreasing linearly in between.

[0039] in That is, the first weight, This is the curvature threshold parameter. The characteristic of this function is: when... hour ;when The time weights decrease linearly; when The time weight is truncated to 0. It can be adaptively determined based on the overall curvature distribution of the data: ,in The preset multiple (preferred value is 3 to 5) is used. Indicates standard deviation, This is an approximation of the second derivative. This adaptive design ensures that the weights are insensitive to the scale of the data itself.

[0040] More preferably, in an exemplary embodiment, in step S43, the first weight is... Substituting the weighted cubic smoothing spline function into the first weighted spline fit yields the first fitted curve and the first fitted values ​​for each point, specifically: Set smoothing parameters (The preferred value is 0.985, which is close to interpolation but retains a certain degree of smoothness), solve the weighted least squares problem: ; Using the variational method, the solution to this problem can be obtained as a cubic smooth spline function. Let this solution be denoted as the first fitted curve. Calculate the first fitted value for each point. In the formula, The function f represents the fitted curve. Regarding the independent variable The second derivative, where t represents the integration variable and the integration range is... .

[0041] It should be noted that the above minimization problem is solved directly using the csaps function (cubic smooth spline) in the MATLAB Curve FittingToolbox.

[0042] The core of the csaps function is to find a cubic spline. In the sense of minimizing the above objective function, the data is smoothly fitted. This function introduces a smoothing parameter. (Values ​​range from 0 to 1), seeking an optimal balance between "fitting accuracy" (the first term) and "curve smoothness" (the second term, i.e., the integral penalty for the square of the second derivative). When the curve degenerates into an interpolation spline passing through all data points, when... When the curve approaches the least squares straight line fit, the curve becomes close to the least squares straight line fit.

[0043] More preferably, in an exemplary embodiment, in step S5, the second linear truncation weight function based on the threshold is: The second weight is one when the residual is zero; the second weight decreases linearly when the residual is between zero and the residual mean; and the second weight is truncated to zero when the residual is greater than or equal to the residual mean. The residual mean is the arithmetic mean of the residuals.

[0044] Specifically, in this exemplary embodiment, after the first weighted fitting, the fitting residuals of most normal frequency points are small and concentrated near the mean, while the residuals of strongly disturbed outliers are significantly larger than the mean. Using the mean as the boundary, frequency points with residuals less than the mean participate in the second fitting, while frequency points with residuals greater than or equal to the mean are excluded (weight 0). This "hard truncation" strategy, compared to common "soft decay" strategies (such as exponential decay), can more thoroughly eliminate outliers. Using the mean instead of the median as the threshold, although the mean is more sensitive to outliers, when the proportion of outliers is low (usually below 20%), the mean can better reflect the average level of normal residuals, and the calculation is simpler. The linear truncation function ensures a smooth transition of weights from 1 to 0, avoiding fitting instability that might result from a step change. Specifically: Step S51: Calculate the residual and determine the threshold Calculate the absolute residual at each point In the formula, This represents the sorted and data-domain-transformed response parameter array y. i The value, Indicates the first fitted curve exist The fitted value at that point. Calculate the mean of the residuals. (The median can also be used instead. The mean is more sensitive to outliers, while the median is more robust. The mean is preferred.) n is the number of points.

[0045] Step S52: Weight assignment based on residuals – Residual weighting function. This exemplary embodiment proposes a second threshold-based linear truncation weighting function:

[0046] in, That is, the second weight. The characteristic of this function is: when... hour, ;when The time weights decrease linearly; when The time weight is truncated to 0. This means that points whose residuals exceed the average residual will not be included in the second fitting (weight is 0), thus completely eliminating outlier points. For a more gentle decay curve, an exponential curve can also be used. However, the linear truncation function is simpler to calculate and has a clear physical meaning.

[0047] It should be noted that both weighting functions (curvature weighting function and residual weighting function) proposed in the above exemplary embodiments adopt... The unified form of the function has the following advantages: (1) Clear physical meaning: the weight is 1 when value=0, and the weight is 0 when value≥T, with a linear transition in the middle; (2) Interpretable threshold parameter: T directly represents "the maximum allowable deviation / curvature", and it is completely excluded if it exceeds the limit; (3) Extremely simple calculation: only one division and two comparisons are required, which is suitable for large-scale batch processing; (4) Easy to adapt: ​​T can be automatically set according to data statistics (standard deviation, mean, etc.).

[0048] More preferably, in an exemplary embodiment, in step S6, the step of identifying the frequency index corresponding to the dead frequency band of the audio magnetotelluric data based on prior knowledge specifically involves: S61: Determine the frequency range of the dead band based on the acquired parameters. The dead frequency band has a frequency range of 1000 Hz to 5000 Hz; a set of frequency indexes falling within this frequency range is determined on the logarithmic frequency axis. For example, when the dead band corresponds to the 46th to 54th frequency points, .

[0049] S62: For Set its weight to a preset minimum value. (The preferred value is 0.001), that is: ; For other frequency points, the third weighting The second weight value is maintained. This minimal weight ensures that data points within the dead band hardly participate in the process of minimizing the spline fitting error, allowing the fitted curve to smoothly "pass through" this region without being pulled by jump points.

[0050] Specifically, in this exemplary embodiment, the dead frequency band of the audio-frequency magnetotelluric method, ranging from 1000 to 5000 Hz, originates from the energy trough of the natural electromagnetic field in this frequency band. This physical phenomenon is objectively real and does not change with the location of the measurement point or geological conditions; therefore, it can be solidified as prior knowledge in the method. Within this frequency band, the measured apparent resistivity and phase curves typically exhibit large jumps or even missing data. If normal weights are used for fitting, the spline curve will be pulled in the wrong direction by these outliers. The consideration for setting the weights to 0.001 (instead of 0) is twofold: firstly, to ensure that the data in this frequency band hardly participates in the fitting (its contribution to the objective function is negligible); and secondly, to maintain the continuity of the spline curve in this region, avoiding instability in the spline interpolation values ​​caused by completely zero weights.

[0051] More preferably, in an exemplary embodiment, in step S7, the weights determined in the third step are utilized. The final reference curve constructed in this invention is fitted.

[0052] Step S71: Second weighted spline fitting Use the updated weights Then call the weighted cubic smoothed spline function again, with the smoothing parameter being... The third fitted curve was obtained by solving the problem. Calculate the fitted value at each point. This curve is the final reference curve constructed in this invention. The calculation method is similar to that in step S3: Set smoothing parameters (The preferred value is 0.985, which is close to interpolation but retains a certain degree of smoothness), solve the weighted least squares problem: ; Using the variational method, the solution to this problem can be obtained as a cubic smooth spline function. Let this solution be denoted as the first fitted curve. Calculate the first fitted value for each point. In the formula, The function f represents the fitted curve. Regarding the independent variable The second derivative, where t represents the integration variable and the integration range is... .

[0053] It should be noted that the above minimization problem is solved directly using the csaps function (cubic smooth spline) in the MATLAB Curve FittingToolbox.

[0054] The core of the csaps function is to find a cubic spline. In the sense of minimizing the above objective function, the data is smoothly fitted. This function introduces a smoothing parameter. (Values ​​range from 0 to 1), seeking an optimal balance between "fitting accuracy" (the first term) and "curve smoothness" (the second term, i.e., the integral penalty for the square of the second derivative). When the curve degenerates into an interpolation spline passing through all data points, when... When the curve approaches the least squares straight line fit, the curve becomes close to the least squares straight line fit.

[0055] Step S72: Inverse Transformation and Reference Curve Output: Perform an inverse transformation based on the response parameter type: if it is apparent resistivity, then If it is an impedance phase, then .

[0056] Restore the original order according to the reverse mapping sorted in step S2, and output the optimized reference response parameter array (including reference apparent resistivity). , and reference phase , This reference curve serves as the benchmark for subsequent quality evaluation.

[0057] More preferably, in an exemplary embodiment, in step S8, after obtaining the reference curve, the data to be evaluated is given a quality rating.

[0058] Step S81: Obtain the data to be evaluated and the reference data: Obtain the response parameter data of the measurement point to be evaluated, including the frequency array. Apparent resistivity array , and impedance phase array , Simultaneously, obtain the reference data output in step S7. .

[0059] S82: Since the frequency sampling of the data to be evaluated and the reference data may be inconsistent, frequency point alignment is required. The step of comparing the data to be evaluated with the reference curve frequency by frequency point includes the frequency point alignment step: Taking the base-10 logarithm of the frequency arrays of the data to be evaluated and the reference data respectively, we obtain... , ; For each frequency point in the data to be evaluated, find the nearest frequency point in the logarithmic frequency sequence of the reference data. ; If the minimum distance is less than the preset tolerance δ (preferred value is 0.02 logarithmic units, corresponding to approximately 4.7% relative frequency error), the matching is accepted; otherwise, the frequency point to be evaluated is discarded, resulting in a set of effective frequency points and the corresponding subset of reference data.

[0060] Specifically, in this exemplary embodiment, in actual audio magnetotelluric data acquisition, the frequency point lists obtained from different measurement points, different acquisition devices, or different processing procedures may not be completely consistent. For example, some data may contain 60 frequency points, while others may only contain 55, and the frequency values ​​may have slight differences. Directly comparing each point will lead to frequency misalignment. The alignment principle of this embodiment is: using the logarithmic frequency as the matching space, for each frequency point of the data to be evaluated, the nearest frequency point is searched in the frequency sequence of the reference data. Only when the distance is less than a preset tolerance is it considered a valid match. A tolerance of 0.02 logarithmic units corresponds to a relative frequency error of approximately 4.7%. This tolerance range is sufficient to cover frequency deviations that occur during normal acquisition and processing, while also preventing incorrect matching of response parameters at different frequency points.

[0061] More preferably, in an exemplary embodiment, step S8 further includes: Step S83: Select the evaluation mode. Based on the user-specified evaluation mode, extract the corresponding response components from the data to be evaluated and the reference data. If the mode is Extract only apparent resistivity in direction and phase If the pattern is Extract only apparent resistivity in direction and phase If the pattern is Then handle them separately. and The final score is a weighted average of the two factors, and the direction is determined by the final score.

[0062] More preferably, in an exemplary embodiment, in step S8, step S84: the comprehensive similarity is calculated in the following manner: Calculate the apparent resistivity similarity separately Phase similarity Then, the apparent resistivity similarity and phase similarity are weighted and fused into a comprehensive similarity. ; The apparent resistivity similarity is obtained by linearly mapping the logarithmic domain relative error to the 0-1 interval. Let... , This represents the absolute value of the difference between the logarithm of the data to be evaluated and the reference data, each taken with respect to 10. This represents the apparent resistivity data to be evaluated. This indicates the apparent resistivity data for reference. Then, the apparent resistivity similarity : ; in The apparent resistivity tolerance parameter is preferably 0.3 (i.e., a logarithmic difference of 0.3 corresponds to approximately 2 times the relative error, at which point the similarity drops to 0).

[0063] The phase similarity is obtained by linearly mapping the absolute value of the phase difference to the 0-1 interval after considering a 360-degree periodic rotation: ; in, Indicates the impedance phase to be evaluated. Indicates the reference impedance phase, The phase similarity is expressed as the absolute difference between the phase of the impedance to be evaluated and the phase of the reference impedance, limited to between -180° and 180°. : ; in The phase tolerance parameter is preferably 30° (i.e., the similarity drops to 0 when the phase difference is 30°).

[0064] Fusion: Equal-weight fusion is adopted. .

[0065] Step S85: Qualified frequency point determination Set similarity threshold (The preferred value is 0.8). When If the frequency point is within acceptable limits, mark it as an acceptable frequency point; otherwise, mark it as an unacceptable frequency point. Count the number of acceptable frequency points. and the total number of effective frequency points Calculate the pass rate: R .

[0066] Specifically, in the above exemplary embodiment, the physical meaning of the relative error calculated by apparent resistivity similarity in the logarithmic domain is as follows: In electromagnetic exploration, apparent resistivity often spans several orders of magnitude (from a few Ω·m to several thousand Ω·m). In the logarithmic domain, a deviation of one logarithmic unit corresponds to a relative error of 10 times, which can more reasonably reflect the degree of relative deviation at different amplitude levels. The tolerance parameter 0.3 means that when the apparent resistivity differs by approximately 2 times (10^0.3≈2.0), the similarity drops to 0, meaning that the frequency point is considered unqualified. The phase similarity considers a 360° periodicity because the impedance phase has a periodicity of ±180° when measured in degrees. For example, the actual difference between -170° and +190° is only 20°, not 360°. The phase tolerance of 30° is derived from industry experience—in high-quality audio-frequency magnetotelluric data, the phase curve is usually relatively smooth, and the phase change between adjacent frequency points is generally within 10° to 20°. A deviation of 30° is already significantly outside the normal range. The reason for the equal weighting of apparent resistivity and phase (each accounting for 50%) is that both independently reflect the electrical information of the underground medium and are equally important.

[0067] More preferably, in an exemplary embodiment, the calculation of the comprehensive score adopts a weighted scoring method: Weights are assigned based on the distribution density of valid frequency points on the logarithmic frequency axis. The weight of each frequency point is determined according to the width of the frequency range represented by that frequency point on the logarithmic frequency axis. The endpoints are spaced on one side, and the intermediate frequency points are spaced by half the distance between adjacent frequency points. After normalizing the weights of each frequency point, the normalized weights are multiplied by the qualified mark of that frequency point and summed to obtain the weighted comprehensive score.

[0068] Specifically, in this exemplary embodiment, the user selects to use a weighted scoring method: If weighted Then, weights are assigned based on the distribution density of frequency points on the logarithmic frequency axis. Specifically, for the logarithmic frequency of effective frequency points... After sorting, calculate the width of the frequency range represented by each frequency point, which is the corresponding weight: ; intermediate frequency weight Normalized weights The weighted composite score is then:

[0069] In yet another exemplary embodiment, if unweighted Then the overall score .

[0070] It should be noted that traditional equal-weighted scoring methods treat the contribution of each effective frequency point equally, which is reasonable when the frequency points are evenly distributed. However, audio magnetotelluric data typically has frequency points arranged at logarithmic intervals along the frequency axis, with sparse frequency points in the low-frequency band (representing a wider frequency range) and dense frequency points in the high-frequency band (representing a narrower frequency range). If equal-weighted scoring is used, the dense frequency points in the high-frequency band will dominate the overall score, diluting the contribution of the low-frequency band, resulting in a score that does not truly reflect the usability of the data across the entire frequency band. In the frequency point weight calculation method used in this exemplary embodiment, the weight of the middle frequency point is equal to the sum of half the distance between its two adjacent frequency points (i.e., the width of the frequency range represented by that point), and the endpoints use unilateral spacing. After normalization, the sum of the weights of all frequency points is 1. This "area-weighted" strategy ensures that the contribution of each frequency point is proportional to the actual bandwidth it represents.

[0071] In other words, in this exemplary embodiment, after obtaining the reference curve, the data to be evaluated is compared with the reference curve frequency by frequency, and the comprehensive similarity of each frequency point is calculated (apparent resistivity similarity and phase similarity are equally weighted). Then, the effective frequency points are weighted according to their actual representative bandwidth on the logarithmic frequency axis, and the weighted similarity ratio is calculated as the comprehensive score. Specifically, the weight of each frequency point... (Endpoints use unilateral spacing), calculated after normalization. This scoring strategy effectively overcomes the distortion of evaluation results caused by an excessive density of frequency points within the dead frequency band.

[0072] More preferably, in an exemplary embodiment, the quality grade includes four levels: Excellent, Good, Intermediate, and Unqualified. The classification rule is as follows: an overall score greater than or equal to 0.95 is rated as Excellent; an overall score less than 0.95 but greater than or equal to 0.80 is rated as Good; an overall score less than 0.80 but greater than or equal to 0.60 is rated as Intermediate; and an overall score less than 0.60 is rated as Unqualified. Numerically, Grade A (Excellent) is: Grade B (Good): Grade C (Medium): Grade D (Unqualified): .

[0073] It should be noted that the entire method, from data input to reference curve construction and quality rating, requires only a few default parameters (smoothing parameter 0.995, dead band index, etc.), and all thresholds can be automatically calculated. Actual testing shows that the total processing time for a single measurement point (approximately 60 frequency points) on a standard PC is less than 1.2 seconds; batch evaluation of 100 measurement points takes approximately 2 minutes. Compared to manual visual evaluation (4 people working for 5 days) and traditional iterative inversion methods (tens of seconds to several minutes per point), the efficiency is improved by hundreds to tens of thousands of times.

[0074] To verify the effectiveness and reliability of this exemplary embodiment, synthetic data testing and practical applications were conducted.

[0075] (1) Synthetic data testing To verify the effectiveness of the method of this invention in constructing audio-frequency magnetotelluric reference curves, a one-dimensional four-layer theoretical model was designed: resistivity is respectively , , , The layer thicknesses are respectively , , (The bottom half-space is infinitely thick). The theoretical apparent resistivity and phase curves were obtained using one-dimensional AMT forward modeling, and 10% Gaussian random noise was added to the entire frequency band. To simulate the "dead band" interference in actual acquisition, jump noise with significantly discontinuous amplitude was introduced in the 1000–5000 Hz frequency band, causing severe distortion of the curves in this range.

[0076] Figure 2 and Figure 3 The original noisy curve is shown. Figure 2 For diamond-shaped red dots, Figure 3 The comparison results of the square red dots, theoretical curves (light blue lines), and reference curves constructed using different methods are all based on the XY model. Among them: (1), where Figure 2 This corresponds to the case where the original curve of the noisy synthetic data is an apparent resistivity curve. Figure 3 The corresponding case is when the original curve of the noisy synthetic data is an impedance phase curve; (2) The methods compared include: unweighted cubic spline interpolation fitting (black line, i.e. Figure 2 and Figure 3 Unweighted spline fitting), five-point mean filtering (dark blue line, i.e.) Figure 2 and Figure 3 The five-point mean filter and the five-point median filter (light green line) are shown in the image. Figure 2 and Figure 3 The 5-point median filtering in the model and the cubic weighted cubic spline fitting method proposed in this invention (pink line, i.e.) Figure 2 and Figure 3 Reference data from the present invention patent). The horizontal axis represents frequency, and the vertical axis represents apparent resistivity (…). Figure 2 Phase () Figure 3 ).

[0077] The results show that all four methods can effectively suppress local glitches caused by single-frequency spikes, and the constructed reference curves match the theoretical curves well. However, in the 1000–5000 Hz frequency band where "dead bands" exist, the unweighted spline fitting curve is significantly skewed by jump points. Although 5-point mean filtering and 5-point median filtering can smooth out some jumps, the overall trend deviates significantly from the theoretical curve and cannot reproduce the true electrical response. In contrast, the method of this invention suppresses glitches through curvature weighting, forces weight reduction in dead bands to avoid being pulled by outliers, and further eliminates outliers through residual weighting. The constructed reference curve is highly consistent with the theoretical curve throughout the entire frequency band (especially in the dead band region).

[0078] In summary, the synthetic data test results show that when facing continuous strong interference in dead bands, which is common in actual AMT data, the method of this invention has significantly better anti-interference ability and curve reconstruction accuracy compared with traditional spline fitting and conventional smoothing filtering, and can provide a more reliable reference benchmark for subsequent data quality evaluation.

[0079] (2) Practical application effect To verify the effectiveness and engineering applicability of this invention in actual audio magnetotelluric (AMT) data, a batch quality assessment was conducted on 96 AMT measurement data points collected in a potash exploration area in a plain. The testing environment was MATLAB R2017b. The program automatically completed the following tasks: plotting the logarithmic coordinates of the original curves and reference curves, calculating the similarity of the XY / YX modes at each measurement point and assigning weighted scores, classifying the quality levels (A / B / C / D), and statistically summarizing the rating results for all measurement points. Finally, the evaluation results (including curves, scores, and levels) for all 96 measurement points were automatically output as a PDF document. The entire processing took 99 seconds, translating to an average processing time of approximately 1 second per measurement point (including plotting and PDF writing).

[0080] Some rating results are as follows Figures 4-19 The horizontal axis represents frequency, and the vertical axis represents apparent resistivity ( ). Figure 4 , Figure 5 , Figure 8 , Figure 9 , Figure 12 , Figure 13 , Figure 16 , Figure 17 Phase () Figure 6 , Figure 7 , Figure 10 , Figure 11 , Figure 14 , Figure 16 , Figure 18 , Figure 19 Specifically: Figures 4-7 The AMT measurement point curve for Class A is shown. Figure 4In XY mode, the original data curve is the apparent resistivity curve, and a Class A AMT measurement point curve of the reference data is constructed. Figure 5 In YX mode, the original data curve is the apparent resistivity curve, and a Class A AMT measurement point curve of the reference data is constructed. Figure 6 In XY mode, with the original data curve as the phase curve, construct the Class A AMT measurement point curve of the reference data; Figure 7 In YX mode, using the original data curve as the phase curve, a Grade A AMT measurement point curve of the reference data was constructed. Evaluation results: Score=0.95, Grade=A.

[0081] Figures 8-11 The AMT measurement point curve for Class B is shown. Figure 8 In XY mode, the original data curve is the apparent resistivity curve, and a B-level AMT measurement point curve of the reference data is constructed. Figure 9 In YX mode, the original data curve is the apparent resistivity curve, and a B-level AMT measurement point curve of the reference data is constructed. Figure 10 In XY mode, with the original data curve as the phase curve, construct a B-level AMT measurement point curve diagram of the reference data; Figure 11 In YX mode, using the original data curve as the phase curve, a B-grade AMT measurement point curve of the reference data was constructed. Evaluation results: Score=0.82, Grade=B.

[0082] Figures 12-15 The AMT measurement point curve for Class C is shown. Figure 12 In XY mode, the original data curve is the apparent resistivity curve, and a C-level AMT measurement point curve of the reference data is constructed. Figure 13 In YX mode, the original data curve is the apparent resistivity curve, and a C-level AMT measurement point curve of the reference data is constructed. Figure 14 In XY mode, the original data curve is the phase curve, and a C-level AMT measurement point curve of the reference data is constructed. Figure 15 In YX mode, using the original data curve as the phase curve, a C-level AMT measurement point curve of the reference data was constructed. Evaluation results: Score=0.72, Grade=C.

[0083] Figures 16-19 The AMT measurement point curve for Class D is shown. Figure 16 In XY mode, the original data curve is the apparent resistivity curve, and a D-level AMT measurement point curve of the reference data is constructed. Figure 17 In YX mode, the original data curve is the apparent resistivity curve, and a D-level AMT measurement point curve of the reference data is constructed. Figure 18 In XY mode, the original data curve is the phase curve, and a D-level AMT measurement point curve of the reference data is constructed. Figure 19In YX mode, using the original data curve as the phase curve, a D-level AMT measurement point curve of the reference data was constructed. Evaluation results: Score=0.30, Grade=D.

[0084] The rating results of the method of this invention were compared with the visual rating results of two senior experts. The experts independently judged the apparent resistivity-phase curves according to industry standards, and gave a final rating after reaching a consensus, taking a total of one day. The comparison results showed that out of 96 measurement points, the rating of this invention was completely consistent with the expert ratings on 93 points, achieving a consistency rate of 96.9%. Only 3 measurement points showed discrepancies—the experts rated them as Grade B (Good), while this invention rated them as Grade C (Medium). Subsequent review revealed that these 3 measurement points were all located in areas with strong interference (near villages, with superposition of power frequency harmonics and dead-band interference), and the phase curves showed multiple discontinuous jumps. The experts' judgments contained some subjective differences, while the method of this invention strictly followed the weighted scoring threshold, resulting in a more objective result.

[0085] It is worth noting that the expert team's one-day evaluation time does not include the results compilation and archiving, while the method of this invention completes all processing, rating, and automatic generation of PDF reports in less than two minutes, improving efficiency by more than a thousand times. Furthermore, the evaluation results have advantages such as traceability and strong consistency. This practical application demonstrates that the method of this invention can significantly improve the efficiency and objectivity of quality evaluation of massive audio-visual magnetotelluric data, and possesses good engineering application value.

[0086] Obviously, the above embodiments are merely illustrative examples for clear explanation and are not intended to limit the implementation. Those skilled in the art can make other variations or modifications based on the above description. It is neither necessary nor possible to exhaustively list all possible implementations. However, obvious variations or modifications derived therefrom are still within the scope of protection of this invention.

Claims

1. A rapid evaluation method for the quality of audio magnetotelluric data, characterized in that, Includes the following steps: Read the raw observation data of the audio magnetotelluric measurement points stored on the hard disk or in memory. The raw observation data includes a frequency array and a response parameter array. The response parameter array includes an apparent resistivity array or an impedance phase array. The frequency array and the response parameter array are preprocessed to convert them into column vectors and sorted in ascending order of frequency to obtain sorted frequency sequences and sorted response parameter sequences. Perform a logarithmic transformation on the sorted frequency array to generate a logarithmic frequency sequence, which serves as the sorted frequency array after data domain transformation. ; and perform data domain transformation on the sorted response parameter sequence according to the type of the response parameter array: when the response parameter array is an apparent resistivity array, perform a logarithmic transformation on the sorted response parameter sequence to generate a logarithmic response sequence as the sorted and data domain transformed response parameter array. ; When the response parameter array is an impedance phase array, the units of the impedance phases in the sorted response parameter sequence are converted to radians, thus obtaining the sorted and data domain-transformed response parameter array. , The array index; The first weighting is a cubic spline fit based on geometric weights of the second derivative: calculating the frequency array after sorting and data domain transformation. With the sorted and domain-transformed response parameter array Approximate value of the second derivative Approximating the second derivative The first weight is calculated by substituting the first linear truncation weight function based on the threshold, and the approximate value of the second derivative in the first linear truncation weight function is obtained. The larger the point, the smaller the initial weight; the second derivative approximation. The smaller the point, the greater its initial weight. The threshold in the first linear cutoff weight function is calculated using an approximation of the second derivative. Multiply the standard deviation by a preset factor; substitute the first weight into the weighted cubic smooth spline function to perform the first weighted spline fitting, and obtain the first fitted curve and the first fitted value of each point; The second weighting is a bias weight based on the fitting residuals: the first fitted values ​​are compared with the sorted and domain-transformed response parameter array. residual The second weight is calculated using a threshold-based second linear truncation weight function. In this function, points with smaller residuals have larger second weights, and points with larger residuals have smaller second weights. The threshold in the second linear truncation weight function is calculated using the residuals... The average or residual the median; The third weighting is a forced weight reduction of dead frequency bands based on prior knowledge: the frequency index corresponding to the dead frequency band of audio magnetotelluric data is identified based on prior knowledge, and the frequency points in the dead frequency band are forcibly set with minimal weights, while the weights of the other frequency points remain unchanged from the second weight. The overall weight is used as the third weight. A second weighted spline fitting is performed using the third weight to obtain the final reference curve; Read the data to be evaluated stored in the hard disk or memory, compare the data to be evaluated with the reference curve point by point, calculate the comprehensive similarity of each effective frequency point, determine the qualified frequency points based on the comprehensive similarity, count the proportion of qualified frequency points, calculate the comprehensive score, divide the quality level based on the comprehensive score and output the rating result.

2. The method according to claim 1, characterized in that, The threshold-based first linear truncation weight function is: When the curvature threshold is zero, the first weight is one; when the absolute value of curvature is between zero and the curvature threshold, the first weight decreases linearly; when the absolute value of curvature is greater than or equal to the curvature threshold, the first weight is truncated to zero. The curvature threshold is adaptively determined based on the overall curvature distribution of the data, that is, the curvature threshold is equal to a preset multiple multiplied by the standard deviation of the absolute value of curvature.

3. The method according to claim 1, characterized in that, The threshold-based second linear truncation weight function is: The second weight is one when the residual is zero; the second weight decreases linearly when the residual is between zero and the residual mean; and the second weight is truncated to zero when the residual is greater than or equal to the residual mean. The residual mean is the arithmetic mean of the residuals.

4. The method according to claim 1, characterized in that, The frequency index corresponding to the dead frequency band of audio magnetotelluric data based on prior knowledge identification is specifically as follows: The frequency range of the dead band is determined based on the collected parameters, and the frequency range of the dead band is from 1000 Hz to 5000 Hz; the set of frequency point indices falling within this frequency range is determined on the logarithmic frequency axis. The minimum weight is 0.

001.

5. The method according to claim 1, characterized in that, The step of comparing the data to be evaluated with the reference curve frequency by frequency point includes a frequency point alignment step: Take the base-10 logarithm of the frequency arrays of the data to be evaluated and the reference data respectively; For each frequency point in the data to be evaluated, the nearest frequency point is found in the logarithmic frequency sequence of the reference data. If the minimum distance is less than the preset tolerance, the match is accepted; otherwise, the frequency point is discarded, resulting in a set of effective frequency points and the corresponding subset of reference data. The preset tolerance is 0.02 logarithmic units.

6. The method according to claim 1, characterized in that, The overall similarity is calculated in the following way: Apparent resistivity similarity and phase similarity are calculated separately, and then the apparent resistivity similarity and phase similarity are fused together with equal weight to form a comprehensive similarity. The apparent resistivity similarity is obtained by linearly mapping the relative error in the logarithmic domain to the interval of 0 to 1, and the phase similarity is obtained by linearly mapping the absolute value of the phase difference to the interval of 0 to 1 after considering the 360-degree periodic rotation. In the process of determining qualified frequency points based on comprehensive similarity, the qualified threshold for comprehensive similarity is 0.

8.

7. The method according to claim 1, characterized in that, The overall score is calculated using a weighted scoring method. Weights are assigned based on the distribution density of valid frequency points on the logarithmic frequency axis. The weight of each frequency point is determined according to the width of the frequency range represented by that frequency point on the logarithmic frequency axis. The endpoints are spaced on one side, and the intermediate frequency points are spaced by half the distance between adjacent frequency points. After normalizing the weights of each frequency point, the normalized weights are multiplied by the qualification mark of that frequency point and summed to obtain the weighted comprehensive score.

8. The method according to claim 1, characterized in that, The quality grades include four levels: excellent, good, medium, and unqualified. The classification rules are as follows: an overall score greater than or equal to 0.95 is rated as excellent; an overall score less than 0.95 but greater than or equal to 0.80 is rated as good; an overall score less than 0.80 but greater than or equal to 0.60 is rated as medium; and an overall score less than 0.60 is rated as unqualified.

9. The method according to claim 1, characterized in that, The smoothing parameter for the first weighted spline fitting is 0.985, and the smoothing parameter for the second weighted spline fitting is 0.955.

Images