CN104007478B - A kind of ground nuclear magnetic resonance inversion method based on harmonic search algorithm - Google Patents

Abstract

The present invention discloses a kind of ground nuclear magnetic resonance inversion method based on harmonic search algorithm, it first establishes constraint R TLS inverse model, under conditions of the initial amplitude value sequence of background area resistivity Distribution Value sequence and measurement signal all exists error, improve inversion accuracy；After propose harmonic search algorithm that IHS improves to solve problem that this model conversation is the nonlinear optimization by constraint by the derivation of equation.IHS algorithm makes standard HS algorithm early stage have good ergodic, the later stage solves and has high accuracy characteristic；And due to the restriction of constraints, this IHS algorithm remains to when inversion matrix equation is underdetermined equation accurately solve, eliminate the restriction that in inverse model, water-bearing layer maximum is divided the number of plies.

Description

A kind of ground nuclear magnetic resonance inversion method based on harmonic search algorithm

Technical field

The present invention relates to ground nuclear magnetic resonance technical field, be specifically related to a kind of ground based on harmonic search algorithm nuclear-magnetism Resonance inversion method.

Background technology

Ground nuclear magnetic resonance (Surface Nuclear Magnetic Resonance is called for short SNMR) technology has operation Simply, the characteristic such as informative and quantitative inversion, and be the geophysical prospecting method of the most unique a kind of direct look for water at present. The principle of SNMR technology is: by launching a series of different size of pulse squares to underground, to excite water in the range of different depth Middle proton so that it is energy jump occurs, when transmitter stops emission signal, has occurred the proton of energy jump to affranchise State, reception antenna can sense that this process proton releases energy number, receiver gathers the signal of reception antenna sensing, logical Cross docking collection of letters inversion interpretation, can determine that the hydrogeology structure in investigative range.This technology Underground water, archaeology, The fields such as underground water pollution detection have obtained certain application.

In recent years, under the effort research of many experts and scholar, SNMR technology has obtained further perfect.Wherein, In terms of SNMR signal inverting, one-dimensional Forward And Inverse Problems is the most ripe, and is furtheing investigate two dimension, three-dimensional FORWARD AND INVERSE PROBLEMS model. SNMR signal inverting all solves for theoretical foundation just drilling equation An=E, wherein, A be kernel matrix, E be each SNMR signal initial amplitude value, n that excitation pulse square is corresponding are each layer geologic body water content.Search coverage background resistivity space Distribution situation determines the characteristic of kernel matrix A, and the resistivity value of measuring point is often by some geophysical prospecting method or surrounding The existing drilling information in region obtains, and this makes kernel matrix A certainly exist certain error.

But, the inversion algorithm of contemporary literature publication is all to assume that search coverage background resistivity spatial distribution is accurate Carrying out under conditions of value, this have impact on the inversion accuracy of related algorithm, especially at background area resistivity value less than 50 Europe Under conditions of meter, when background resistivity exists bigger estimation difference, inversion result can be caused invalid.

Summary of the invention

The technical problem to be solved is to background for SNMR inversion algorithm in existing layered geology electric model The estimation difference that the distribution of zone resistance rate exists, it is provided that a kind of ground nuclear magnetic resonance inverting side based on harmonic search algorithm Method.

For solving the problems referred to above, the present invention is achieved by the following technical solutions:

A kind of ground nuclear magnetic resonance inversion method based on harmonic search algorithm, comprises the steps:

Step 1, the object function of the constraint total least square model that the just stage of drilling is set up, abstract solve for inverting Model

In formula, | | | |₂For two dimension norm；A is kernel matrix；E is that the ground nuclear-magnetism that each excitation pulse square is corresponding is total to Shake signal initial amplitude value；n_iFor each layer geologic body water content；L is the partial derivative matrix of limited model flatness；ξ is for constraint just The then regularization factors of change-total least squares journey；M_NFor water-bearing layer number in model；

Step 2, first randomly generates number equal to the solution vector of data base size HMS, and the length of each solution vector all with Water-bearing layer number M in model_NEqual；After these solution vectors are put in harmony data base HM, i.e.

