CN104008249A - Dynamic contour model based surface nuclear magnetic resonance inversion method - Google Patents

Abstract

The invention discloses a dynamic contour model based surface nuclear magnetic resonance inversion method. The dynamic contour model based surface nuclear magnetic resonance inversion method comprises setting the number of water content layers in a detection area and initializing the thickness and a water content value of every layer in the SNMR (Surface Nuclear Magnetic Resonance) layered inversion solution process; dividing a sub-water-content vertical distribution graph into MN micro elements to meet a matrix equation of a dynamic model; performing iteration solution on the matrix equation and performing dynamic adjustment on the thickness and the water content value of every water content layer in the iteration process to search an optimal solution which meets the matrix equation. The dynamic adjustment is constantly performed on a contour of the vertical distribution graph of the water content values in the detection area in the integral solution process and accordingly the dynamic contour model based surface nuclear magnetic resonance inversion method is called the dynamic contour model and an SGD (Stochastic Gradient Descent) method is designed to solve the model. According to the dynamic contour model based surface nuclear magnetic resonance inversion method, the convergence speed is high, the inversion result is high in accuracy and stable, the performance is excellent in comparison with a regularization method, and the high-accuracy solution on an SNMR inversion problem can be achieved.

Description

Ground nuclear magnetic resonance inversion method based on contours model

Technical field

The present invention relates to ground nuclear magnetic resonance (Surface Nuclear Magnetic Resonance is called for short SNMR) field, be specifically related to a kind of ground nuclear magnetic resonance inversion method based on contours model.

Background technology

Ground nuclear magnetic resonance (Surface Nuclear Magnetic Resonance, abbreviation SNMR) technology is the current geophysical prospecting method of unique a kind of direct water detection in the world, and this technology has obtained certain application in fields such as Underground water, archaeology, groundwater contamination detections.

In recent years, along with expert and scholars' further investigation gradually, SNMR technology has obtained further perfect.Inversion Calculation water percentage is the key link in this technical research process, and the accuracy of inversion result and the operation efficiency of refutation process are the key indexs of weighing inversion algorithm performance.Wherein, one dimension Forward And Inverse Problems is comparatively ripe, has in succession published out multiple efficient algorithm, as improved simulated annealing inverting, has improved degree of stability and the speed of convergence of existing inversion algorithm; QT inversion algorithm, utilize whole sampling numbers certificates that each excitation pulse square is corresponding to carry out inverting, fully excavated reception signal message, improved to a certain extent inversion accuracy, but, owing to receiving signal, present approximate exponential damping, late period, Signal-to-Noise was very low, and the method is only applicable to high s/n ratio environment; Also having afterwards scholar to adopt integration gate technique to receive signal, and improved the precision of each sampling number certificate, and carried out full attenuation inverting, is a kind of improvement to QT inverting.

SNMR technology is aspect two dimensional inversion, Boucher, Girard and Legchenko etc. have studied in two dimensional cross-section direction E0-q curve with the variation tendency of underground water-bearing structure, but they have only done qualitative examination to two dimensional inversion, do not provide concrete two dimensional inversion formula.Legchenko etc. have done certain research to 3-d inversion, although can be finally inversed by three dimensions the water-bearing structure of model, but because the size of mesh opening of setting at three dimensions is larger, can only go out underground water-bearing structure by "ball-park" estimate, its resolution of inversion has much room for improvement.Because two dimension, 3-d inversion algorithm exist the problems such as operand is large, variable number to be solved is many, non-linear, unique commercial version Inversion Software NUMISPLUS still adopts one-dimensional inversion in the world at present.

Ground nuclear magnetic resonance inverting energy is abstract is a problem that solves ill-condition matrix equation, in existing document, the regularization methods that adopt solve it more, but, existing Standard Regularization method and improving one's methods are on its basis deposited deficiency both ways, respectively: regularization parameter is chosen difficulty; The approximate solution of inverse problem is forced to add slickness constraint, weakened the difference of water cut between adjacent geological stratification.

Summary of the invention

Technical matters to be solved by this invention is to provide a kind of ground nuclear magnetic resonance inversion method based on contours model, and it can be realized the high precision of ground nuclear magnetic resonance inverse problem is solved.

