CN114460654B - Semi-aviation transient electromagnetic data inversion method and device based on L1L2 mixed norm - Google Patents

Abstract

The invention discloses a semi-aviation transient electromagnetic data inversion method based on L1L2 mixed norm, which adopts ground to transmit transient electromagnetic signals and utilizes an unmanned aerial vehicle carrying and receiving device to obtain electromagnetic response data, and comprises the following steps: collecting electromagnetic response data of a region to be detected by using a receiving device; constructing a stratum model of a region to be detected, and presetting an iteration termination condition; obtaining a transverse constraint weighting matrix corresponding to an L2 norm constraint term of the electromagnetic response data; solving or updating a longitudinal constraint matrix corresponding to the L1 norm constraint term of the electromagnetic response data; forming a target function by the data fitting gradient term, the L1 norm constraint term and the L2 norm constraint term, and solving an expansion formula after the incremental derivative of the stratum model is equal to zero; and (4) calculating the increment of the stratum model by using an expansion formula through inversion calculation, correcting the stratum model, and obtaining an inversion imaging image after the iteration termination condition is reached.

Description

Method and device for inversion of semi-aviation transient electromagnetic data based on L1L2 hybrid norm

technical field

The invention relates to the technical field of geophysical aviation electromagnetic exploration, in particular to a semi-aerial transient electromagnetic data inversion method and device based on L1L2 mixed norm.

Background technique

At present, the geophysical electromagnetic exploration in the prior art mostly adopts the aerial electromagnetic method, which uses a helicopter or fixed-wing aircraft to carry a launch and observation system, which cannot be used in urban underground space exploration, ground geological survey and other scenarios. In addition, the airborne electromagnetic method in the prior art adopts a one-dimensional single-point inversion algorithm, which has no constraints in the lateral direction, and the continuity of the profile resistivity is poor when the data noise is strong; regularization is usually used in the longitudinal constraints on the resistivity. The main disadvantage of this regularization constraint is that it makes the model too smooth and cannot effectively describe the electrical mutation interface information.

For example, the patent application number is "202110820825.8" and the Chinese invention patent titled "A Semi-Aeronautical Transient Electromagnetic Conductivity-Depth Imaging Method and Equipment", which pre-simulates the electromagnetic response according to the actual situation to query the database, and then The detected electromagnetic response data is searched in the "library", and there is no need to perform multiple iterative calculations like inversion, so that the semi-aviation transient electromagnetic method can perform fast imaging, and quickly obtain preliminary imaging results and inversion initial models.

Another example is the Chinese invention patent with the patent application number of "202010751396.9" and the title of "A Method for Analysis and Interpretation of UAV Semi-Aeronautical Time Domain Electromagnetic Detection Data", which includes the following steps: a. Preprocessing the survey line data, eliminating Motion noise; b. Perform spectrum analysis and digital filtering on the secondary field data; c. Form data to be imaged or interpreted by inversion; d. Perform one-dimensional fast inversion of the secondary field data after stacking and extraction, Establish a subsurface electrical structural profile, and perform geological interpretation of the inversion imaging results in combination with geological data to form a comprehensive interpretation result.

The above techniques have the following problems: the quick imaging results of the above techniques are relatively rough, and cannot accurately reflect the formation resistivity value and layer interface information; for the noisy data, the single-point inversion does not consider the influence of the formation information of the adjacent points of the measurement point. , the inversion results are horizontally discontinuous and appear striped.

Therefore, it is urgent to propose a simple, accurate and reliable semi-aerial transient electromagnetic data inversion method based on L1L2 mixed norm.

SUMMARY OF THE INVENTION

In view of the above problems, the purpose of the present invention is to provide a semi-aviation transient electromagnetic data inversion method based on L1L2 mixed norm, and the technical scheme adopted in the present invention is as follows:

In the first part, the present invention provides a semi-aviation transient electromagnetic data inversion method based on L1L2 hybrid norm, which adopts ground-transmitting transient electromagnetic signals, and utilizes a UAV-mounted receiving device to obtain electromagnetic response data, which includes the following steps:

Use the receiving device to collect the electromagnetic response data of the area to be measured;

Build the stratigraphic model of the area to be detected, and preset the iteration termination conditions;

Obtain the transverse constraint weighting matrix corresponding to the L2 norm constraint term of the electromagnetic response data;