In formula, HM is harmony data base；HMS is data base size；M_NFor water-bearing layer number in model；

Step 3, substitutes into each solution vector in harmony data base HM the inverting solving model of step 1 gained, asks this anti- Drill the fitness function value vector F=(F of solving model₁,F₂,…,F_HMS), and record the maximum in fitness function value vector F Value is F_max；

Step 4, the harmony data base that regeneration one is newWherein T representing matrix turn Put；This new harmony data base HM^*In each solution vectorWith data base choose probability HMCR from Sound memory storehouse HM'sIn randomly select a value；

Step 5, is adjusted probability P AR to obtained new harmony data base HM with tone^*Carry out the disturbance of bw step-length, i.e.It is updated toOrItsMore new formula be

In formula,For new harmony data base HM^*In solution vector, bw is disturbance step-length；

Step 6, by new harmony data base HM^*In each solution vector substitute into step 1 gained inverting solving model, ask The fitness function value vector of this inverting solving modelAnd record fitness function value vector F^*In Maximum be F_max ^*；

Step 7, if F_max ^*＜ F_max, then F is used_max ^*Go to update F_max, and use F_max ^*Corresponding new harmony data base HM^*Solution vector goes to update F in harmony data base HM_maxCorresponding solution vector；

Step 8, if iterations reaches maximum iteration time N set_max, or current adaptive optimal control degree functional value Less than inversion accuracy threshold value Threshold set, then stopping iteration, inversion result is minimum in fitness function value vector F Value F_minThe solution vector of the harmony data base corresponding to；Otherwise, step 4 is returned.

Described step 1 particularly as follows:

Step 1.1, sets up constraint total least square model in the stage of just drilling, and its object function is

In formula, A is kernel matrix；E is the ground nuclear magnetic resonance signal initial amplitude value that each excitation pulse square is corresponding；A^* For the kernel matrix without error；E^*For the ground nuclear magnetic resonance signal initial amplitude value without error；||·||_FFor F model Number；||·||₂For two dimension norm；n_iFor each layer geologic body water content；L is the partial derivative matrix of limited model flatness；δ is one Individual normal number；M_NFor water-bearing layer number in model；

Step 1.2, the constraint Lagrange's equation corresponding to above-mentioned object function is

In formula, A is kernel matrix；E is the ground nuclear magnetic resonance signal initial amplitude value that each excitation pulse square is corresponding；A^* For the kernel matrix without error；E^*For the ground nuclear magnetic resonance signal initial amplitude value without error；||·||_FFor F model Number；||·||₂For two dimension norm；n_iFor each layer geologic body water content；L is the partial derivative matrix of limited model flatness；δ is one Individual normal number；M_NFor water-bearing layer number in model；μ is Lagrange factor；

Step 1.3, allows δ tend to 0, and above-mentioned for step constraint Lagrange's equation is converted to inverting solving model

In formula, | | | |₂For two dimension norm；A is kernel matrix；E is that the ground nuclear-magnetism that each excitation pulse square is corresponding is total to Shake signal initial amplitude value；n_iFor each layer geologic body water content；L is the partial derivative matrix of limited model flatness；ξ is for constraint just The then regularization factors of change-total least squares journey；M_NFor water-bearing layer number in model.

In described step 1, the partial derivative matrix L of limited model flatness takes 1 time or 2 partial derivatives.

Element n in described step 2, in harmony data base HM_i,jIt it is the random number between 0 to 1.

In described step 4, data base is chosen probability HMCR and is dynamically adjusted according to iterations, i.e.

In formula, HMCR is that data base chooses probability, and Num is iterations, N_maxFor maximum iteration time.

In described step 5, tone is adjusted probability P AR dynamically to adjust according to iterations, i.e.

In formula, PAR is that tone is adjusted probability, and Num is iterations, N_maxFor maximum iteration time.

In described step 5, disturbance step-length bw dynamically adjusts according to iterations, i.e.

In formula, bw is disturbance step-length, and Num is iterations, N_maxFor maximum iteration time.

