CN108761547A - A kind of fixed-wing aviation electromagnetic data conductivity Depth Imaging method with system auxiliary parameter - Google Patents
A kind of fixed-wing aviation electromagnetic data conductivity Depth Imaging method with system auxiliary parameter Download PDFInfo
- Publication number
- CN108761547A CN108761547A CN201810916545.5A CN201810916545A CN108761547A CN 108761547 A CN108761547 A CN 108761547A CN 201810916545 A CN201810916545 A CN 201810916545A CN 108761547 A CN108761547 A CN 108761547A
- Authority
- CN
- China
- Prior art keywords
- data
- conductivity
- component
- point
- height
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V3/00—Electric or magnetic prospecting or detecting; Measuring magnetic field characteristics of the earth, e.g. declination, deviation
- G01V3/38—Processing data, e.g. for analysis, for interpretation, for correction
Landscapes
- Life Sciences & Earth Sciences (AREA)
- Engineering & Computer Science (AREA)
- Environmental & Geological Engineering (AREA)
- Geology (AREA)
- Remote Sensing (AREA)
- Physics & Mathematics (AREA)
- General Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Geophysics (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
The invention belongs to time domain aviation electromagnetic data processing method domain variabilities to disclose a kind of fixed-wing aviation electromagnetic data conductivity Depth Imaging method with system auxiliary parameter, establish each time track data table under multiple fixed altitudes, fixed swing angle, fixed yaw angle, fixed pitch angle；The data and tables of data in road of corresponding time are extracted successively；It determines position of the measurement data in table, apparent conductivity and each auxiliary parameter value is determined according to linear interpolation algorithm, if meeting setting threshold value, output conductance rate and auxiliary parameter；Judge whether that the conductivity for completing All Time road is inquired；Imaging depth is obtained by skin depth formula.Algorithm image taking speed disclosed by the invention is fast, efficient, meets the requirement of aviation transient electromagnetic processing mass data, and the fast prospecting ability to improve Mineral Resources in China establishes important science and technology basis.
Description
Technical field
The present invention relates to time domain aviation electromagnetic data processing method technical fields, and in particular to a kind of band system auxiliary ginseng
Several fixed-wing aviation electromagnetic data conductivity Depth Imaging methods.
Background technology
The detection of fixed-wing aviation electromagnetic is that one kind is adopted based on Faraday's electromagnetic induction law using fixed wing aircraft as carrier
With the geophysical prospecting method of dipole-dipole R-T unit, with exploration efficiency, the advantages such as high, at low cost, are widely used in
The fields such as geologic prospect, hydrocarbon exploration, hydrology generaI investigation.Fixed-wing aviation electromagnetic data volume is huge, if being carried out to total data anti-
Drill explanation, time and effort consuming.It is generally fast using conductivity depth imaging technique (CDI, Conductivity-Depth Imaging)
Speed processing magnanimity aviation electromagnetic data, obtain subsurface conductivity distribution map, so that it is determined that subsurface anomaly distribution situation.Traditional aviation
Electromagnetic surveying generally in the case that it is assumed that system aiding information parameter constant or known progress conductivity Depth Imaging, and
Airborne electromagnetic system can be influenced by factors such as wind speed, air pressures in practical flight measurement process, cause coil pitching, wave,
Yaw rotation and gondola swing, cause systematic parameter to change, seriously affect observation data consistency, to influence at
As effect.
Zhu Kaiguang (fixed-wing electromagnetic data double component joint conductivity Depth Imaging [J] Jilin University journal (),
2015) a kind of joint conductivity Depth Imaging algorithm tabled look-up based on the double components (Bx, Bz) in magnetic field is proposed, flying height is made
Conductivity Depth Imaging algorithm is introduced for parameter, combines to table look-up using the double components in magnetic field and determines apparent conductivity with interpolation algorithm, is kept away
Exempt from the apparent conductivity uncertain problem caused by electromagnetic data two-value, achieves preferable imaging effect；Dou Mei (bases
In fixed-wing time domain aviation electromagnetic data conductivity Depth Imaging research [D] the Jilin University of iteration look-up table, 2017)
On the basis of double component look-up tables, attitude angle pitch and flying height height is introduced into aviation electricity as parameter simultaneously
In conductance Depth Imaging, it is proposed that a kind of iteration based on pitch and height is tabled look-up CDI algorithms, and imaging is further improved
Precision；(the coil posture and gondola of fixed-wing airborne electromagnetic system are swung influences research and correction [J] geophysics to Wang Qi
Report, 2013) according to transmitting, the geometrical relationship of receiving coil relative position, swing Green tensor is acquired, and derived any attitude
Fixed-wing aviation electromagnetic in the case of angle and arbitrary swing angle responds three-component calculation expression, for band system auxiliary ginseng
Several fixed-wing aviation electromagnetic data conductivity Depth Imaging algorithms provide theoretical foundation；Application No. is 201410061917.2
Chinese invention patent：A kind of joint conductivity Depth Imaging method of the double component look-up tables of fixed-wing aviation electromagnetic data is first
Conductivity-altitude data table that temporally road divides is established, secondly by the normalization induced electromotive force response of double components
Combine and table look-up and interpolation algorithm, determine apparent conductivity, finally obtain diffusion depth with skin depth formula, exports double components joints
Conductivity Depth Imaging result.But in above-mentioned conductivity Depth Imaging algorithm, the part for having only introduced electromagnetic system flies
Row parameter (height, pitch) carries out CDI imagings, is still unavoidable from since the inconsistency of observation value influences imaging knot
The problem of fruit accuracy.Up to the present, there is not yet the fixed-wing with system auxiliary parameter (height, pitch, rall, yaw)
Aviation electromagnetic data conductivity Depth Imaging method.It quickly, simply, accurately and with system is assisted in consideration of it, how to provide one kind
The fixed-wing aviation electromagnetic data conductivity Depth Imaging method of parameter is that those skilled in the art need the technical barrier solved.
Invention content
The present invention causes systematic parameter to change for the factors such as aspect in the prior art, to seriously affect imaging essence
The problem of degree, and provide it is a kind of quickly, simple, the accurate fixed-wing aviation electromagnetic data conductivity depth with system auxiliary parameter
Imaging method.
The present invention using following technical scheme in order to solve the above technical problems, realized：
A kind of fixed-wing aviation electromagnetic data conductivity Depth Imaging method with system auxiliary parameter is designed, including as follows
Step：
Step 1：It is assumed that flying height height, swing angle roll, yaw angle yaw are fixed value,
The half space detection model that multiple conductivity Earth models are corresponded under multiple pitch angle pitch is chosen, it is one-dimensional using fixed-wing
Positive algorithm calculates B the corresponding x-component of each each model in time road, z-component electromagnetic responses, establishes what temporally road divided
Bx-Bz- σ-pitch data forms；
Step 2：According to the method described in step 1, the B that temporally road divides is established respectivelyy-Bz- σ-roll tables of data
Lattice, Bx-By- σ-yaw data forms, Bx-Bz- σ-height data forms；When transmitting coil, receiving coil posture are with respective
When attitude angle rotates, wherein calculating x-component, B electromagnetic response formula of y-component and z-component are respectively：
Wherein,
For system transmitting-receiving away from, (x, y, z) be receiving coil coordinate, J0For zero-order Bessel letter
Number, J1For first-order bessel function, λ is integration variable, ARXFor receiving coil effective area, m is transmitting magnetic moment, and h is emission lines
Circle height, z are receiving coil height, reflectance factor
Step 3：First, respectively by formula Bx(w), formula By(w), formula Bz(w) transform in the domains s, and divided by s obtain x-component,
The magnetic field transient response expression formula of y-component and z-component：
Then, respectively to formulaIt does Laplace Transform and obtains B sound of system three-component
It answers, as the tables of data source of conductivity Depth Imaging, is shown below：
Step 4:B typing aviation electromagnetic x-component, y-component and z-component electromagnetic responses, determine position of the data point in table
It sets, interpolation obtains apparent conductivity：
First, i=1 is enabled；X-component and the B number of responses strong points of z-component in i-th of time road are extracted successively；
Then, iterations k=1 is enabled；The corresponding height=height of extractionik, roll=rollik, yaw=yawikUnder
The B in road of corresponding timex-Bz- σ-pitch data forms；
Then, range formula is utilizedCalculate data point and table
The distance of middle each point, wherein Bxi、BziRespectively correspond to x-component and the B number of responses strong points of z-component in i-th of time road, Bxmk、
BzmkFor road of corresponding i-th of time, kth time iteration tables of data in any point；
Finally, it finds apart from closest approach in table, if closest approach, in the inside of tables of data, closest approach is by four quadrangles
Including the as vertex of quadrangle, calculates the corresponding linear equation in four sides of each quadrangle successively, by the x of data point points
Amount or z-component response bring into linear equation and can obtain the point corresponding z-component or x-component response on straight line, by with
Corresponding data point coordinate value compares, then can obtain the data point and the relative position of straight line, be judged successively, can enclose
The specific location at fixed number strong point；If closest approach, in the outer edge line of tables of data, closest approach is by one or two quadrangle
Including, the method for judging position is consistent with situation of the data point in table；
Step 5：After step 4 determines specific location of the data point in tables of data, tables of data partial enlargement is adopted
The apparent conductivity σ in the time road is obtained with linear interpolationik1And pitch angle pitchik, formula is as follows：At this point, the conductivity on four vertex of quadrangle is denoted as respectively The pitch angle on four vertex is denoted as p respectively1, p2,
p3, p4, p1=p4, p2=p3；Cross data point do one with point 1, the parallel straight line of 4 lines of point respectively with both sides phase above and below quadrangle
It hands over, d1It is data point at a distance from the intersection point of top, d2It is data point at a distance from following intersection point；Cross data point do one with point 1,
The parallel straight line of 2 lines of point intersects with quadrangle the right and left respectively, d3It is data point at a distance from the intersection point of the left side, d4For data
Point is at a distance from the intersection point of the right.
Step 6：Next extracting i-th of time road successively corresponds to height=heightik, pitch=pitchik,
Yaw=yawikUnder By-Bz- σ-roll data forms table look-up with reference to step 3 to step 5 and obtain the time with interpolation
The apparent conductivity σ in roadik2And swing angle rolli(k+1)；I-th of time road is extracted successively corresponds to height=heightik,
Pitch=pitchik, roll=rolli(k+1)Under Bx-By- σ-yaw data forms, are looked into reference to step 3 to step 5
Table obtains the apparent conductivity σ in the time road with interpolationik3And yaw angle yawi (k+1)；I-th of time road is extracted successively to correspond to
Pitch=pitchik, roll=rolli(k+1), yaw=yawi(k+1)Under Bx-Bz- σ-height data forms, with reference to step
Three table look-up to step 5 obtains the apparent conductivity σ in the time road with interpolationik4And flying height heighti(k+1)；
Step 7：Compare heightikWith | heighti(k+1)Size, if precision be less than L, iteration stopping, output regard electricity
ConductanceIf precision is not less than L, kth=k+1 iteration is carried out, then return
It returns to step 4 and above-mentioned steps is repeated；
Step 8：Judge whether that the conductivity inquiry for completing All Time road enables i=i+1, be back to if not completing
Above-mentioned steps are repeated in step 4；If completing to get to the corresponding apparent conductivity in electromagnetic response institute's having time road.
Step 9：Imaging depth is obtained by skin depth formula, exports the fixed-wing aviation electricity with system auxiliary parameter
The imaging results of magnetic data conductivity Depth Imaging algorithm.
Preferably, the reflection R0Y in formula1Pass through recurrence formula：
Obtained by being calculated, wherein k=n-1, n-2 ..., 1, Yn=Nn, σk,
μk, tkRespectively multi-layered earth kth layer conductivity, magnetic conductivity and thickness, ω are angular frequency.
Preferably, the unlimited integral of zero and first order Bessel function utilizes in the step 3 Laplace transform formula
Numerical value filter algorithm calculates.
Preferably, the value of the L is 1.
Preferably, the conductivity information σ to be tabled look-up using iterationiEach time is calculated using skin depth formula
The corresponding imaging depth d in roadi(1≤i≤n) is：Wherein n is time road road number, tiFor i-th of time
Road corresponding time, magnetic permeability μ in air0=4 π 10-7H/m, σiIt is corresponding for i-th of time road being tabled look-up by iteration
Conductivity.
A kind of fixed-wing aviation electromagnetic data conductivity Depth Imaging method with system auxiliary parameter proposed by the present invention,
Advantageous effect is：
(1) electromagnetic response that aviation electromagnetic detects three components is utilized in algorithm synthesis disclosed by the invention, by aircraft
Flight attitude angle (pitch, roll, yaw) and flying height (height) while introducing aviation conductivity as parameter
In Depth Imaging, establishes the iteration with height, pitch, roll, yaw and table look-up CDI algorithms；
(2) algorithm disclosed by the invention not only solves the problem of conventional method loss multi-component measurements information, simultaneously also
Avoid the imaging effect offset issue caused by state of flight inaccuracy；
(3) iteration look-up table conductivity Depth Imaging process of the algorithm disclosed by the invention with auxiliary parameter is only simple
One-Dimensional Half-Space forward modelling, matching search, iteration are tabled look-up, and linear interpolation, image taking speed is fast, efficient, meets aviation transition
The requirement of Electromagnetic Treatment mass data so that possibility is treated as in real time to aviation electromagnetic multi-component data, to improve China
The fast prospecting ability of mineral resources establishes important science and technology basis.
Description of the drawings
The present invention is described in further detail for embodiment in below in conjunction with the accompanying drawings, but does not constitute to the present invention's
Any restrictions.
Fig. 1 is the flow diagram of the method for the present invention；
Fig. 2 is the 8th corresponding Bx-Bz- σ-pitch tables of data schematic diagrames in the embodiment of the present invention；
Fig. 3 be in the embodiment of the present invention data point in Bx-BzPartial enlarged view in-σ-pitch tables of data；Wherein * is data
Point, are in table apart from data point closest approach；
Fig. 4 is the imaging results of the fixed-wing aviation electromagnetic data CDI of the embodiment of the present invention；Figure dotted line and solid line difference
Imaging results for the fixed-wing aviation electromagnetic data CDI with system auxiliary parameter and practical Earth model；
Fig. 5 is the CDI imaging results under fixed wing aircraft flight auxiliary parameter unknown situation in the embodiment of the present invention；In figure
Dotted line and solid line are respectively CDI imaging results and practical Earth model under fixed wing aircraft flight auxiliary parameter unknown situation；
Fig. 6 is measurement data attenuation curve figure in the embodiment of the present invention, and wherein Fig. 6 a are that the fields B of x-component respond；Fig. 6 b are
B responses of y-component；Fig. 6 c are B responses of z-component.
Specific implementation mode
Following will be combined with the drawings in the embodiments of the present invention, and technical solution in the embodiment of the present invention carries out clear, complete
Site preparation describes, it is clear that described embodiments are only a part of the embodiments of the present invention, instead of all the embodiments.
Shown in attached drawing 1, first, setting flying height h=120m, roll=0 ° of swing angle, yaw angle yaw=
In the case of 0 °, pitch angle pitch chooses 41 at -20 ° -20 ° at equal intervals, the conductivity of half space model 0.0025~
20 closely equal log intervals choose 64, form 41 × 64 groups of half space detection models, utilize the one-dimensional forward modeling of fixed-wing aviation electromagnetic
Algorithm, the fields B for calculating the corresponding x-component of 14 each models of time road (0.01~10ms) and z-component respond, and temporally road is drawn
Point, establish 14 Bx-BzConductivity-pitch angle data form, the 8th B is given refering to attached drawing 2x-BzConductivity-is bowed
Elevation scale is according to table, and horizontal axis is the fields the B response of x-component in table, and the longitudinal axis is that fields B of z-component responds, by each pitch angle
The point of corresponding multiple conductivity connects, and shape is into a line, and the corresponding multiple pitch angles of each conductivity are connected
Get up, also shape is into a line, one reticular structure of interlaced composition.Remaining By-Bz- σ-roll data forms, Bx-By-σ-
Yaw tables of data, Bx-Bz- σ-height tables of data and so on.When transmitting coil, receiving coil posture are with respective attitude angle
When rotating, the secondary field expression formula that three-component receiving coil receives is, wherein calculating x-component, y-component and z-component B
Electromagnetic response formula is respectively：
Wherein,
For system transmitting-receiving away from, (x, y, z) be receiving coil coordinate.J0For zero Bessel function, J1It is one
Rank Bessel function, λ are integration variable, ARXFor receiving coil effective area, m is transmitting magnetic moment, and h is transmitting coil height, and z is
Receiving coil height, reflectance factorAir permeability μ0=4 π × 10-7H/
M, Y1Pass through recurrence formula：Calculate k=n-1, n ..., 1, Yn=Nn, σk, μk, tkRespectively multi-layered earth the conductivity, magnetic conductivity and thickness
Degree, ω is angular frequency.Formula (1), formula (2), formula (3) are transformed in domain (complex frequency domain) respectively, and divided by obtain component, component and
The magnetic field transient response expression formula of component：
Laplace Transform is done to formula (4)-(6) respectively and obtains the B responses of system three-component, as conductivity Depth Imaging
Tables of data source, as shown in formula (7)-(9).
The unlimited integral (Hankel transform) of zero and first order Bessel function utilizes numerical value filter algorithm in formula
(Guptasarma120 points and 140 points) calculate；
Secondly, typing measurement data extracts B aviation electromagnetic x-component, y-component and the z-component number of responses in road of corresponding time
According to table, determines position of the data point in table, apparent conductivity is obtained by interpolation；
(1) it is 110m that one group of typing, which uses above system device, transmitting coil height, and pitch angle is 2 °, swing angle
It it is 2 °, yaw angle is 2 °, and Earth model is three layers of multi-layered earth, and conductivity is respectively σ1=0.02S/m, σ2=0.2S/m, σ3
=0.002S/m, the earth depth are respectively d1=100m, d2B the x-component of=20m, y-component and z-component responses, response are bent
Line as shown in fig. 6, illustrated by taking wherein the 8th electromagnetic data as an example,
(2) i=1 is enabled；
(3) x-component of extraction the 8th and z-component field number of responses strong point (the * points in Fig. 3)；
(4) iterations k=1 is enabled；
(5) corresponding height=120m, the B under roll=0 °, yaw=0 ° are extractedx-Bz- σ-pitch data forms；Specifically
Data are refering to shown in attached drawing 2；
(6) matching search is first with range formula；
Utilize formulaCalculate each point in data point and table away from
From wherein Bxi、BziRespectively correspond to component and the component fields number of responses strong point in the 8th time road, Bxmk、BzmkIt is the corresponding 8th
A time road, the 1st iteration tables of data in any point.Be at that time apart from minimum point (refering to the o points in attached drawing 3),
The nearest position of distance measurement data i.e. in table, closest approach is by four quadrangles in the inside of tables of data, then closest approach at this time
Including, since tables of data interval is smaller, there is almost linear feature, it is straight line to be approximately considered in table between every 2 points, according to
It is secondary to calculate the corresponding linear equation in each four sides of quadrangle.
The corresponding equation in four sides of upper left side quadrangle is followed successively by by clock-wise order in attached drawing 3：
z1=-0.441x1+2.06×10-6, z2=1.5618x2+0.5708×10-6, z3=-0.461x3+2.3209×
10-6, z4=1.504x4+0.5552×10-6；The corresponding equation in four sides of upper right side quadrangle presses clock-wise order successively
For：z1=1.504x1+0.5552×10-6, z2=-0.4792x2+2.3371×10-6, z3=1.448x3-0.3602×10-6,
z4=-0.4579x4+2.0736×10-6；；The corresponding equation in four sides of lower left quadrangle is followed successively by by clock-wise order：z1
=1.659x1+0.4598×10-6, z2=-0.441x2+2.06×10-6, z3=1.596x3+0.4839×10-6, z4=
0.4171x4+1.2693×10-6；The corresponding equation in four sides of lower right quadrangle is followed successively by by clock-wise order：z1=
1.596x1+0.4839×10-6, z2=-0.4579x2+2.0736×10-6, z3=1.5356x3-2.3056×10-7, z4=-
0.4356x4+1.7989×10-6.At this time by data point (7.4935 × 10-7, 1.6931 × 10-6) x-component or z-component response
Bring z in the corresponding equation of upper left side quadrangle into1=1.7295 × 10-6, z3=1.9754 × 10-6, x2=7.1859 × 10-7,
x4=7.5658 × 10-7,
z-z1< 0, z-z3< 0, x-x2> 0, x-x4< 0, i.e. (z-z1)(z-z3) > 0, (x-x2)(x-x4) < 0；Therefore,
Data point is not in upper left quadrangle；By data point (7.4935 × 10-7, 1.6931 × 10-6) x-component or z-component ring
Z in the corresponding equation of upper right side quadrangle should be brought into2=1.978 × 10-6, z4=1.7305 × 10-6, x1=7.5658 × 10-7,
x-x1< 0, x-x3< 0, i.e. (z-z2)(z-z4) > 0, (x-x1)(x-x3) > 0, therefore, data point is not in the quadrangle in upper right side
It is interior；By data point (7.4935 × 10-7, 1.6931 × 10-6) x-component or z-component response bring lower-left z into2=1.7295 × 10-6, z4=1.5819 × 10-6, x1=7.2007 × 10-7,
x3=7.5789 × 10-7, z-z2< 0, z-z4> 0, x-x1> 0, x-x3< 0, i.e. (z-z2)(z-z4) < 0, (x-
x1)(x-x3) < 0, therefore, data point (as shown in Fig. 3) in the quadrangle of lower left；By data point (7.4935 × 10-7,
1.6931×10-6) x-component or z-component response bring into the corresponding equation of lower right quadrangle
z2=1.7305 × 10-6, z4=1.4725 × 10-6,
x1=7.5764 × 10-7, x3=1.2527 × 10-6, z-z2< 0, z-z4> 0,
x-x1< 0, x-x3< 0, i.e. (z-z2)(z-z4) < 0,
(x-x1)(x-x3) > 0, therefore, data point is not in the quadrangle of lower right.
(7) since the interval of tables of data is smaller, there is almost linear feature, determining that data point is specific in tables of data
After position, by tables of data partial enlargement, shown in attached drawing 3, the apparent conductivity that linear interpolation obtains the time road can be used
σik1.Shown in its interpolation formula such as formula (10).At this point, the conductivity difference on four vertex (point 1, point 2, point 3, point 4) of quadrangle
For WhereinCross data point do one with point 1, the parallel straight line of 4 lines of point respectively with both sides above and below quadrangle
Intersection, d1=2.3345 × 10-8It is data point at a distance from the intersection point of top, d2=3.6425 × 10-8It is handed over following for data point
The distance of point；The pitch angle on four vertex (point 1, point 2, point 3, point 4) of quadrangle is respectively p1=1 °, p2=2 °, p3=
2 °, p4=1 °, wherein p1=p4, p2=p3.Cross data point do one with point 1, the parallel straight line of 2 lines of point respectively with a quadrangle left side
Right both sides intersection, d3=1.7771 × 10-7It is data point at a distance from the intersection point of the left side, d4=3.366 × 10-8For data point and the right side
The distance of side intersection point；
(8) next the 8th time road of extraction corresponds to height=120m, under pitch=1.8408 °, yaw=0 °
By-Bz- σ-roll data forms table look-up and obtain the apparent conductivity σ in the time road with interpolationik2It=0.0315S/m and waves
Roll=2.4961 ° of angle；It extracts the 8th time road and corresponds to height=120m, pitch=1.8408 °, roll=
2.4961 the B under °x-By- σ-yaw data forms table look-up and obtain the apparent conductivity σ in the time road with interpolationik3=
0.0352S/m and yaw angle yawi(k+1)=1.6563 °；It is pitch=1.8408 ° corresponding to extract the 8th time road, roll
B under=2.4961 °, yaw=1.6563 °x-Bz- σ-height data forms, table look-up and obtain the time road with interpolation
Apparent conductivity σik4=0.1075S/m and flying height heighti(k+1)=129.48m.
(9) compare heightik=110m and heighti(k+1)The size of=129.48m, precision are more than 1 meter, then carry out the 2nd
Above-mentioned steps are repeated back to (5) in secondary iteration, and until precision is less than 1 meter, iteration terminates, and exports apparent conductivity
(10) judge whether that the conductivity inquiry for completing All Time road is back to (2) and is repeated if not completing
State step；If completing to get to the corresponding apparent conductivity in electromagnetic response institute's having time road.
Third walks, and by skin depth formula and calculates imaging depth, exports the fixed-wing aviation with system complex parameter
Electromagnetic data conductivity Depth Imaging result.
The conductivity information σ to be tabled look-up using iterationiEach time road is calculated using skin depth formula to correspond to
Diffusion depth di(1≤i≤n) is：
Wherein magnetic permeability μ in air0=4 π 10-7H/m, σiWhen=0.1075S/m is i-th to be tabled look-up by iteration
Between the corresponding conductivity in road, to CDI of the final output with system auxiliary parameter as a result, result refering to shown in attached drawing 4；Attached drawing 5
Give CDI imaging contrast's results under fixed wing aircraft flight auxiliary parameter unknown situation.
The foregoing is only a preferred embodiment of the present invention, but scope of protection of the present invention is not limited thereto,
Any one skilled in the art in the technical scope disclosed by the present invention, according to the technique and scheme of the present invention and its
Inventive concept is subject to equivalent substitution or change, should be covered by the protection scope of the present invention.
Claims (5)
1. a kind of fixed-wing aviation electromagnetic data conductivity Depth Imaging method with system auxiliary parameter, which is characterized in that packet
Include following steps：
Step 1：It is assumed that flying height height, swing angle roll, yaw angle yaw are fixed value, choose
The half space detection model that multiple conductivity Earth models are corresponded under multiple pitch angle pitch utilizes the one-dimensional forward modeling of fixed-wing
Algorithm calculates B the corresponding x-component of each each model in time road, z-component electromagnetic responses, establishes the B that temporally road dividesx-Bz-
σ-pitch data forms；
Step 2：According to the method described in step 1, the B that temporally road divides is established respectivelyy-Bz- σ-roll data forms, Bx-
By- σ-yaw data forms, Bx-Bz- σ-height data forms；When transmitting coil, receiving coil posture are with respective attitude angle
When degree rotates, wherein calculating x-component, B electromagnetic response formula of y-component and z-component are respectively：
Wherein,
For system transmitting-receiving away from, (x, y, z) be receiving coil coordinate, J0For zero Bessel function, J1It is one
Rank Bessel function, λ are integration variable, ARXFor receiving coil effective area, m is transmitting magnetic moment, and h is transmitting coil height, and z is
Receiving coil height, reflectance factor
Step 3：First, respectively by formula Bx(w), formula By(w), formula Bz(w) it transforms in the domains s, and divided by s obtains x-component, y divides
The magnetic field transient response expression formula of amount and z-component：
Then, respectively to formulaIt does Laplace Transform and obtains the B responses of system three-component, make
For the tables of data source of conductivity Depth Imaging, it is shown below：
Step 4:B typing aviation electromagnetic x-component, y-component and z-component electromagnetic responses, determine position of the data point in table,
Interpolation obtains apparent conductivity：
First, i=1 is enabled；X-component and the B number of responses strong points of z-component in i-th of time road are extracted successively；
Then, iterations k=1 is enabled；The corresponding height=height of extractionik, roll=rollik, yaw=yawikLower correspondence
The B in time roadx-Bz- σ-pitch data forms；
Then, range formula is utilizedIt calculates each in data point and table
The distance of point, wherein Bxi、BziRespectively correspond to x-component and the B number of responses strong points of z-component in i-th of time road, Bxmk、Bzmk
For road of corresponding i-th of time, kth time iteration tables of data in any point；
Finally, it finds apart from closest approach in table, if closest approach, in the inside of tables of data, closest approach is wrapped by four quadrangles
Contain, as the vertex of quadrangle, the corresponding linear equation in four sides of each quadrangle is calculated successively, by the x-component or z of data point
Component response, which is brought into linear equation, can obtain the point corresponding z-component or x-component response on straight line, by with it is corresponding
Data point coordinate value compares, then can obtain the data point and the relative position of straight line, be judged successively, can draw a circle to approve data
The specific location of point；If closest approach, in the outer edge line of tables of data, closest approach is to be included by one or two quadrangle
, the method for judging position is consistent with situation of the data point in table；
Step 5：After step 4 determines specific location of the data point in tables of data, by tables of data partial enlargement, using line
Property interpolation obtains the apparent conductivity σ in the time roadik1And pitch angle pitchik, formula is as follows：At this point, the conductivity on four vertex of quadrangle is denoted as respectively The pitch angle on four vertex is denoted as p respectively1, p2,
p3, p4, p1=p4, p2=p3；Cross data point do one with point 1, the parallel straight line of 4 lines of point respectively with both sides phase above and below quadrangle
It hands over, d1It is data point at a distance from the intersection point of top, d2It is data point at a distance from following intersection point；Cross data point do one with point 1,
The parallel straight line of 2 lines of point intersects with quadrangle the right and left respectively, d3It is data point at a distance from the intersection point of the left side, d4For data
Point is at a distance from the intersection point of the right.
Step 6：Next extracting i-th of time road successively corresponds to height=heightik, pitch=pitchik, yaw=
yawikUnder By-Bz- σ-roll data forms table look-up with reference to step 3 to step 5 and obtain regarding for the time road with interpolation
Conductivityσik2And swing angle rolli(k+1)；I-th of time road is extracted successively corresponds to height=heightik, pitch=
pitchik, roll=rolli(k+1)Under Bx-By- σ-yaw data forms, are tabled look-up with reference to step 3 to step 5 and interpolation
Obtain the apparent conductivity σ in the time roadik3And yaw angle yawi(k+1)；I-th of time road is extracted successively corresponds to pitch=
pitchik, roll=rolli(k+1), yaw=yawi(k+1)Under Bx-Bz- σ-height data forms, with reference to step 3 to step
Five, which table look-up, obtains the apparent conductivity σ in the time road with interpolationik4And flying height heighti(k+1)；
Step 7：Compare heightikAnd heighti(k+1)Size, if precision be less than L, iteration stopping, export apparent conductivityIf precision is not less than L, kth=k+1 iteration is carried out, is again returned to
Above-mentioned steps are repeated in step 4；
Step 8：Judge whether that the conductivity inquiry for completing All Time road enables i=i+1, be back to step if not completing
Four are repeated above-mentioned steps；If completing to get to the corresponding apparent conductivity in electromagnetic response institute's having time road.
Step 9：Imaging depth is obtained by skin depth formula, exports the fixed-wing aviation electromagnetic number with system auxiliary parameter
According to the imaging results of conductivity Depth Imaging algorithm.
2. a kind of fixed-wing aviation electromagnetic data conductivity Depth Imaging with system auxiliary parameter according to claim 1
Method, which is characterized in that the reflection R0Y in formula1Pass through recurrence formula：
Obtained by being calculated, wherein k=n-1, n-2 ..., 1, Yn=Nn,
σk, μk, tkRespectively multi-layered earth kth layer conductivity, magnetic conductivity and thickness, ω are angular frequency.
3. a kind of fixed-wing aviation electromagnetic data conductivity Depth Imaging with system auxiliary parameter according to claim 1
Method, which is characterized in that the unlimited integral profit of zero and first order Bessel function in the step 3 Laplace transform formula
It is calculated with numerical value filter algorithm.
4. a kind of fixed-wing aviation electromagnetic data conductivity Depth Imaging with system auxiliary parameter according to claim 1
Method, which is characterized in that the value of the L is 1.
5. a kind of fixed-wing aviation electromagnetic data conductivity Depth Imaging with system auxiliary parameter according to claim 1
Method, which is characterized in that the conductivity information σ to be tabled look-up using iterationiPer a period of time is calculated using skin depth formula
Between the corresponding imaging depth d in roadi(1≤i≤n) is：Wherein n is time road road number, tiWhen being i-th
Between the road corresponding time, magnetic permeability μ in air0=4 π 10-7H/m, σiI-th of time road to be tabled look-up by iteration corresponds to
Conductivity.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810916545.5A CN108761547A (en) | 2018-08-13 | 2018-08-13 | A kind of fixed-wing aviation electromagnetic data conductivity Depth Imaging method with system auxiliary parameter |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810916545.5A CN108761547A (en) | 2018-08-13 | 2018-08-13 | A kind of fixed-wing aviation electromagnetic data conductivity Depth Imaging method with system auxiliary parameter |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN108761547A true CN108761547A (en) | 2018-11-06 |
Family
ID=63969835
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201810916545.5A Pending CN108761547A (en) | 2018-08-13 | 2018-08-13 | A kind of fixed-wing aviation electromagnetic data conductivity Depth Imaging method with system auxiliary parameter |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN108761547A (en) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110244367A (en) * | 2019-06-17 | 2019-09-17 | 吉林大学 | A kind of ZTEM posture compensation method based on the more base stations in ground |
Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US5471056A (en) * | 1992-09-25 | 1995-11-28 | Texaco Inc. | Airborne scanner image spectrometer |
| CN101710187A (en) * | 2009-12-17 | 2010-05-19 | 成都理工大学 | Method for calibrating time domain aviation electromagnetic altitude |
| CN103675927A (en) * | 2013-12-20 | 2014-03-26 | 吉林大学 | Correction method for pendulum angle of receiving pod of airborne electromagnetic system in fixed wing aircraft |
| CN103675926A (en) * | 2012-09-24 | 2014-03-26 | 成都理工大学 | Conductivity-depth conversion method for aviation transient electromagnetic data |
| CN103792587A (en) * | 2014-02-24 | 2014-05-14 | 吉林大学 | Fixed wing aviation electromagnetic data double-component look-up table method combination conductivity depth imaging method |
| US20140312905A1 (en) * | 2013-04-22 | 2014-10-23 | Brent D. Wheelock | Reverse Semi-Airborne Electromagnetic Prospecting |
| CN108120439A (en) * | 2017-12-21 | 2018-06-05 | 北华航天工业学院 | A kind of three-component induction coil attitude measurement method and device |
-
2018
- 2018-08-13 CN CN201810916545.5A patent/CN108761547A/en active Pending
Patent Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US5471056A (en) * | 1992-09-25 | 1995-11-28 | Texaco Inc. | Airborne scanner image spectrometer |
| CN101710187A (en) * | 2009-12-17 | 2010-05-19 | 成都理工大学 | Method for calibrating time domain aviation electromagnetic altitude |
| CN103675926A (en) * | 2012-09-24 | 2014-03-26 | 成都理工大学 | Conductivity-depth conversion method for aviation transient electromagnetic data |
| US20140312905A1 (en) * | 2013-04-22 | 2014-10-23 | Brent D. Wheelock | Reverse Semi-Airborne Electromagnetic Prospecting |
| CN103675927A (en) * | 2013-12-20 | 2014-03-26 | 吉林大学 | Correction method for pendulum angle of receiving pod of airborne electromagnetic system in fixed wing aircraft |
| CN103792587A (en) * | 2014-02-24 | 2014-05-14 | 吉林大学 | Fixed wing aviation electromagnetic data double-component look-up table method combination conductivity depth imaging method |
| CN108120439A (en) * | 2017-12-21 | 2018-06-05 | 北华航天工业学院 | A kind of three-component induction coil attitude measurement method and device |
Non-Patent Citations (4)
| Title |
|---|
| 李冰冰: ""基于多分量测量的固定翼航空电磁数据电导率深度成像研究"", 《中国优秀硕士学位论文全文数据库 基础科学辑》 * |
| 王琦: ""固定翼时间域航空电磁数据整体反演"", 《中国优秀博士学位论文全文数据库 基础科学辑》 * |
| 窦梅: ""基于迭代查表法的固定翼时间域航空电磁数据电导率深度成像"", 《中国优秀硕士学位论文全文数据库 基础科学辑》 * |
| 西永在 等: ""固定翼时间域航空电磁系统的线圈姿态影响研究"", 《工程地球物理学报》 * |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110244367A (en) * | 2019-06-17 | 2019-09-17 | 吉林大学 | A kind of ZTEM posture compensation method based on the more base stations in ground |
| CN110244367B (en) * | 2019-06-17 | 2020-05-29 | 吉林大学 | Attitude compensation method of ZTEM system based on multiple ground base stations |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Reid et al. | Airborne electromagnetic footprints in 1D earths | |
| CN108254780A (en) | A kind of microseism positioning and anisotropic velocity structure tomographic imaging method | |
| CN108984818A (en) | Fixed-wing time domain aviation electromagnetic data intend restricted by three-dimensional space entirety inversion method | |
| CN105279330B (en) | The numerical value emulation method of sea moving ship turbulent wake | |
| CN107992676A (en) | A kind of high-speed simulation modeling method of moving target time domain scatter echo | |
| CN107315173A (en) | A kind of GPR and differential GPS method for synchronizing time and system | |
| CN110146846A (en) | A kind of sound source position estimation method, readable storage medium storing program for executing and computer equipment | |
| CN106842080B (en) | A kind of magnetic field measuring device posture swing interference minimizing technology | |
| CN103513235B (en) | Clear sky aircraft wake stable section radar scattering characteristic computing method | |
| CN105136073B (en) | A kind of meteorological calibration model in deformation of slope monitoring | |
| CN106168682B (en) | A kind of moving target body monitoring method based on rotational gravity field | |
| CN106443803B (en) | Ocean controllable source electromagnetic response computational methods based on emitter actual measurement morphological data | |
| CN106526596A (en) | Data processing method of synthetic aperture radar ocean wind-field retrieval variation model | |
| CN107607936A (en) | A kind of high frequency day earthwave Radar Sea ocean surface flow inversion method | |
| Zuo et al. | Surface turbulent flux measurements over the Loess Plateau for a semi-arid climate change study | |
| CN109975880A (en) | A kind of orientation method based on characteristic vector, apparatus and system | |
| CN110399680A (en) | A kind of shallow sea elastic construction radiated sound field calculation method | |
| CN108761547A (en) | A kind of fixed-wing aviation electromagnetic data conductivity Depth Imaging method with system auxiliary parameter | |
| CN107607998A (en) | A kind of nuclear magnetic resonance water detector Magnetic Resonance parameter extracting method and system | |
| CN109581492A (en) | Petrophysical parameter calculation method and system based on Simulating Seismic Wave | |
| Liou et al. | On the initialization of tropical cyclones with a three dimensional variational analysis | |
| CN105093300B (en) | A kind of boundary recognition of geological body method and device | |
| CN109490978A (en) | A kind of frequency domain quick high accuracy forward modeling method on fluctuating stratum | |
| CN108549102A (en) | The earth formation Curvature Estimation method of the more window analyses of joint gradient-structure tensor sum | |
| CN109143362A (en) | Scattered wave separation method based on total scattering angle gathers |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| RJ01 | Rejection of invention patent application after publication |
Application publication date: 20181106 |
|
| RJ01 | Rejection of invention patent application after publication |