Obtain or update the longitudinal constraint matrix corresponding to the L1 norm constraint term of the electromagnetic response data;

The data fitting gradient term, the L1 norm constraint term, and the L2 norm constraint term constitute the objective function, and the expansion formula after the incremental derivation of the formation model is equal to zero is obtained;

The formation model increment is obtained by inversion calculation using the expansion formula, and the formation model is corrected, and the inversion imaging image is obtained after the iteration termination condition is reached.

The second part, the present invention provides a semi-aviation transient electromagnetic data inversion device, which includes:

The electromagnetic signal excitation module is installed in the geological area to be measured, and stimulates and emits transient electromagnetic signals;

The electromagnetic signal receiving device is coupled with the electromagnetic signal excitation module, mounted on the drone, and collects electromagnetic response data;

Stratigraphic model building module, obtain stratigraphic information and build stratigraphic model;

The stratigraphic model correction module is connected to the stratigraphic model building module, and uses the data fitting gradient term, the L1 norm constraint term, and the L2 norm constraint term for correction, and uses the corrected stratigraphic model for inversion.

The third part, the present invention provides a computer-readable storage medium, which stores a computer program, and when the computer program is executed by a processor, realizes the steps of a semi-aviation transient electromagnetic data inversion method based on L1L2 hybrid norm .

Compared with the prior art, the present invention has the following beneficial effects:

(1) The present invention adopts the ground-transmitting transient electromagnetic signal, and uses the UAV to carry the receiving device to obtain the electromagnetic response data, which has the advantages of higher precision, convenient implementation, lower cost, and good safety; it also has the advantages of fast exploration speed , the advantages of exploration that can overcome obstacles. The invention has broad application prospects in the fields of urban underground space detection, ground geological survey, mineral resource exploration and environmental monitoring.

(2) The present invention adopts the L1 norm for the vertical space constraint (regular term) of the resistivity distribution of the measuring point, and the inversion method of the L2 norm for the lateral space constraint. Its objective function includes two constraints, the L1 norm and the L2 norm, which respectively consider the parameter mutual constraints between the layers in the vertical direction of the stratum, that is, the vertical constraints and the lateral constraints between the adjacent measurement points. Compared with the results of the general method, the inversion results are more continuous in the horizontal direction of the formation, and the identification of the formation in the vertical direction is more accurate.

(3) The present invention adopts an adaptive regularization factor to adjust the weight of the constraint term, so that the inversion has higher fitting precision and more stable iteration.

(4) The invention is iteratively stable for the inversion of constrained data, can clearly distinguish the electrical interface of the formation, and makes up for the defects of the current semi-aviation time domain electromagnetic detection data inversion and interpretation system.

To sum up, the present invention has the advantages of simple logic, accuracy and reliability, high detection efficiency, etc., and has high practical value and promotion value in the technical field of geophysical aviation electromagnetic exploration.

Description of drawings

In order to illustrate the technical solutions of the embodiments of the present invention more clearly, the following briefly introduces the accompanying drawings used in the embodiments. It should be understood that the following drawings only show some embodiments of the present invention, and therefore should not be It is regarded as a limitation on the protection scope, and for those skilled in the art, other related drawings can also be obtained according to these drawings without any creative effort.

FIG. 1 is a logic flow diagram of the present invention.

Figure 2 shows the inversion results of groundwater low-resistance anomalies in a certain area in the present invention.

FIG. 3 is a comparison diagram of the theoretical model and inversion results of a certain region in the present invention.

FIG. 4 is an inversion imaging diagram of measured data in a certain area in the present invention.

Detailed ways

In order to make the purpose, technical solutions and advantages of the present application clearer, the present invention will be further described below with reference to the accompanying drawings and examples. The embodiments of the present invention include but are not limited to the following examples. Based on the embodiments in the present application, all other embodiments obtained by those of ordinary skill in the art without creative work shall fall within the protection scope of the present application.

Example 1