Described maximum iteration time N_maxSpan between 4800～5100；Described inversion accuracy threshold value The span of Threshold is between 0.5 × 10^-8～1.2 × 10^-8Between.

Compared with prior art, the present invention has a characteristic that

1, establish constraint R-TLS inverse model, background area resistivity Distribution Value sequence and measure signal initial Under conditions of amplitude sequence all exists error, improve inversion accuracy；

2, propose IHS improve harmonic search algorithm (the improved harmony search algorithm, It is called for short IHS) to solve problem that this model conversation is the nonlinear optimization by constraint by the derivation of equation；

3, the IHS algorithm put forward makes standard HS algorithm early stage have good ergodic, the later stage solve have high-precision Degree characteristic；

4, due to the restriction of constraints, this IHS algorithm remains to when inversion matrix equation is underdetermined equation accurately solve, Eliminate the restriction that in inverse model, water-bearing layer maximum is divided the number of plies.

Accompanying drawing explanation

Fig. 1 is E0-q correlation curve figure.

Fig. 2 is spatial resistivity scatter chart.

Fig. 3 is water content comparison figure.

Detailed description of the invention

A kind of ground nuclear magnetic resonance inversion method based on harmonic search algorithm, comprises the steps:

Step 1: the object function of the constraint total least square model set up in the just stage of drilling, abstract solves for inverting Model.

In layered medium model, by Maxwell equation group can derive point around Loop source (r, z) place, i-th layer The vertical magnetic field component H of medium_izWith radial component H_irMeet expression formula:

$H_{ir}(r,z)=\frac{I_{0}a}{2}{\∫}_0^{\mathrm{\∞}}\[a_i{e}^{-u_{i}z}-b_i{e}^{u_{i}z}\]u_i{J}_1\left(\mathrm{\λ}r\right)J_1\left(\mathrm{\λ}a\right)d\mathrm{\λ}---\left(1\right)$

$H_{iz}(r,z)=\frac{I_{0}a}{2}{\∫}_0^{\mathrm{\∞}}\[a_i{e}^{-u_{i}z}+b_i{e}^{u_{i}z}\]{\mathrm{\λJ}}_0\left(\mathrm{\λ}r\right)J_1\left(\mathrm{\λ}a\right)d\mathrm{\λ}---\left(2\right)$

Wherein, I₀For excitation current intensity； For Larmor frequency, u₀For permeability of vacuum；ρ_i It it is the i-th layer resistivity value；J₀() is the first kind 0 rank Bessel functions；J₁() is the first kind 1 rank Bessel functions；a_iAnd b_i For the parameter relevant to geological structure.

In HORIZONTAL LAYERED MEDIUM WITH HIGH ACCURACY, NMR signal initial amplitude discretization model is:

$E_0\left(q_i\right)=\underset{j=1}{\overset{\mathrm{\∞}}{\Σ}}K(q_i,z_j)n_j{\mathrm{\Δz}}_j---\left(3\right)$

Wherein, q_iFor i-th excitation pulse square, n_jWith Δ z_jIt is respectively water content and thickness, the K (q of jth layer_i,z_j) for round Kernel function after discretization under cylindrical coordinate.

Wherein, M₀For the intensity of magnetization；b_1⊥=u₀H_1⊥/I₀, H_1⊥For exciting magnetic field intensity to be perpendicular to the component in earth's magnetic field； θ_imnjFor differentiation element (r_m,Φ_n,z_j) place pull angle down；Δr_mAnd ΔΦ_nIt is respectively differential unit (r_m,Φ_n) place radial direction and Lateral length；M and N is respectively radially and the number of horizontal differentiation element.

In ground nuclear magnetic resonance detects, by launching different size of excitation pulse square q_i, to excite at different depth Hydrion in water, is converted into matrix equation by (3) formula:

An=E (5)

Wherein,I_NFor excitation pulse square number；A_ij=K (q_i,z_j)Δz_j；M_NFor water-bearing layer number in model.