For addressing the above problem, the present invention is achieved by the following technical solutions:

A ground nuclear magnetic resonance inversion method based on contours model, comprises the steps:

Step 1, the contours model of the ground nuclear magnetic resonance inverting of foundation,

An＝E ①

In formula, A is kernel matrix; N is water cut vector to be solved, m _ninfinitesimal number for water-bearing zone in model; E is ground nuclear magnetic resonance signal initial amplitude value vector, $E={(E\left(q_1\right),E\left(q_2\right),...,E\left(q_{I_N}\right))}^T,$ I _nfor excitation pulse square number;

Step 2, the moisture number of plies in initialization search coverage is L, thickness and the moisture content value of each layer are respectively H _l=(h ₁, h ₂..., h _l) ^tand N _l=(n ₁, n ₂..., n _l) ^t; Maximum iteration time is N _max; Primary iteration number of times is N _i=0; Inversion accuracy Threshold is θ;

Step 3, is divided into M by each sub-water cut vertical distribution figure of a current iteration _nindividual infinitesimal, each infinitesimal thickness is △ h, and generates the kernel matrix A of contours model;

Step 4, calculates the searching route of future generation of contours model, and the searching route of renewal is considered as to the searching route of a current iteration; Wherein searching route of future generation more new formula be

$\left\{\begin{array}{c}{h}_{k+1}^{*}=h_k+w_h\mathrm{\&Delta;}{u}_k\\ h_{k+1}=\left\{\begin{array}{c}1,h_{k+1}^{*}<1\\ h_{k+1}^{*},1\&le;h_{k+1}^{*}\&le;20\\ 20,h_{k+1}^{*}>20\end{array}\right.\\ n_{k+1}^{*}=n_k+w_n\mathrm{\&Delta;}{v}_k\\ n_{k+1}=\left\{\begin{array}{c}0,n_{k+1}^{*}<0\\ n_{k+1}^{*},0\&le;n_{k+1}^{*}\&le;1\\ 1,n_{k+1}^{*}>1\end{array}\right.\end{array}\right.$ ②

In formula, h _kand n _kfor thickness and the water cut of this each layer of iteration, h _k+1and n _k+1thickness and water cut for each layer of next iteration; w _hand w _nbe respectively the thickness of each layer and the mobile weight of moisture content value; △ u _kwith △ v _kbe respectively the thickness of each layer and the moving direction of moisture content value vector;

Above-mentioned △ u _kfor the k row in each layer thickness moving direction matrix △ u, △ v _kfor the k row in each layer of water cut moving direction matrix △ v, $\mathrm{\&Delta;}{v}_k={(\mathrm{\&Delta;}{v}_{1k},{\mathrm{\&Delta;v}}_{2k},...,{\mathrm{\&Delta;v}}_{M_{N}k})}^T;$ Wherein

In formula, the element △ u in each layer thickness moving direction matrix △ u _ij∈ { x|-1,0,1}, the element △ v in each layer of water cut moving direction matrix △ v _ij∈ { x|-1,0,1}, and each iteration △ u and all generations at random of △ v.

Step 5, the fitness function of calculating contours model, upgrades optimum solution;

Fitness function is defined as:

f _min＝||An _k+1-E ₀|| ₂ ⑤

In formula, A is kernel matrix, n _k+1for sub-water cut vector of a current generation, E ₀for initial amplitude value corresponding to each excitation pulse square in actual measurement ground nuclear magnetic resonance signal;

Calculate the optimum solution in each sub-water cut vector in a current generation; If when containing layer thickness in current generation optimum solution and being the water-bearing zone of △ h rice, delete this layer, and the layer thickness of its last layer is increased to △ h rice, the moisture number of plies L=L-1 in update detection region; When if ground floor layer thickness is △ h rice, itself and the second layer are merged;

Step 6: judge whether iteration stops; If iterations reaches maximum iteration time N _max, or current optimal-adaptive degree functional value is less than inversion accuracy threshold value θ, stops iteration, and inversion result is current optimum solution; Otherwise, return to step 4, carry out next iteration.

In described step 4, the mobile weight w of the thickness of each layer and moisture content value _hand w _nadopt variable step-size search,

The mobile weight w of the thickness of each layer _hfor

w _h＝(1-N _I/N _max)△w _h ⑥

In formula, N _maxfor maximum iteration time, N _ifor current iteration number of times, △ w _hfixed step size for thickness.

The mobile weight w of the water cut of each layer _nfor

w _n＝(1-N _I/N _max)△w _n ⑦

In formula, N _maxfor maximum iteration time, N _ifor current iteration number of times, △ w _nfixed step size for water cut.

The present invention, in SNMR stratiform inverting solution procedure, first sets the moisture number of plies in search coverage, and the thickness of each layer of initialization and moisture content value; Then, be divided into M _nindividual infinitesimal, makes it meet the matrix equation of dynamic model; Finally, to matrix equation iterative, in iterative process, the thickness in each water-bearing zone and moisture content value are dynamically adjusted, to search for the optimum solution that meets matrix equation.Whole solution procedure, constantly dynamically adjusts the profile of the vertical distribution figure of moisture content value in search coverage, therefore the method is called " contours model ", and has designed random gradient descent method (SGD) and solves this model; The present invention has fast convergence rate, inversion result precision is high and stable, and its performance is better than regularization method, can realize the high precision of SNMR inverse problem is solved.

Accompanying drawing explanation

Fig. 1 is ground nuclear magnetic resonance inversion method based on contours model and ground nuclear magnetic resonance inversion method error of fitting convergence curve figure in solution procedure based on regularization model.

Fig. 2 is initial amplitude value E in the ground nuclear magnetic resonance inversion method based on contours model ₀correlation curve with actual value, measured value and the inversion result value of excitation pulse square q.

Fig. 3 is the water cut vertical distribution contrast figure of realistic model and inversion result in the ground nuclear magnetic resonance inversion method based on contours model.

Embodiment

In order to verify that contours model of the present invention is for the validity of SNMR inverting, model water-bearing zone vertical distribution model, then by forward model, calculate the theoretical value of the initial amplitude that a series of excitation pulse squares are corresponding, on the basis of initial amplitude theoretical value, add the noise of certain signal to noise ratio (S/N ratio) as observation signal, finally adopt respectively ACM model (Advanced Computer Modeling, advanced computer model) and TRM model (the Tikhonov regularization method, regularization model) carry out inverting.Meanwhile, for implementation model contrast, adopt random gradient descent method to ask for respectively the optimum solution of above-mentioned two models.Definition signal to noise ratio snr formula is:

$\mathrm{SNR}=20{\log }_{10}\frac{{\left|\right|d\left|\right|}_2}{{\left|\right|d-d^{*}\left|\right|}_2},$

Wherein, d ^*be respectively observation signal and theory signal with d.E corresponding to each excitation pulse square obtaining with inverting ₀the error of fitting of value is evaluated inversion result performance, and error of fitting formula is:

$f_E={\left|\right|E_0-E_0^{*}\left|\right|}_2=\sqrt{\underset{i=1}{\overset{M_N}{\mathrm{\&Sigma;}}}{(E_{0i}-E_{0i}^{*})}^2},$

Wherein, E _0iwith be respectively the initial amplitude value corresponding to each excitation pulse square of theoretical model and inversion result.

Suppose that theoretical model background resistivity is 200 Europe rice, geomagnetic field intensity is 48000nT, magnetic declination is 40 °, maximum excitation pulse square is 20As, employing radius is that the round coil of 100 meters transmits and receives signal, in vertical depth 0-150 rice interval, contain three water-bearing zones, be respectively: between vertical depth 20-35 rice, water cut is 20%; Between vertical depth 60-75 rice, water cut is 50%; Between vertical depth 105-120 rice, water cut is 30%, and other degree of depth water cut are 1%.During calculating, first adopt theoretical model to calculate the initial amplitude value of the SNMR signal that each excitation pulse square is corresponding, and on the basis of theoretical initial amplitude value, to add signal to noise ratio (S/N ratio) be that the random noise of 20dB is as observation signal.Adopt respectively ACM model and TRM model to carry out inverting to above-mentioned moisture tectonic structure.

The ground nuclear magnetic resonance inversion method that the present invention is based on contours model, comprises the steps:

Step 1: the contours model of the ground nuclear magnetic resonance inverting of foundation.

In horizontal layer model, SNMR (ground nuclear magnetic resonance) signal initial amplitude discretization model is:

$E_0\left(q_i\right)=\underset{j=1}{\overset{\&infin;}{\mathrm{\&Sigma;}}}K(q_i,z_j)n_j{\mathrm{\&Delta;z}}_j,---\left(1\right)$

Wherein, q _ibe i excitation pulse square, n _jwith △ z _jbe respectively water cut and the thickness of j layer, K (q _i, z _j) be the kernel function after discretize under cylindrical-coordinate system;

Wherein, M ₀for the magnetization; b _{1 ⊥}=u ₀h _{1 ⊥}/ I ₀, H _{1 ⊥}for exciting magnetic field intensity perpendicular to the component of terrestrial magnetic field, u ₀be permeability of vacuum, it equals 4 π * 10 ^-7h/m, I ₀for reflected current intensity; θ _imnjfor differentiation element (r _m, Φ _n, z _j) locate pull angle down; △ r _mwith △ Φ _nbe respectively (the r of differential unit _m, Φ _n) locate radially and lateral length; M and N are respectively radially and the horizontal number of differentiation element;

In ground nuclear magnetic resonance is surveyed, by launching the excitation pulse square q of different sizes _i, to excite hydrogen ion in different depth place water, (1) formula is converted into matrix equation:

An＝E, (3)

Wherein, $E={(E\left(q_1\right),E\left(q_2\right),...,E\left(q_{I_N}\right))}^T,$ I _nfor excitation pulse square number; A _ij=K (q _i, z _j) △ z _j; m _ninfinitesimal number for water-bearing zone in model.

The concrete enforcement solution procedure of the contours model of the ground nuclear magnetic resonance inverting of above-mentioned foundation is as follows:

Step 2: algorithm initialization.Suppose that the maximum probe degree of depth is H _max, the kernel matrix A of growth equation (3), adopting layer thickness is the infinitesimal of △ h, in iterative optimum solution process, only need to calculate kernel matrix A one time like this, can greatly improve the solution efficiency of algorithm, in the present embodiment, △ h=1m; The moisture number of plies in initialization search coverage is L, and thickness and the moisture content value of each layer are respectively H _l=(h ₁, h ₂..., h _l) ^tand N _l=(n ₁, n ₂..., n _l) ^t; Maximum iteration time is N _max; Primary iteration number of times is N _i=0; The sub-number of individuals of every generation is M _n; Inversion accuracy threshold value is θ, and in the present embodiment, θ value is 10 ^-8.The value of layer thickness △ h is 1 meter.

Step 3: calculate searching route of future generation.Searching route of future generation more new formula is:

$\left\{\begin{array}{c}{h}_{k+1}^{*}=h_k+w_h\mathrm{\&Delta;}{u}_k\\ h_{k+1}=\left\{\begin{array}{c}1,h_{k+1}^{*}<1\\ h_{k+1}^{*},1\&le;h_{k+1}^{*}\&le;20\\ 20,h_{k+1}^{*}>20\end{array}\right.\\ n_{k+1}^{*}=n_k+w_n\mathrm{\&Delta;}{v}_k\\ n_{k+1}=\left\{\begin{array}{c}0,n_{k+1}^{*}<0\\ n_{k+1}^{*},0\&le;n_{k+1}^{*}\&le;1\\ 1,n_{k+1}^{*}>1\end{array}\right.\end{array}\right.---\left(4\right)$

Wherein, w _hand w _nbe respectively the thickness of each layer and the mobile weight of moisture content value, △ u _kwith △ v _kbe respectively the thickness of each layer and the moving direction of moisture content value vector.

Each layer thickness moving direction matrix is wherein, △ u _ij∈ x|-1,0,1}, and each iteration △ u all generates at random.△ u vector is turned to △ u=(△ u ₁, △ u ₂..., △ u _l), wherein, ${\mathrm{\&Delta;u}}_k={({\mathrm{\&Delta;u}}_{1k},{\mathrm{\&Delta;u}}_{2k},...,{\mathrm{\&Delta;u}}_{M_{n}k})}^T.$

Each layer of water cut moving direction matrix is wherein, △ v _ij∈ x|-1,0,1}, and each iteration △ v all generates at random.△ v vector is turned to △ v=(△ v ₁, △ v ₂..., △ v _l), wherein, $\mathrm{\&Delta;}{v}_k={(\mathrm{\&Delta;}{v}_{1k},{\mathrm{\&Delta;v}}_{2k},...,{\mathrm{\&Delta;v}}_{M_{N}k})}^T.$

For convergence speedup speed, guarantee computational accuracy simultaneously, adopt variable step-size search, w _h=(1-N _i/ N _max) △ w _h, △ w _hfor fixed step size, in this example, value is 4, w _n=(1-N _i/ N _max) △ w _n, △ w _nfor fixed step size, in this example, value is 0.05.

Step 4: water cut vertical distribution figure infinitesimalization is processed.To work as each sub-water cut vertical distribution figure of former generation and be divided into M _nindividual infinitesimal, each infinitesimal thickness is △ h, makes it meet equation (3).

Step 5: calculate fitness function, upgrade optimum solution.Fitness function is defined as:

f _min＝||An _k+1-E ₀|| ₂ (5)

Wherein, A is the kernel matrix of trying to achieve in step 2, n _k+1for a sub-water cut vector of a current generation of trying to achieve in step 4, E ₀for initial amplitude value corresponding to each excitation pulse square in actual measurement SNMR signal.Calculate the optimum solution in each sub-water cut vector in a current generation.

If when containing layer thickness in current generation optimum solution and being the water-bearing zone of △ h rice, delete this layer, and the layer thickness of its last layer is increased to △ h rice, the moisture number of plies L=L-1 in update detection region; When if ground floor layer thickness is △ h rice, itself and the second layer are merged.

Step 6: judge whether iteration stops.If iterations reaches N _max, or current optimal-adaptive degree functional value is less than θ, stops iteration, and inversion result is current optimum solution; Otherwise, return to execution step 3.

Fig. 1 is ground nuclear magnetic resonance inversion method based on contours model and ground nuclear magnetic resonance inversion method error of fitting convergence curve figure in solution procedure based on regularization model.Fig. 2 is initial amplitude value E in the ground nuclear magnetic resonance inversion method based on contours model ₀correlation curve with actual value, measured value and the inversion result value of excitation pulse square q.Fig. 3 is the water cut vertical distribution contrast figure of realistic model and inversion result in the ground nuclear magnetic resonance inversion method based on contours model.Two kinds of final fitting precisions of model are similar as can be seen from Figure 1, but the speed of convergence of ACM model is faster than TRM model, and ACM model only needs 20 iteration just can obtain stable solution, and TRM model needs 200 iteration.As can be seen from Figure 3 the inversion result of ACM model can present three water-bearing zones corresponding with theoretical model, does not have false water-bearing zone, although the inversion result error of deep regions increases to some extent, this error is in the scope of physical prospecting permission; And the distribution of the inversion result water-bearing zone of TRM model is chaotic, differ greatly with theoretical model.

Above-described embodiment, is only the specific case that object of the present invention, technical scheme and beneficial effect are further described, and the present invention is not defined in this.All any modifications of making within scope of disclosure of the present invention, are equal to replacement, improve etc., within being all included in protection scope of the present invention.

Claims (2)

1. the ground nuclear magnetic resonance inversion method based on contours model, is characterized in that comprising the steps:

Step 1, the contours model of the ground nuclear magnetic resonance inverting of foundation,

An＝E ①

In formula, A is kernel matrix; N is water cut vector to be solved, m _ninfinitesimal number for water-bearing zone in model; E is ground nuclear magnetic resonance signal initial amplitude value vector, $E={(E\left(q_1\right),E\left(q_2\right),...,E\left(q_{I_N}\right))}^T,$ I _nfor excitation pulse square number;

Step 2, the moisture number of plies in initialization search coverage is L, thickness and the moisture content value of each layer are respectively H _l=(h ₁, h ₂..., h _l) ^tand N _l=(n ₁, n ₂..., n _l) ^t; Maximum iteration time is N _max; Primary iteration number of times is N _i=0; Inversion accuracy Threshold is θ;

Step 3, is divided into M by each sub-water cut vertical distribution figure of a current iteration _nindividual infinitesimal, each infinitesimal thickness is △ h, and generates the kernel matrix A of contours model;

Step 4, calculates the searching route of future generation of contours model, and the searching route of renewal is considered as to the searching route of a current iteration; Wherein searching route of future generation more new formula be

$\left\{\begin{array}{c}{h}_{k+1}^{*}=h_k+w_h\mathrm{\&Delta;}{u}_k\\ h_{k+1}=\left\{\begin{array}{c}1,h_{k+1}^{*}<1\\ h_{k+1}^{*},1\&le;h_{k+1}^{*}\&le;20\\ 20,h_{k+1}^{*}>20\end{array}\right.\\ n_{k+1}^{*}=n_k+w_n\mathrm{\&Delta;}{v}_k\\ n_{k+1}=\left\{\begin{array}{c}0,n_{k+1}^{*}<0\\ n_{k+1}^{*},0\&le;n_{k+1}^{*}\&le;1\\ 1,n_{k+1}^{*}>1\end{array}\right.\end{array}\right.$ ②

In formula, h _kand n _kfor thickness and the water cut of this each layer of iteration, h _k+1and n _k+1thickness and water cut for each layer of next iteration; w _hand w _nbe respectively the thickness of each layer and the mobile weight of moisture content value; △ u _kwith △ v _kbe respectively the thickness of each layer and the moving direction of moisture content value vector;

Above-mentioned △ u _kfor the k row in each layer thickness moving direction matrix △ u, △ v _kfor the k row in each layer of water cut moving direction matrix △ v, $\mathrm{\&Delta;}{v}_k={(\mathrm{\&Delta;}{v}_{1k},{\mathrm{\&Delta;v}}_{2k},...,{\mathrm{\&Delta;v}}_{M_{N}k})}^T;$ Wherein

In formula, the element △ u in each layer thickness moving direction matrix △ u _ij∈ { x|-1,0,1}, the element △ v in each layer of water cut moving direction matrix △ v _ij∈ { x|-1,0,1}, and each iteration △ u and all generations at random of △ v.

Step 5, the fitness function of calculating contours model, upgrades optimum solution;

Fitness function is defined as:

f _min＝||An _k+1-E ₀|| ₂ ⑤

In formula, A is kernel matrix, n _k+1for sub-water cut vector of a current generation, E ₀for initial amplitude value corresponding to each excitation pulse square in actual measurement ground nuclear magnetic resonance signal;

Calculate the optimum solution in each sub-water cut vector in a current generation; If when containing layer thickness in current generation optimum solution and being the water-bearing zone of △ h rice, delete this layer, and the layer thickness of its last layer is increased to △ h rice, the moisture number of plies L=L-1 in update detection region; When if ground floor layer thickness is △ h rice, itself and the second layer are merged;

Step 6: judge whether iteration stops; If iterations reaches maximum iteration time N _max, or current optimal-adaptive degree functional value is less than inversion accuracy threshold value θ, stops iteration, and inversion result is current optimum solution; Otherwise, return to step 4, carry out next iteration.

2. the ground nuclear magnetic resonance inversion method based on contours model according to claim 1, is characterized in that, in described step 4, and the mobile weight w of the thickness of each layer and moisture content value _hand w _nadopt variable step-size search,

The mobile weight w of the thickness of each layer _hfor

w _h＝(1-N _I/N _max)△w _h ⑥

In formula, N _maxfor maximum iteration time, N _ifor current iteration number of times, △ w _hfixed step size for thickness.

The mobile weight w of the water cut of each layer _nfor

w _n＝(1-N _I/N _max)△w _n ⑦

In formula, N _maxfor maximum iteration time, N _ifor current iteration number of times, △ w _nfixed step size for water cut.