As shown in FIG. 1 , this embodiment provides a method and device for inversion of semi-aviation transient electromagnetic data based on L1L2 hybrid norm. Specifically, the device in this embodiment includes: an electromagnetic signal excitation module, an electromagnetic signal receiving device, A stratigraphic model building module and a stratigraphic model correction module; wherein, the electromagnetic signal excitation module is installed in the geological region to be measured, and stimulates and emits transient electromagnetic signals. In addition, the electromagnetic signal receiving device, coupled with the electromagnetic signal excitation module, is mounted on the UAV and collects electromagnetic response data; the formation model building module obtains formation information and constructs the formation model; the formation model correction module is connected with the formation model building module , and use the data fitting gradient term, L1 norm constraint term, L2 norm constraint term for correction, and use the corrected stratigraphic model for inversion.

In this embodiment, the semi-aviation transient electromagnetic data inversion method includes the following steps:

The first step is to use the receiving device to collect the electromagnetic response data of the area to be measured, which also includes the following basic information: the area of the receiving coil, the length of the line source, the number of measurement points of a single measurement line, the total number of measurement lines, the current peak value, and the extraction time. number and channel extraction time; sequentially obtain the x, y coordinates, elevations, and corresponding response values of each time channel of all measuring points involved in the operation.

In the second step, the stratigraphic model of the area to be detected is constructed according to the actual workplace, and the iteration termination conditions are preset. The parameters of the formation model include the total number of formations, the total formation thickness, the initial formation thickness, and the resistivity of any layer. In addition, the iteration termination conditions include the maximum number of iterations and the minimum fitting error.

The third step is to obtain the transverse constraint weighting matrix corresponding to the L2 norm constraint term of the electromagnetic response data; obtain or update the longitudinal constraint matrix corresponding to the L1 norm constraint term of the electromagnetic response data; fit the data to the gradient term, L1 norm term The number constraint term and the L2 norm constraint term constitute the objective function, and after the model increment is derived to be equal to zero, the formula is expanded, and the model increment is obtained by inversion calculation using the expansion formula and the model is corrected. After reaching the iteration termination condition, the inversion imaging image is obtained.

Specifically:

(1) Forward Response Calculation

In the semi-aviation transient electromagnetic survey, the receiving coil mainly receives the magnetic field _Hz component. In this embodiment, the forward response is based on analyzing the horizontal layered vertical component _Hz :

Perform frequency-time conversion to get the rate of change of the magnetic field strength with respect to time:

To obtain the induced electromotive force, its expression is:

In the formula, H _z (w) represents the frequency domain magnetic field response of the vertical component (A/m); w represents the sampling angular frequency (rad/s); I represents the current intensity (A); the center of the wire source is set at the coordinate The origin, extending to L and -L on both sides of the x-axis; R represents the offset distance; r _TE represents the reflection coefficient in TE mode; k ₀ represents the air medium wave number; z represents the coil receiving height (m); J ₁ represents The first-order Bessel function is different from the previous Jacobian matrix; λ represents the integral variable;

In the formula, ReHz(w _j ) represents the real part of Hz(w _j ); N represents the number of turns of the coil; S represents the effective area of the coil (m ² ); μ ₀ represents the magnetic permeability under vacuum (4π×10 ^{− 7} H/m).

The frequency domain response H _z (w) is transformed into the time domain response V _z (t) using a cosine transform.

(2) Inversion imaging method

The semi-airborne transient electromagnetic data inversion method based on hybrid norm is an optimal inversion algorithm based on fitting observation data and obtaining the most ideal constraint model. The basic idea is to divide the underground medium into multi-layered layered structures. , the thickness of each layer is added to the model constraint function, that is, the roughness matrix, and the resistivity of each layer is obtained through calculation, thereby constructing the underground geoelectric structure. Its overall objective function can be summed up as:

expands to

In the _formula : φ( _m ) represents the total objective function; λ represents the regularization factor; φm(m) represents the model constraint objective function; φd(m) represents the observation data objective function; M represents the model vector, which can be expressed as formula (3 ); σ represents resistivity.

M=[σ ₁₁ ,σ ₁₂ …σ _1p …σ _i1 ,σ _i2 ,…σ _ip …σ _n1 ,σ _n2 …σ _np ] ^T (6)

The first item on the right side of formula (5) is the data fitting item, F(M) represents the forward response of the model vector M, W _d is generally a main diagonal matrix, and its elements are the reciprocal of the observed data noise; the second item is the model The vertical constraint term is also called the regular penalty term, W _m is the model vertical constraint matrix, and the model constraint term used in this embodiment is:

That is, the unit first-order difference operator; the third item is the lateral constraint item of the model. In this paper, four adjacent measurement point constraint items are considered, W _{h, i} are the transverse constraint matrix, q represents the number of constraint items, and the same The first-order difference operator is expressed as:

Since the L1 norm is non-derivable, the iteration will be unstable, and the direction of the fastest descent of the regular term cannot be obtained. A small value ξ can be taken, and formula (5) can be rewritten as follows:

Using iterative reweighted least squares (IRLS), the L1 regular term can be rewritten:

This embodiment improves this algorithm, and expresses the L1 regular term as:

Among them, V is the diagonal weighting matrix, and the elements on the main diagonal are:

表示纵向约束项粗糙度矩阵转置；W_m表示纵向约束粗糙度矩阵。Among them, V represents the diagonal weighting matrix, _xi represents the weighted model vector element ξ represents a small value; X ^T represents the transpose of the weighted model vector matrix; X represents the weighted model vector;

represents the transposition of the longitudinal constraint item roughness matrix; W _m represents the longitudinal constraint roughness matrix.

The weight matrix D in the L2 norm constraint is expanded as:

Among them, r represents the diagonal element of the transverse constraint weighting matrix; n represents the nth measurement point; p represents the pth adjacent point of the measurement point n.

这里x和y是所计算的测点到参与计算的相临点之间的相对坐标，这里计算的两点之间的距离取倒数。若不考虑加权，r_i取对角元素为1。最后计算

项并加入目标函数。Different diagonal weighting elements can be considered according to the actual situation, when distance weighting is used

Here x and y are the relative coordinates between the calculated measuring point and the adjacent points involved in the calculation, and the distance between the two points calculated here is the reciprocal. If weighting is not considered, _ri takes the diagonal element as 1. final calculation

term and add the objective function.

Take the derivation of the model correction amount Δm _k , and set the derived formula equal to 0 to obtain the solution formula of the model correction amount:

表示雅可比矩阵转置；

表示数据拟合项粗糙度矩阵转置；W_d表示数据拟合项粗糙度矩阵；d^obs表示实测响应数据；F(M_k)表示模型向量正演计算项；λ₁表示第一正则化因子；

表示纵向约束项粗糙度矩阵转置；V_k表示纵向约束项加权矩阵；W_m表示纵向约束项粗糙度矩阵；M_k表示模型向量；λ₂表示第二正则化因子；

表示横向约束粗糙度矩阵转置；D_i表示横向约束加权矩阵；W_h,i表示横向约束粗糙度矩阵。Among them, the first term on the right is the data fitting gradient term, the second term is the L1 norm vertical space constraint term, and the third term is the L2 norm lateral space constraint term (regularized gradient term). in addition,

represents the transpose of the Jacobian matrix;

represents the transposition of the data fitting item roughness matrix; W _d represents the data fitting item roughness matrix; d ^obs represents the measured response data; F(M _k ) represents the model vector forward calculation term; λ ₁ represents the first regularization factor ;

represents the longitudinal constraint term roughness matrix transposition; V _k represents the longitudinal constraint term weighting matrix; W _m represents the longitudinal constraint term roughness matrix; M _k represents the model vector; λ ₂ represents the second regularization factor;

represents the transposition of the transverse constraint roughness matrix; D _i represents the transverse constraint weight matrix; W _h,i represents the transverse constraint roughness matrix.

In the fourth step, the adaptive regularization factor is used to adjust the weight of the L1 norm constraint term and the L2 norm constraint term in the inversion process; and the modified stratigraphic model is used for inversion, and the inversion imaging image is obtained.

This embodiment provides that the regularization factor adopts an adaptive adjustment method, the data gradient term is used as the numerator, and the five model constraint gradient terms are used as the denominator. A parameter a is used to adjust the weights of λ ₁ and λ ₂ in the objective function respectively, which are used to adjust the weights of the vertical and horizontal constraints in the iterative process. The specific formula is as follows:

表示拟合数据梯度项；

表示纵向约束梯度项；

表示横向约束梯度项。Among them, a represents the weight ratio parameter;

represents the fitted data gradient term;

represents the longitudinal constraint gradient term;

represents the laterally constrained gradient term.

Example 2

This embodiment provides a computer-readable storage medium, which stores a computer program, and when the computer program is executed by a processor, implements the steps of the semi-aviation transient electromagnetic data inversion method in Embodiment 1.