In ground nuclear magnetic resonance (SNMR) inverting, kernel function precision depends on choosing of the resistivity value of each ground electric layer, And resistivity value obtains often by geophysical prospecting method detections such as vertical electrical methods, due to measurement error and the interference of various noise, Making resistivity value there is certain error, and then have influence on the precision of kernel function, actual measurement SNMR signal initial amplitude value is deposited equally In error.In order to improve SNMR inversion accuracy, the present invention establishes constraint total least square model, and its object function is:

$\begin{array}{c}\min {||(A,E)-(A^{*},E^{*})||}_F\\ \left\{\begin{array}{c}{A}^{*}n=E^{*}\\ {||Ln||}_2\≤\mathrm{\δ}\\ 0\≤n_i\≤1,i=1,2,\mathrm{...},M_N\end{array}\right.\end{array}---\left(6\right)$

Wherein, A is kernel matrix；E is the ground nuclear magnetic resonance signal initial amplitude value that each excitation pulse square is corresponding；A^* For the kernel matrix without error；E^*For the ground nuclear magnetic resonance signal initial amplitude value without error；||·||_FFor F model Number；N is each layer geologic body water content；L is the partial derivative matrix of limited model flatness, generally takes 1 time or 2 partial derivatives；δ It it is a normal number；M_NFor water-bearing layer number in model.Above-mentioned A, E are measured value, A^*、E^*For previously known theoretical value.

(6) constraint Lagrange (Lagrange) equation that formula is corresponding is:

$\{\begin{array}{c}L(A^{*},n,\mathrm{\μ})={||(A,E)-(A^{*},A^{*}n)||}_F^2+\mathrm{\μ}({||Ln||}_2^2-{\mathrm{\δ}}^2)\\ 0\≤n_i\≤1,i=1,2,\mathrm{...},M_N\end{array}---\left(7\right)$

Wherein, μ is the Lagrange factor.In order to make his result milder, it is to avoid equation is absorbed in morbid state solution, when δ tends to 0 Time, (7) formula is converted into:

$\{\begin{array}{c}\min \{\frac{{||An-E||}_2^2}{1+{||n||}_2^2}+\mathrm{\ξ}{||Ln||}_2^2\}\\ 0\≤n_i\≤1,i=1,2,\mathrm{...},M_N\end{array}---\left(8\right)$

Wherein, ξ is constraint regularization-total least square (regularization-total least square, letter Claim R-TLS) regularization factors of equation.

In order to solve the moisture content value in (8) formula object function, it is as follows that the core of IHS algorithm of the present invention realizes step:

Step 2: algorithm initialization.Randomly generating number is that (Harmony memory size is called for short data base size HMS) solution vector is put in harmony data base (Harmony memory is called for short HM), solution vector length and water-bearing layer in model Number is equal for M_N, i.e.

In formula, HM harmony data base (Harmony memory is called for short HM) is harmony data base；HMS randomly generates number It is data base size for data base size (Harmony memory size is called for short HMS)；n_i,jIt it is the random number between 0 to 1.

In a preferred embodiment of the invention, HMS=5, n_i,jIt it is the random number between 0 to 1；By each solution vector generation in HM Enter (8) formula, seek its fitness value vector Fitness=(F₁,F₂,…,F_HMS), in Fitness, maximum is F_max.In the present invention In preferred embodiment, maximum iteration time is N_max=5000；Primary iteration number of times is Num=0；Inversion accuracy threshold value is Threshold=10^-8。

Step 3: generate a new harmonyHM^*In each element With HMCR (data base chooses probability) from HM'sIn randomly select a value, be selected from probability 1-HMCR [0,1] random number between.IfTake from HM, then adjusted probability P AR it to be carried out the disturbance of bw step-length, i.e. with toneUpdate ForOrOwing to water content span is between 0 to 1, soMore new formula is

$n_i^{*}=\left\{\begin{array}{cc}0,& n_i^{*}-bw<0\\ n_i^{*}-bw,& n_i^{*}-bw\&GreaterEqual;0\\ n_i^{*}+bw,& n_i^{*}+bw\&le;1\\ 1,& n_i^{*}+bw>1\end{array}\right.---\left(9\right)$

In order to make HS algorithm early stage have good ergodic, the later stage solve and there is high accuracy characteristic, the present invention is also to mark Quasi-HS algorithm improves, and makes parameter HMCR, PAR and bw dynamically adjust according to iterations.I.e.

Step 4: calculate fitness function, updates HM.By harmony HM new in step 2^*Substitute into formula (8) and calculate fitness letter Numerical value F^*If, F^*＜ F_max, then HM is used^*Update F in HM_maxCorresponding solution vector, and use F^*Update F_max。

Step 5: judge whether iteration stops.If iterations reaches N_max, or current adaptive optimal control degree functional value is little In Threshold, then stopping iteration, inversion result is the solution vector that Fitness minima is corresponding；Otherwise, execution step is returned II。

The present invention is the ground nuclear magnetic resonance inverting of harmonic search algorithm, its propose the harmonic search algorithm of improvement in the hope of SolveMoisture content value in object function, the harmonic search algorithm of improvement makes standard and sonar surveillance system Rope algorithm early stage has good ergodic, the later stage solves and has high accuracy characteristic.Due to the restriction of constraints, this IHS calculates Method remains to when inversion matrix equation is underdetermined equation accurately solve, and eliminates in inverse model and water-bearing layer maximum is divided the number of plies Restriction.In order to improve inversion accuracy, often the thickness in water-bearing layer is set as smaller value so that the aqueous number of plies is more than exciting arteries and veins Rushing square number, inversion equation An=E becomes underdetermined equation, and in this case, the harmonic search algorithm of above-mentioned improvement is still suitable for.? In the present invention, in order to improve inversion accuracy, water-bearing layer thickness is set as 1 meter.

In the experiment of measured data inverting in the wild, inverting data come from Laos capital everything basin (Vientiane Basin, Laos) a ground water detection result of detection, choose the Site1 measuring point surveying district 1.This detection uses ground nuclear magnetic resonance With vertical electrical sounding (Vertical Electrical Sounding, be called for short VES) joint inversion, and combine probing evidence obtaining, with Verify underground hydrological geological structure.With the r.m.s. RMS evaluation inversion result performance of the moisture content value that inverting obtains, RMS formula For:

$RMS=\sqrt{\frac{1}{M_N}\underset{i=1}{\overset{M_N}{\&Sigma;}}{(n_i-n_i^{*})}^2}\&times;100\%---\left(10\right)$

Site1 measuring point region geomagnetic field intensity is 43770nT, uses the Larmor frequency of 1864Hz, and magnetic dip angle is 24 °, Laying the square coil that the length of side is 100 to detect, measured data signal to noise ratio is 6.9dB, surveys SNMR signal initial amplitude value Curve is as shown in "-the o-" in Fig. 1.VES inversion result shows, point position underground vertical depth 0～4 meters of zone resistance rate values are 550 Europe rice, vertical depth 4～11 meters of zone resistance rate values are 5000 Europe rice, and vertical depth 11～70 meters of zone resistance rate values are 20 Europe rice, electricity Resistance rate spatial distribution is as shown in Figure 2.Results of drilling shows this region vertical depth 0～38 meters of alluviation rocks formed for abandoned channel, its Main component is gravel, sand and clay；Vertical depth 38～70 meters are rock salt, and its aqueous distribution is as shown in Figure 3.In inverting, in order to Improving inversion accuracy, water-bearing layer thickness is set as 1 meter.

The inversion result of the present invention and measured data matching is preferable as can be seen from Figure 1；From figure 3, it can be seen that this method Inversion result, its r.m.s. RMS is 3.12%, and France Samovar v6.2 Inversion Software inversion result, its r.m.s. RMS It is 3.65%.It can be seen that the inversion result of the present invention is better than France's Samovar v6.2 Inversion Software, with probing knot Really the goodness of fit is higher, can reflect underground hydrological geological structure more accurately.

Claims (8)

1. a ground nuclear magnetic resonance inversion method based on harmonic search algorithm, is characterized in that comprising the steps:

Step 1, the object function of the constraint total least square model that the just stage of drilling is set up, abstract for inverting solving model

$\left\{\begin{array}{c}min\left\{\frac{\left|\right|{\mathrm{An}}_i-E|{|}_2^2}{1+\left|\right|n_i|{|}_2^2}+\mathrm{\&xi;}\left|\right|{\mathrm{Ln}}_i|{|}_2^2\right\}\\ 0\&le;n_i\&le;1,i=1,2,\mathrm{...},M_N\end{array}\right.$ ①

In formula, | | | |₂For two dimension norm；A is kernel matrix；E is the ground nuclear magnetic resonance signal that each excitation pulse square is corresponding Initial amplitude value；n_iFor each layer geologic body water content；L is the partial derivative matrix of limited model flatness；ξ for constraint regularization- The regularization factors of total least squares journey；M_NFor water-bearing layer number in model；

Step 2, first randomly generates number equal to the solution vector of data base size HMS, and the length of each solution vector is all and model Middle water-bearing layer number M_NEqual；After these solution vectors are put in harmony data base HM, i.e.

In formula, HM is harmony data base；HMS is data base size；M_NFor water-bearing layer number in model；

Step 3, substitutes into each solution vector in harmony data base HM the inverting solving model of step 1 gained, asks this inverting to ask Solve the fitness function value vector F=(F of model₁,F₂,…,F_HMS), and the maximum recorded in fitness function value vector F is F_max；

Step 4, the harmony data base that regeneration one is newThe wherein transposition of T representing matrix；Should New harmony data base HM^*In each solution vector(i=1,2 ..., M_N) with data base choose probability HMCR from and sound memory Storehouse HM'sIn randomly select a value；

Step 5, is adjusted probability P AR obtained new harmony data base HM* to be carried out the disturbance of bw step-length, i.e. with toneMore It is newlyOrItsMore new formula be

$n_i^{*}=\left\{\begin{array}{cc}0,& n_i^{*}-bw<0\\ n_i^{*}-bw,& n_i^{*}-bw\&GreaterEqual;0\\ n_i^{*}+bw,& n_i^{*}+bw\&le;1\\ 1,& n_i^{*}+bw>1\end{array}\right.$ ③

In formula,For new harmony data base HM^*In solution vector, bw is disturbance step-length；

Step 6, by new harmony data base HM^*In each solution vector substitute into step 1 gained inverting solving model, ask this anti- Drill the fitness function value vector of solving modelAnd record fitness function value vector F^*In Big value is F_max ^*；

Step 7, if F_max ^*<F_max, then F is used_max ^*Go to update F_max, and use F_max ^*Corresponding new harmony data base HM^*Xie Xiang Amount goes to update F in harmony data base HM_maxCorresponding solution vector；

Step 8, if iterations reaches maximum iteration time N set_max, or currently adaptive optimal control degree functional value is less than Inversion accuracy threshold value Threshold set, then stop iteration, and inversion result is minima F in fitness function value vector F_min The solution vector of the harmony data base corresponding to；Otherwise, step 4 is returned.

A kind of ground nuclear magnetic resonance inversion method based on harmonic search algorithm the most according to claim 1, is characterized in that Step 1 particularly as follows:

Step 1.1, sets up constraint total least square model in the stage of just drilling, and its object function is

min||(A,E)-(A^*,E^*)||_F

$\left\{\begin{array}{c}{A}^{*}{n}_i=E^{*}\\ \left|\right|L\&times;n_i|{|}_2\&le;\mathrm{\&delta;}\\ 0\&le;n_i\&le;1,i=1,2,\mathrm{...},M_N\end{array}\right.$ ④

In formula, A is kernel matrix；E is the ground nuclear magnetic resonance signal initial amplitude value that each excitation pulse square is corresponding；A^*For not Kernel matrix containing error；E^*For the ground nuclear magnetic resonance signal initial amplitude value without error；||·||_FFor F norm；| |·||₂For two dimension norm；n_iFor each layer geologic body water content；L is the partial derivative matrix of limited model flatness；δ be one just Constant；M_NFor water-bearing layer number in model；

Step 1.2, the constraint Lagrange's equation corresponding to above-mentioned object function is

$\left\{\begin{array}{c}L\left(A^{*},n_i,\mathrm{\&mu;}\right)=\left|\right|\left(A,E\right)-\left(A^{*},A^{*}{n}_i\right)|{|}_F^2+\mathrm{\&mu;}\left(\left|\right|{\mathrm{Ln}}_i|{|}_2^2-{\mathrm{\&delta;}}^2\right)\\ 0\&le;n_i\&le;1,i=1,2,\mathrm{...},M_N\end{array}\right.$ ⑤

In formula, A is kernel matrix；E is the ground nuclear magnetic resonance signal initial amplitude value that each excitation pulse square is corresponding；A^*For not Kernel matrix containing error；E^*For the ground nuclear magnetic resonance signal initial amplitude value without error；||·||_FFor F norm；| |·||₂For two dimension norm；n_iFor each layer geologic body water content；L is the partial derivative matrix of limited model flatness；δ be one just Constant；M_NFor water-bearing layer number in model；μ is Lagrange factor；

Step 1.3, allows δ tend to 0, and above-mentioned for step constraint Lagrange's equation is converted to inverting solving model

$\left\{\begin{array}{c}min\left\{\frac{\left|\right|{\mathrm{An}}_i-E|{|}_2^2}{1+\left|\right|n_i|{|}_2^2}+\mathrm{\&xi;}\left|\right|{\mathrm{Ln}}_i|{|}_2^2\right\}\\ 0\&le;n_i\&le;1,i=1,2,\mathrm{...},M_N\end{array}\right.$ ⑥

In formula, | | | |₂For two dimension norm；A is kernel matrix；E is the ground nuclear magnetic resonance signal that each excitation pulse square is corresponding Initial amplitude value；n_iFor each layer geologic body water content；L is the partial derivative matrix of limited model flatness；ξ for constraint regularization- The regularization factors of total least squares journey；M_NFor water-bearing layer number in model.

A kind of ground nuclear magnetic resonance inversion method based on harmonic search algorithm the most according to claim 2, is characterized in that In step 1, the partial derivative matrix L of limited model flatness takes 1 time or 2 partial derivatives.

A kind of ground nuclear magnetic resonance inversion method based on harmonic search algorithm the most according to claim 1, is characterized in that, In step 2, the element ni, j in harmony data base HM is the random number between 0 to 1.

A kind of ground nuclear magnetic resonance inversion method based on harmonic search algorithm the most according to claim 1, is characterized in that, In step 4, data base is chosen probability HMCR and is dynamically adjusted according to iterations, i.e.

$HMCR=05+0.47\frac{Num}{N_{\max }}$ ⑦

In formula, HMCR is that data base chooses probability, and Num is iterations, N_maxFor maximum iteration time.

A kind of ground nuclear magnetic resonance inversion method based on harmonic search algorithm the most according to claim 1, is characterized in that, In step 5, tone is adjusted probability P AR dynamically to adjust according to iterations, i.e.

$PAR=0.2+0.7(1-\frac{Num}{N_{\max }})$ ⑧

In formula, PAR is that tone is adjusted probability, and Num is iterations, N_maxFor maximum iteration time.

A kind of ground nuclear magnetic resonance inversion method based on harmonic search algorithm the most according to claim 1, is characterized in that, In step 5, disturbance step-length bw dynamically adjusts according to iterations, i.e.

$bw=0.01+0.1(1-\frac{Num}{N_{\max }})$ ⑨

In formula, bw is disturbance step-length, and Num is iterations, N_maxFor maximum iteration time.

A kind of ground nuclear magnetic resonance inversion method based on harmonic search algorithm the most according to claim 1, is characterized in that: Described maximum iteration time N_maxSpan between 4800～5100；Taking of described inversion accuracy threshold value Threshold Value scope is between 0.5 × 10^-8～1.2 × 10^-8Between.