Example 3

As shown in Figure 2 to Figure 3, this embodiment takes the data inversion of an underground low-resistance model with a noise of 5% as an example:

The semi-aviation transient electromagnetic technology of unmanned aerial vehicle is used to simulate the search for groundwater, and a groundwater low-resistance layer (thickness 60m) of 10Ω·m is set in a high-resistance layer (thickness 240m) with a background resistivity of 100Ω·m. In the inversion, the thickness of each layer is set with an exponential growth factor of 1.05, and the initial model of the inversion is a uniform half-space of 50Ω·m. The inversion results are shown in Figure 2. It can be seen from the figure that the inversion result is generally divided into three layers of "high-low-high", in which the high resistance is distributed between 100Ω·m-110Ω·m, and the inversion resistivity is about 10Ω·m. Very close to the real model, the UAV semi-aviation transient electromagnetic inversion method can effectively invert low-resistance anomalies, inversion of noisy data is stable, and the electrical interface is clearly depicted. This shows that the semi-aviation transient electromagnetic is sensitive to low-resistance anomalies.

Example 4

As shown in Figure 4, in this embodiment, the measured data of a certain place is used for analysis, and after simple processing of the flight detection data, the uniform half-space is used as the initial model, and the results of inversion of multiple survey lines are shown in the figure. The overall resistivity shows a low-high-low distribution trend, and the overall inversion effect is better. The inversion speed of this method is faster than that of the one-dimensional single-point inversion. It can be seen from the figure that the continuity of the inversion results is also better. Well, the characterization of the electrical interface is very clear, and the inversion iteration process for noisy data is more stable.

The above-mentioned embodiments are only the preferred embodiments of the present invention, and are not intended to limit the protection scope of the present invention. Any changes made by adopting the design principles of the present invention and non-creative work on this basis shall belong to the protection scope of the present invention. within.

Claims (9)

1. The semi-aviation transient electromagnetic data inversion method based on the L1L2 mixed norm adopts a ground transmitted transient electromagnetic signal and utilizes an unmanned aerial vehicle carrying and receiving device to obtain electromagnetic response data, and is characterized by comprising the following steps of:

collecting electromagnetic response data of a region to be detected by using a receiving device;

constructing a stratum model of a region to be detected, and presetting an iteration termination condition;

obtaining a transverse constraint weighting matrix corresponding to an L2 norm constraint term of the electromagnetic response data;

solving or updating a longitudinal constraint matrix corresponding to the L1 norm constraint term of the electromagnetic response data;

forming a target function by the data fitting gradient term, the L1 norm constraint term and the L2 norm constraint term, and solving an expansion formula after the incremental derivation of the stratum model is equal to zero;

obtaining stratum model increment and correcting the stratum model by using an expansion formula for inversion calculation, and obtaining an inversion imaging image after an iteration termination condition is reached; adjusting the weights of the L1 norm constraint term and the L2 norm constraint term in the inversion process by adopting an adaptive regularization factor;

fitting a gradient term, an L1 norm constraint term and an L2 norm constraint term according to the data to form an objective function, performing model increment derivation equal to zero and expanding, and then correcting the stratum model, wherein the expression is as follows:

is unfolded into

Wherein,

representing the fitting of the data to the gradient term,

represents the L1 norm constraint term;

represents the L2 norm constraint term;

representing a Jacobian matrix transpose;

representing data fitting item roughness matrix transposition; w _d Representing a data fitting term roughness matrix; d is a radical of ^obs Representing measured response data; f (M) _k ) Forward calculation terms representing model vectors; lambda ₁ Representing a first regularization factor;

a roughness matrix representing longitudinal constraint terms; v _k Representing a longitudinal constraint weighting matrix; w _m A roughness matrix representing longitudinal constraint terms; m _k Representing a model vector; lambda ₂ Representing a second regularization factor;

representing a transverse constraint term roughness matrix; d _i A weighting matrix representing the laterally constrained terms; w is a group of _h,i Roughness moment representing transverse constraint termArraying; m represents a model vector; f (M) represents the forward response of the model vector M; q represents the number of constraint terms.

2. The L1L2 mixed norm based semi-aviation transient electromagnetic data inversion method of claim 1, further comprising:

acquiring basic information of the receiving device; the basic information comprises the area of a receiving coil, the flying height, the number of track extraction time and the track extraction time;

acquiring basic information of a transmitting device; the basic information of the emitting device comprises a line source length and a current peak value;

and acquiring corresponding values corresponding to the number of measurement points of a single measurement line in the area to be measured, the total number of measurement lines, the coordinates of any measurement point, the offset distance, the elevation and the time channel.

3. The L1L2 mixed norm-based semi-aviation transient electromagnetic data inversion method as recited in claim 1, wherein the iteration termination condition comprises a maximum number of iterations and a minimum fitting error.

4. The method for semi-aeronautical transient electromagnetic data inversion based on the L1L2 mixed norm according to claim 1, 2 or 3, characterized in that the parameters of the formation model include total number of formation layers, total formation thickness, starting formation thickness, resistivity of any layer.

5. The method for inverting semi-aviation transient electromagnetic data based on L1L2 mixed norm as claimed in claim 1, 2 or 3 wherein the expression of the transverse constraint weighting matrix is:

wherein r represents the diagonal element of the transverse constraint weighting matrix, n represents the nth measuring point, and p represents the p adjacent point of the n measuring point.

6. The L1L2 mixed-norm based semi-aviation transient electromagnetic data inversion method as claimed in claim 5, wherein the expression of the regular term of the L1 norm is:

wherein V represents a diagonal weighting matrix; x is the number of _i Representing a small value of model vector element xi; x ^T Representing a transpose of a weighted model matrix; x weights the model matrix;

transposing the roughness matrix; w is a group of _m Representing a roughness matrix; m represents a model vector;

the main diagonal elements of the diagonal weighting matrix are:

7. the semi-aviation transient electromagnetic data inversion method based on the L1L2 mixed norm as claimed in claim 1, characterized in that the weights of the L1 norm constraint term and the L2 norm constraint term in the inversion process are adjusted by using an adaptive regularization factor, and the expression is as follows:

wherein, a represents a weight proportion parameter;

representing fitted data gradient terms;

representing longitudinal constrained gradientsAn item;

representing the transverse constraint gradient term.

8. A semi-airborne transient electromagnetic data inversion apparatus, comprising:

the electromagnetic signal excitation module is arranged in a geological region to be detected and used for exciting and transmitting transient electromagnetic signals;

the electromagnetic signal receiving device is coupled with the electromagnetic signal excitation module, carried on the unmanned aerial vehicle and used for collecting electromagnetic response data;

the stratum model building module is used for obtaining stratum information and building a stratum model;

the stratum model correction module is connected with the stratum model construction module, corrects the stratum model by adopting a data fitting gradient term, an L1 norm constraint term and an L2 norm constraint term, and performs inversion by using the corrected stratum model; adjusting the weight of the L1 norm constraint term and the L2 norm constraint term in the inversion process by adopting a self-adaptive regularization factor;

fitting a gradient term, an L1 norm constraint term and an L2 norm constraint term according to the data to form an objective function, solving the incremental derivative of the model to be equal to zero, expanding, and then correcting the stratum model, wherein the expression is as follows:

is unfolded into

Wherein,

representing the fitting of the data to the gradient term,

represents the L1 norm constraint term;

represents the L2 norm constraint term;

representing a Jacobian matrix transpose;

representing data fitting item roughness matrix transposition; w is a group of _d Representing a data fitting item roughness matrix; d ^obs Representing measured response data; f (M) _k ) Forward calculation terms representing model vectors; lambda [ alpha ] ₁ Representing a first regularization factor;

a roughness matrix representing longitudinal constraint terms; v _k Representing a longitudinal constraint weighting matrix; w is a group of _m A roughness matrix representing longitudinal constraint terms; m is a group of _k Representing a model vector; lambda [ alpha ] ₂ Representing a second regularization factor;

representing a transverse constraint term roughness matrix; d _i A weighting matrix representing the laterally constrained terms; w _h,i A roughness matrix representing lateral constraint terms; m represents a model vector; f (M) represents the forward response of the model vector M; q represents the number of constraint terms.

9. A computer readable storage medium storing a computer program, wherein the computer program, when executed by a processor, performs the steps of the L1L2 mixed norm based semi-airborne transient electromagnetic data inversion method of any one of claims 1 to 7.

Images