WO2021068527A1 - 一种足迹引导的高效航空电磁法数值模拟方法 - Google Patents
一种足迹引导的高效航空电磁法数值模拟方法 Download PDFInfo
- Publication number
- WO2021068527A1 WO2021068527A1 PCT/CN2020/093495 CN2020093495W WO2021068527A1 WO 2021068527 A1 WO2021068527 A1 WO 2021068527A1 CN 2020093495 W CN2020093495 W CN 2020093495W WO 2021068527 A1 WO2021068527 A1 WO 2021068527A1
- Authority
- WO
- WIPO (PCT)
- Prior art keywords
- footprint
- area
- station
- calculation
- green
- Prior art date
Links
- 238000000034 method Methods 0.000 title claims abstract description 146
- 238000004088 simulation Methods 0.000 title claims abstract description 21
- 238000004364 calculation method Methods 0.000 claims abstract description 110
- 230000005672 electromagnetic field Effects 0.000 claims abstract description 25
- 230000004044 response Effects 0.000 claims abstract description 9
- 230000014509 gene expression Effects 0.000 claims description 57
- 230000005684 electric field Effects 0.000 claims description 42
- 239000011159 matrix material Substances 0.000 claims description 28
- 230000005284 excitation Effects 0.000 claims description 16
- 230000002159 abnormal effect Effects 0.000 claims description 9
- 230000008569 process Effects 0.000 claims description 8
- 238000000638 solvent extraction Methods 0.000 abstract 2
- 230000006870 function Effects 0.000 description 19
- 238000004613 tight binding model Methods 0.000 description 15
- 238000010586 diagram Methods 0.000 description 12
- 238000005259 measurement Methods 0.000 description 12
- 238000001514 detection method Methods 0.000 description 10
- 238000009826 distribution Methods 0.000 description 7
- 230000002547 anomalous effect Effects 0.000 description 2
- 238000009434 installation Methods 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 230000035699 permeability Effects 0.000 description 2
- 238000013459 approach Methods 0.000 description 1
- 230000008901 benefit Effects 0.000 description 1
- 239000001354 calcium citrate Substances 0.000 description 1
- 238000004422 calculation algorithm Methods 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 238000002224 dissection Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000011156 evaluation Methods 0.000 description 1
- 238000012854 evaluation process Methods 0.000 description 1
- 239000007787 solid Substances 0.000 description 1
- 238000009827 uniform distribution Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
Definitions
- the invention relates to the technical field of aviation electromagnetic forward modeling, in particular to a footprint-guided high-efficiency aviation electromagnetic method numerical simulation method.
- the aerial electromagnetic observation device uses an airborne platform, which is characterized by high efficiency, and the detection area often covers tens of square kilometers to thousands of square kilometers.
- Traditional numerical simulation methods require that the entire detection area be used as the calculation area. This large-area grid calculation will consume a lot of calculation memory and time.
- the truncated boundary vector finite element method only requires the abnormal body in the detection area as the calculation area, which has the characteristics of high efficiency. However, in the face of the existence of a large range of abnormal bodies, a large number of Green's function calculations are required, which makes the efficiency of this method questionable.
- the present invention provides a footprint-guided high-efficiency aviation electromagnetic method numerical simulation method, which combines the aviation electromagnetic method observation device with a limited footprint area (that is, the aviation electromagnetic method footprint Features) and the truncated boundary vector finite element method, the entire aerial electromagnetic detection area is subdivided into multiple sub-areas according to the station to calculate the electromagnetic field of the station measurement point one by one, which has the characteristics of proper calculation area and high calculation efficiency.
- the specific technical solutions are as follows:
- the present invention provides a footprint-guided high-efficiency aviation electromagnetic method numerical simulation method, which includes the following steps:
- Step S100 Use a uniform grid to divide the aerial electromagnetic method to detect the area
- Step S200 Calculate the extension scales of the footprint area of the aviation electromagnetic observation device used in the x, y, and z directions; the number of grids in the x and y directions of the footprint area is equal, and the number of grids in the z direction is x or Between 1/4 and 1/2 of the number of grids in the y direction; the specific number of grids is determined by the ratio between the analytical solution and the numerical solution of the electromagnetic field at the measurement point of the uniform half-space model;
- Step S300 Use the truncated boundary vector finite element method to take the footprint corresponding to the first station as an abnormal body, calculate the scattered current in the footprint, and store the Green function generated during the calculation process for future use;
- the Green function is generated when the truncated boundary vector finite element is used, and is used to calculate the electromagnetic field generated by the boundary of the truncated boundary vector finite element method in the calculation of the scattered current of the discrete elements in the footprint;
- Step S400 Calculate the electromagnetic field response of the measuring point by using the scattered current of the discrete unit in the footprint of the first measuring station, and store the Green function generated in the calculation process for backup;
- the Green's function is used to calculate the electromagnetic field generated by the scattered current of the discrete unit in the footprint at the observation point of the station;
- Step S500 Repeat the Green's function stored in steps S300 and S400 at the first station for subsequent stations, first calculate the scattered current in the station footprint according to S300, and then calculate the electromagnetic field of the station footprint at the measurement point according to S400.
- Green's function recorded in steps S300 and S400 of the present invention it can be seen that they essentially belong to two Green's functions with different purposes. Therefore, in order to better distinguish, the Green's function recorded in step S300 is named the first Green's function.
- the Green function recorded in step S400 is named the second Green function.
- the grid in the step S100 is a regular hexahedral grid.
- step S200 specifically requires that the footprint area is centered on the projection of the emission source on the ground, and the number of grids in the x, y, and z directions is gradually increased at a ratio of 4:4:1 until the scattered current in the footprint is measured
- the error between the electromagnetic field generated by the point and the analytical solution of the measurement point is less than 5%.
- step S200 is specifically:
- Step S210 Calculate the analytical solution of the secondary electromagnetic field generated by the aviation electromagnetic emission source at the receiving point;
- the area A is an area of 400 ⁇ 400 ⁇ 100m 3 ;
- Step S220 the relative error between the calculated numerical solution of the secondary magnetic field and the obtained analytical solution of the secondary magnetic field
- Step S230 Determine the magnitude of the relative error. If the relative error between the real and imaginary parts of the magnetic field does not exceed 5%, the volume of area A in Step S210 is defined as the footprint of the aviation electromagnetic device;
- the area A in S210 will be redefined according to the 4:4:1 ratio expansion and the excitation current in the area and the secondary excitation current generated by the aeronautical electromagnetic measuring device in the area will be calculated. Field value solution, return to step S220.
- the use of the truncated boundary vector finite element method in the step S300 includes:
- Step S310 The calculation area is defined as the footprint and an element thickness wrapper, the electric field is defined at the center point of the edge of the division element of the calculation area, and the vector finite element theory is used to establish a set of equations to be solved that the electric field satisfies;
- Step S320 Calculate the Green function to connect the central current of the unit in the footprint with the boundary electric field of the calculation area, arrange the Green function in a matrix form, the number of rows is the number of edges of the calculation area, and the number of columns is the number of edges of the elements in the footprint area;
- Step S330 Store the coefficient matrix of the relational equation group in step S320 for use;
- Step S340 forming a relational equation group for the electric field at the edge center of the calculation area boundary and the electric field at the edge center of the internal element of the footprint;
- Step S350 Replace the calculated area boundary electric field in step S310 with the expression in step S320 to form a truncated boundary vector finite element method control equation set;
- Step S360 Solve the equations to obtain the electric field value at the center of the edge of the element in the footprint, and use linear interpolation to obtain the scattered current at the center of the element.
- step S400 specifically includes the following steps:
- Step S410 Calculate and connect the measuring point of the first measuring station with the scattered current Green's function of the discrete unit in the footprint of the first measuring station, and arrange it into a row matrix.
- the number of matrix columns is the number of scattered currents, that is, the number of discrete units of the footprint;
- Step S420 store the row matrix formed in step S410 for future use
- Step S430 Multiply the row matrix of Green's function by the column matrix formed by the scattered current to obtain the magnetic field value of the measuring point.
- the step S500 includes the following steps:
- Step S510 Repeat step S310, call the Green function stored in step S330, repeat steps S340-S360, and obtain the scattered current of the unit in the footprint of the subsequent station;
- Step S520 Call the Green function stored in step S420, repeat step S430, and obtain the magnetic field value of the subsequent measuring station.
- the footprint-guided high-efficiency aviation electromagnetic method numerical simulation method provided by the present invention abandons the traditional numerical simulation method that requires all detection areas as calculation areas, and chooses to calculate the detection areas one by one according to the footprint of a single station, which greatly reduces the calculation area.
- the footprint-guided high-efficiency aviation electromagnetic method numerical simulation method provided by the present invention uses uniform grid division, and all station footprints are composed of the same grid.
- the Green function stored in the first station can be reused, which greatly improves the calculation efficiency of subsequent stations.
- Fig. 1 is a schematic diagram of the numerical simulation method of the high-efficiency aviation electromagnetic method guided by the footprint in this embodiment
- Figure 2 is a schematic diagram of the form of the station device and the uniform dissection of the uniform underground medium below (using a uniform grid to divide the detection area of the aerial electromagnetic method);
- Figure 3a is the layout of the real part of the current in the underground homogeneous medium in the form of a coaxial device
- Figure 3b is the layout of the imaginary part of the current in the underground homogeneous medium in the form of a coaxial device
- Figure 4a is the layout of the real part of the current in the underground homogeneous medium in the form of a coplanar device
- Figure 4b is the layout of the imaginary part of the current in the underground homogeneous medium in the form of a coplanar device
- Figure 5 is a diagram of the evaluation process of the underground footprint area scale (the solid line box is the 400 ⁇ 400 ⁇ 100m 3 calculation area used in the first evaluation, and the dashed box is the subsequent calculation area after expansion according to 4:4:1);
- Figure 6 is a schematic diagram of the calculation area required by different methods to calculate the large horizontally extending thin plate abnormal body (Box 1 is the calculation area required by the traditional differential method, Box 2 is the calculation area required by the traditional truncated boundary vector finite element method, and Box 3 is the calculation area required by the traditional truncated boundary vector finite element method. Footprint guided truncated boundary vector finite element method required calculation area);
- Figure 7 is a schematic diagram of the footprint area of the subsequent station (the calculation area of the first station is shown by the solid black box, and the calculation area of the second station is shown by the dashed black box);
- Figure 8 is a schematic diagram of a zigzag complex three-dimensional model (Box No. 1 is the calculation area of the first station in this embodiment, Box No. 2 is the calculation area of the station No. 45 in this embodiment, and Box No. 3 is the traditional truncated boundary vector finite element method. The calculation area required by the station);
- Figure 9a is a comparison diagram of the imaginary part of the magnetic field numerical solution of the zigzag complex three-dimensional model using the method of this embodiment and the traditional truncated boundary vector finite element method when using a coaxial device (the black dot is the imaginary part of the calculation result of the method of this embodiment, and the black line is the traditional The imaginary part of the calculation result of the truncated boundary vector finite element method);
- Fig. 9b is a comparison diagram of the numerical solution real part of the magnetic field of the zigzag complex three-dimensional model between the method of this embodiment and the traditional truncated boundary vector finite element method when a coaxial device is used (the circle is the real part of the calculation result of the method of this embodiment, and the dashed line is the traditional truncated boundary Real part of vector finite element method calculation results);
- Figure 10a is a comparison diagram of the imaginary part of the magnetic field numerical solution of the complex three-dimensional model of the sawtooth between the method of this embodiment and the traditional truncated boundary vector finite element method when the coplanar device is used (the black dot is the imaginary part of the calculation result of the method of this embodiment, and the black line is the traditional The imaginary part of the calculation result of the truncated boundary vector finite element method);
- Figure 10b is a comparison diagram of the numerical solution real part of the magnetic field of the zigzag complex three-dimensional model between the method of this embodiment and the traditional truncated boundary vector finite element method when a coplanar device is used (the circle is the real part of the calculation result of the method of this embodiment, and the dashed line is the traditional truncated boundary The real part of the calculation result of the vector finite element method).
- a footprint-guided high-efficiency aviation electromagnetic method numerical simulation method includes the following steps:
- Step S100 Use a uniform grid to subdivide the airborne electromagnetic method to detect the area.
- the uniform grid used here is a regular hexahedral grid with the same scale in the x, y and z directions, and the size is equal to 5, 10 and Choose between 25 three numbers, the smaller the scale, the higher the accuracy.
- Step S200 Calculate the extension scale of the footprint area of the aviation electromagnetic observation device used in the x, y, and z directions, which specifically includes the following steps:
- the first step Calculate the uniform half-space analytical solution and the full-space analytical solution generated by the emission source of a single station at the measurement point, and use the half-space analytical solution to subtract the full-space analytical solution to obtain the secondary current generated by the uniform underground medium at the measurement point Analytical solution of magnetic field;
- Step 2 Calculate the excitation current generated by the emission source of the aviation electromagnetic device at the center of the uniform grid unit in the uniform underground A area (400 ⁇ 400 ⁇ 100m 3 area) directly below it;
- the third step Calculate the secondary field generated by the excitation current at the measuring device in the uniform grid unit of the defined volume area (400 ⁇ 400 ⁇ 100m 3 area) to superimpose, and obtain the numerical solution of the secondary magnetic field at the measuring point;
- Step 4 Calculate the relative error between the numerical solution of the secondary magnetic field obtained in the third step and the analytical solution of the secondary magnetic field obtained in the first step;
- Step 5 Determine the relative error size. If the relative error between the real and imaginary parts of the magnetic field does not exceed 5%, the volume of area A in the second step of the step is defined as the footprint of the aviation electromagnetic device;
- Step S300 Use the truncated boundary vector finite element method to take the footprint corresponding to the first station as the target area, calculate the scattered current in the footprint, and store the Green function generated during the calculation process for use;
- Step S310 The calculation area is defined as the footprint and an element thickness wrapper, the electric field is defined at the center point of the edge of the division element of the calculation area, and the vector finite element theory is used to establish a set of equations to be solved that the electric field satisfies;
- Step S320 Calculate the Green function to connect the central current of the unit in the footprint with the boundary electric field of the calculation area, arrange the Green function in a matrix form, the number of rows is the number of edges of the calculation area, and the number of columns is the number of edges of the elements in the footprint area;
- Step S330 Store the coefficient matrix of the relational equation group in step S320 for use;
- Step S340 forming a relational equation group for the electric field at the edge center of the calculation area boundary and the electric field at the edge center of the internal element of the footprint;
- Step S350 Replace the calculated area boundary electric field in step S310 with the expression in step S320 to form a truncated boundary vector finite element method control equation set;
- Step S360 Solve the equations to obtain the electric field value at the center of the edge of the element in the footprint, and use linear interpolation to obtain the scattered current at the center of the element.
- Step S400 Calculate the electromagnetic field response of the measuring point by using the scattered current of the discrete elements in the footprint of the measuring station. The steps are as follows:
- Step S410 Calculate and connect the measuring point of the first measuring station with the scattered current Green's function of the discrete unit in the footprint of the first measuring station, and arrange it into a row matrix.
- the number of matrix columns is the number of scattered currents, that is, the number of discrete units of the footprint;
- Step S420 store the row matrix formed in step S410 for future use
- Step S430 Multiply the row matrix of Green's function by the column matrix formed by the scattered current to obtain the magnetic field value of the measuring point.
- Step S500 Repeat the Green's function stored in steps S300 and S400 at the first station for subsequent stations, first calculate the scattered current in the station footprint according to S300, and then calculate the electromagnetic field of the station footprint at the measurement point according to S400.
- the specific S500 includes the following steps:
- Step S510 Repeat step S310, call the Green function stored in step S330, repeat steps S340-S360, and obtain the scattered current of the unit in the footprint of the subsequent station;
- Step S520 Call the Green function stored in step S420, repeat step S430, and obtain the magnetic field value of the subsequent measuring station.
- the aviation electromagnetic method uses an airborne platform to carry an observation device, which is composed of a transmitting coil and a receiving coil.
- an observation device which is composed of a transmitting coil and a receiving coil.
- the transmitting coil and the receiving coil are horizontal coils, it is called a coplanar device form, and when the transmitting coil and the receiving coil are vertical coils, it is called a coaxial device form.
- the altitude of the airborne electromagnetic method is between 40m and the route from right to left.
- the source of the transmitting coil is located 7.5m in front of the receiving coil, and the transmitting frequency of the transmitting coil is 10000hz.
- this embodiment uses a uniform grid underground medium.
- the vertical coil (coaxial device) used in the aerial electromagnetic transmitter can be treated as a horizontal magnetic dipole source, and the horizontal coil used in the transmitter (coplanar device) can be treated as a vertical magnetic dipole source.
- the time-harmonic factor e iwt set the rectangular coordinates, the ground surface is the xy plane, and the z coordinate is downward.
- the emission source in the air is a unit horizontal magnetic dipole source (coaxial device).
- the electric field expression generated at the receiving point (the center of the underground grid unit) satisfies the expression 1)-3):
- i is the unit imaginary number
- w is the angular velocity corresponding to the frequency of the emission source
- ⁇ x and ⁇ y are the difference between the x and y coordinates of the receiving point (the center point of the underground grid) and the aviation electromagnetic emission source, respectively
- r and ⁇ are the distance between the receiving point (the center point of the underground grid) and the dipole source and the horizontal projection distance
- ⁇ is the spatial wave number
- J 0 and J 1 are the first type 0-order and 1st-order Bessels, respectively Er function
- ⁇ is the air permeability
- ⁇ is the conductivity of the underground homogeneous medium.
- the current distribution generated by the uniform underground emission source is as follows:
- Figure 3a, Figure 3b, Figure 4a and Figure 4b respectively show the horizontal slice diagram and the vertical slice diagram of the current amplitude distribution of the horizontal magnetic dipole source (coaxial device) and the vertical magnetic dipole source (coplanar device) .
- the second column of Fig. 3b and Fig. 4b shows the current distribution of the vertical section in the direction of the source horizontal dipole. The current is greater than 10 -9.5 A/m mainly in the depth of 100m underground.
- Figure 3a, Figure 3b, Figure 4a and Figure 4b show that the underground current caused by the aviation electromagnetic emission source is mainly concentrated in the 400 ⁇ 400 ⁇ 100m 3 area directly below the emission source. In the process of footprint assessment, this area is first selected for calculation Numerical solution of the secondary magnetic field at the measuring point.
- the uniform half-space model is selected, and the magnetic field generated by the underground medium at the measuring point is compared by comparing two different methods to define the size of the footprint area of the aviation electromagnetic method device.
- the specifics are: First, the calculation of the emission source caused by the selected underground area After the current distribution; then, calculate the secondary magnetic field superimposition caused by the current in the selected area grid at the aviation electromagnetic observation point to obtain the numerical solution of the secondary magnetic field at the measurement point; finally, the obtained secondary magnetic field numerical solution is compared with the analytical solution to determine Whether the selected current area has reached the required scale.
- the judgment criterion is that the secondary magnetic field analytical solution is the standard, and the numerical solution error of the secondary magnetic field is less than 5%.
- the analytical solution is obtained by subtracting the full-space analytical solution from the uniform half-space analytical solution of the magnetic field at the measuring point.
- the internal current of the grid cell is regarded as a uniform distribution, and when calculating the magnetic field generated at the measuring point, the grid cell is regarded as an electric dipole source.
- the magnetic field generated by the unit electric dipole source in the whole space is the expression 7)-9):
- the magnetic field generated by the unit electric dipole source in the whole space is the expression 10)-12):
- i is the unit imaginary number
- ⁇ 0 is the air conductivity
- w is the angular velocity corresponding to the frequency of the emission source
- ⁇ is the air permeability
- the H subscript indicates the direction of the electric dipole
- the superscript indicates the direction of the generated magnetic field
- ⁇ x, ⁇ y , ⁇ z are the x, y, z coordinate differences between the grid center point of the underground unit and the measurement point of the aeronautical device
- r is the distance between the grid center point of the underground unit and the receiving point of the aeronautical device.
- the transmitting device uses the x-direction magnetic dipole source, and the receiving device measures the x-direction magnetic field.
- the unit magnetic dipole transmitting source is in full space conditions.
- the following expression for the magnetic field generated at the measuring point is expression 13):
- Each parameter in expression 13) is consistent with each parameter in expression 7)-12).
- Each parameter in expression 14) is consistent with each parameter in expressions 1)-6).
- Use expression 14) Minus expression 15) The analytical solution of the secondary magnetic field generated by the emission source of the aerial electromagnetic method coaxial device at the measuring point can be obtained.
- the transmitting device uses the z-direction magnetic dipole source, and the receiving device measures the z-direction magnetic field.
- the magnetic field generated by the unit magnetic dipole transmitting source at the measuring point under the full space condition is expression 15):
- Each parameter in expression 15) is consistent with each parameter in expression 7)-12).
- Each parameter in expression 16) is consistent with each parameter in expression 1)-6).
- Use expression 16) Minus expression 15) The analytical solution of the secondary magnetic field generated by the emission source of the airborne electromagnetic coplanar device at the measuring point can be obtained.
- Figure 5 shows the process of gradually expanding the grid area so that the numerical solution of the secondary magnetic field at the measuring point gradually approaches the analytical solution, and obtaining the footprint area of the aviation electromagnetic device.
- the first numerical solution calculation uses a 400 ⁇ 400 ⁇ 100m 3 area grid.
- the network is expanded according to the ratio of 4:4:1 Grid area.
- This embodiment is based on the quadratic electric field vector finite element method, and the finite element discrete grid element edge center electric field satisfies expression 17):
- E p is the initial field emission sources uniformly generated in the ground, can be calculated by Expression 1) -6);
- K is the overall stiffness matrix
- S is the item related to the initial electric field generated by the emission source underground
- SH is the Poynting vector item generated by the electromagnetic field at the boundary of the calculation area
- K II is symmetric and sparse
- K IB is asymmetric
- the truncated boundary vector finite element method cuts the vector finite element method calculation area into anomalous body and its element thickness wrap.
- the relationship between E Bs and E Is in expression 20) is established through Green's function.
- the truncated boundary vector finite element method calculates the secondary electric field E Bs at the boundary of the region, which can be expressed by the internal element scattering current as expression 21):
- g ee is the connection between the scattered current J I at the center of the internal unit and the electric field E Bs at the boundary of the calculation area.
- E IP is the initial electric field at the edge of the internal unit
- V is the internal unit volume
- Ne is the linear interpolation operator
- ⁇ * is the conductivity of the uniform underground medium
- ⁇ is the conductivity of the underground medium.
- G ee is the expression 24):
- G ee g ee V( ⁇ - ⁇ * )n e 24);
- the magnetic field generated by the center current J I of the discrete element of the abnormal body in the footprint at the measuring point can be regarded as the magnetic field generated in the air by a uniform underground electric dipole source.
- the secondary magnetic field generated by the underground unit electric dipole source at the measuring point is expressed as expression 26)-28):
- g em is the Green's function linking the underground current in a uniform half-space with the magnetic field at the measuring point of the aeronautical device.
- the calculation of matrices g ee and g em in Expression 21) and Expression 32) includes Bessel integrals of the first kind, and such integral calculation requires a lot of time. Fortunately, the truncated boundary vector finite element method guided by the footprint of the present invention can be reused in subsequent stations after the matrix g ee and g em are calculated for the first station. As shown in Figure 7, the calculation area of the first station is a solid line frame, and the calculation area of the second station is a dashed frame.
- the relative positional relationship between the discrete units within the calculation area of the second station and the boundary of the calculation area can be determined from the The relative positional relationship between the discrete unit within the calculation area of a station and the boundary of the calculation area can be found accordingly.
- the second station can use the g ee of the first station.
- the positional relationship between the discrete unit in the calculation area of the second station and the measurement point of the second station can be found in the position relationship between the discrete unit and the measurement point of the first station in the calculation area of the first station.
- the second station can use the g em of the first station.
- subsequent stations can reuse the matrices g ee and g em generated in the calculation process of the first station, which greatly speeds up the calculation of the secondary magnetic field response of the subsequent station.
- the invention calculates a one-dimensional layered model, which has verified the calculation accuracy of the invention.
- the aviation electromagnetic method uses the device shown in Figure 2, the transmitting coil and the receiving coil are separated by 7.5m, and the transmitting source receiving device is 40m away from the ground.
- First use Figure 5 first use a uniform cubic grid to divide a 400 ⁇ 400 ⁇ 100m 3 calculation area, evaluate the secondary magnetic field generated by the discrete unit current in the area at the measuring point and the analytical solution error, so as to determine the area Whether to meet the footprint area requirements.
- the calculation results show that when using a coaxial device, the excitation current in a volume of 400 ⁇ 400 ⁇ 100m 3 in the x direction of the measuring point has a secondary magnetic field response with a numerical solution error of 0% for the real part and 2% for the imaginary part; when using the same-plane device , The excitation current in the volume of 400 ⁇ 400 ⁇ 100m 3 in the z direction of the measuring point, the real part error of the numerical solution of the secondary magnetic field response is 3%, and the imaginary part error is 2%. Therefore, it is judged that the volume of 400 ⁇ 400 ⁇ 100m 3 is the footprint area of the aerial electromagnetic method device.
- the dashed frame in Figure 6 is the footprint area of the aerial electromagnetic method device, and the internal line frame No.
- Table 3 is the calculation area of the layered medium model by the truncated boundary vector finite element method guided by the footprint.
- Table 1 shows the accuracy comparison between the numerical solution and the analytical solution of the truncated boundary vector finite element method (the method of the present invention) guided by the layered model footprint under different device forms. From the data in Table 1, it can be seen that under different device forms, the errors of the real and imaginary parts of the values are kept within 5%, which meets the accuracy requirements of the aviation electromagnetic method.
- Hsx is the calculation result of the magnetic field in the x direction
- Hsz is the calculation result of the magnetic field in the z direction.
- the present invention calculates a relatively complex three-dimensional model to verify the high efficiency and accuracy of the present invention for the calculation of complex geological models.
- the complex model is shown in Figure 8. There are two jagged low-resistance anomalies with different electrical conductivity underground.
- the calculation area is shown in box No. 1 in Fig. 8.
- the area size is 405 ⁇ 21 ⁇ 20m 3 , and the 3 ⁇ 3 ⁇ 4m 3 network is selected.
- the grid is divided into the area.
- the second station is located at 727.5m.
- the mesh of the first station is moved 4 units to the left to obtain the calculation area of the second station.
- subsequent stations move the grid one after another.
- the calculation area is shown in the second box in Figure 8.
- the calculation area of the model includes two anomalous conductive bodies and the wrapping layer.
- the calculation area is shown in the third box in Figure 8.
- the calculation area size is 834 ⁇ 21 ⁇ 20m 3
- Table 2 compares the calculation efficiency of the method of the present invention and the traditional truncated boundary vector finite element method in calculating the data of the single station of the model.
- the data in Table 2 show that the degree of freedom of the equations formed by the truncated boundary vector finite element method guided by the footprint, the required Green's function, and the calculation time are significantly less than the traditional truncated boundary vector finite element method.
- the Green's function between the electric field of the element inside the contact area and the electric field of the boundary node in the two methods can be stored in a matrix form.
- the truncated boundary vector finite element method guided by the footprints calculates the relative position of the area and the receiving device every time, it only needs to calculate the Green's function linking the electric field of the unit inside the calculation area with the magnetic field of the receiving point, which greatly reduces the calculation time.
- Figures 9a, 9b, 10a and 10b show the numerical comparison between the truncated boundary vector finite element method based on a relatively complex three-dimensional model and the footprint guided truncated boundary vector finite element method (the present invention).
- the calculated results of the secondary magnetic field response of the measuring point are plotted in parts per million (ppm) of the full-space magnetic field of the emission source, and the two methods are compared.
- the real and imaginary parts of the horizontal magnetic field (Hsx) and vertical magnetic field (Hsz) calculated by the footprint-guided truncated boundary vector finite element method are the same as those of the traditional truncated boundary vector finite element method.
- the results are consistent. This result further verifies that the calculation result of the present invention is correct and has universal applicability.
- Step 1 Take a cube of the same size with a side length of 4m as the unit body, and divide the grid into 50 grids in the x and y axis directions and 25 grids in the z axis direction with the underground projection of the emission source as the center. 100 ⁇ 100 ⁇ 25 scale grid, the calculation area size is 400 ⁇ 400 ⁇ 100m 3 .
- Step 2 Calculate the excitation current of the emission source in the center of the discrete grid within 400 ⁇ 400 ⁇ 100m 3 ; calculate the secondary magnetic field generated by the excitation current of the emission source in the center of each discrete grid at the receiving point, and combine the grids The magnetic field generated by the internal current at the measuring point is superimposed to form a numerical solution of the magnetic field at the measuring point.
- Step 3 Calculate the analytical solution of the magnetic field generated by the uniform half-space underground medium at the measuring point, and use the analytical solution as the standard to calculate the numerical solution error.
- the error is within 5%, and the footprint area of the aviation electromagnetic device is determined to be 400 ⁇ 400 ⁇ 100m 3 .
- Step 4 Cut the infinitely extending low-resistance plate-shaped model into a finite-volume plate-shaped body of 400 ⁇ 400 ⁇ 4m 3.
- the footprint-guided truncated boundary vector finite element method defines the calculation area as 400 ⁇ 400 ⁇ 12m 3 , and uses 4 ⁇ 4 ⁇ 4m 3 grid division.
- Step 5 Use the truncated boundary vector finite element method to calculate the scattered current at the center of the discrete element in the area of 400 ⁇ 400 ⁇ 12m 3.
- Step 6 Calculate the magnetic field generated by the scattered current from the center of the discrete unit in the 400 ⁇ 400 ⁇ 12m 3 area at the measuring point, and superimpose the magnetic field generated by each unit to obtain the numerical solution of the truncated boundary vector finite element method guided by the layered model footprint.
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
本发明提供一种足迹引导的高效航空电磁法数值模拟方法，该方法包括：使用均匀网格剖分航空电磁法探测区域；逐步增加均匀单元数量，确定航空电磁法观测装置的足迹大小；对第一个测站使用截断边界矢量有限元法计算其足迹内存在的散射电流，并将此步骤中计算得到的格林函数进行存储；计算上一步得到的足迹内散射电流在第一个测点的产生的电磁场响应，并将并将此步骤中计算得到的格林函数进行存储；对后续测站重复前两个步骤。本发明选择将探测区域按单个测站足迹逐次计算，极大地缩小了计算区域；使用均匀网格剖分，在后续测站使用截断边界矢量有限元法时能重复使用格林函数，很大程度上提高了后续测站的计算效率。
Description
本发明涉及航空电磁法正演技术领域，具体的涉及一种足迹引导的高效航空电磁法数值模拟方法。
航空电磁法观测装置使用机载平台，具有高效率的特点，探测区域往往覆盖几十平方公里到上千平方公里。传统数值模拟方法要求将全部探测区域作为计算区域，这种大区域网格计算将耗费大量的计算内存与时间。
截断边界矢量有限元法只要求将探测区域内异常体作为计算区域，具有效率高的特点。但是在面对大范围异常体存在时，需要进行大量的格林函数计算，使得该方法效率备受质疑。
因此，开发一种计算区域适当且计算效率高的模拟方法具有重要意义。
发明内容
为了解决航空电磁法探测区域范围大带来的计算困难，本发明提供了一种足迹引导的高效航空电磁法数值模拟方法，该方法结合航空电磁法观测装置具有有限足迹区域(即航空电磁法足迹特征)与截断边界矢量有限元法，将整个航空电磁法探测区域按测站细分为多个子区域逐一计算测站测点电磁场，具有计算区域适当且计算效率高的特征，具体技术方案如下：
本发明提供一种足迹引导的高效航空电磁法数值模拟方法，包括以下步骤：
步骤S100：使用均匀网格剖分航空电磁法探测区域；
步骤S200：计算所使用航空电磁法观测装置的足迹区域在x、y、z三个方向的延伸尺度；所述足迹区域在x与y方向网格数量相等，且z方向网格数量为x或y方向网格数量的1/4至1/2之间；网格具体数量由均匀半空间模型测点电磁场解析解与数值解比率决定；
步骤S300：使用截断边界矢量有限元法将第一个测站对应的足迹作为异常体，计算足迹内存在的散射电流，并将计算过程中产生的格林函数进行存储备用；
所述格林函数在使用截断边界矢量有限元时产生，是用来计算足迹内离散单元散射电流在截断边界矢量有限元法计算边界产生的电磁场；
步骤S400：使用第一测站足迹内离散单元散射电流计算测点电磁场响应，并对计算过程中产生的格林函数进行存储备用；
所述格林函数用来计算足迹内离散单元散射电流在测站观测点产生的电磁场；
步骤S500：对后续测站重复使用第一测站在S300与S400步骤所存储的格林函数，先按照S300计算测站足迹内散射电流，再按照S400计算站足迹内散射电流在测点电磁场。
对于上述本发明步骤S300和S400记载的格林函数可知，实质上属于两个不同用途的格林函数，因此，为了更好的区分，对于步骤S300中记载的格林函数，命名为第一格林函数，对于步骤S400中记载的格林函数，命名为第二格林函数。
以上技术方案中优选的，所述步骤S100中网格为规则的六面体网格。
以上技术方案中优选的，步骤S200中具体要求足迹区域以发射源在地表投影为中心，按4:4:1比例逐步增加x、y、z方向的网格数量，直到足迹内散射电流在测点产生的电磁场与测点解析解误差小于5％。
以上技术方案中优选的，所述步骤S200具体是：
步骤S210：计算航空电磁发射源在接收点产生的二次电磁场解析解；
计算航空电磁装置发射源正下方均匀半空间地下A区域的激发电流，并计算区域内激发电流在航空电磁测量装置产生的二次场数值解；此处A区域为400×400×100m
3的区域；
步骤S220：计算得到的二次磁场数值解与得到的二次磁场解析解的相对误差；
步骤S230：相对误差大小判断，若磁场实部与虚部相对误差均不超过5％，则步骤S210中的A区域体积定义为该航空电磁装置足迹；
若磁场实部与虚部相对误差均超过5％，则按照4:4:1比例扩展重新定义S210中A区域并计算区域内的激发电流和区域内激发电流在航空电磁测量装置产生的二次场数值解，返回步骤S220。
以上技术方案中优选的，步骤S300中不论测站下方异常电导率规模，只需要将足迹区域作为计算目标区域，使用截断边界矢量有限元法计算足迹内存在的散射电流。
以上技术方案中优选的，所述步骤S300中截断边界矢量有限元法的使用包括：
步骤S310：将计算区域定义为足迹及其一个单元厚度包裹层，电场定义在计算区域剖分单元棱边中心点，使用矢量有限元理论建立电场所满足的待求解方程组；
步骤S320：计算格林函数联系足迹内单元中心电流与计算区域边界电场，将格林函数排列为矩阵形式，行数为计算区域边界棱边数，列数为足迹区域内部单元棱边数；
步骤S330：将步骤S320关系方程组系数矩阵存储备用；
步骤S340：形成计算区域边界棱边中心电场由足迹内部单元棱边中心电场关系式方程组；
步骤S350：将步骤S310中计算区域边界电场使用步骤S320表达式替换，形成截断边界矢量有限元法控制方程组；
步骤S360：求解方程组，得到足迹内单元棱边中心电场值，使用线性插值得到单元中心散射电流。
以上技术方案中优选的，所述步骤S400具体包括以下步骤：
步骤S410：计算联系第一测站测点与第一测站足迹内离散单元散射电流格林函数，并排列为行矩阵，矩阵列数为散射电流个数，即足迹离散单元个数；
步骤S420：将步骤S410形成的行矩阵存储备用；
步骤S430：将格林函数行矩阵乘以散射电流形成的列矩阵，得到测点磁场值。
以上技术方案中优选的，所述步骤S500包括以下步骤：
步骤S510：重复步骤S310，调用步骤S330存储的格林函数，重复步骤S340-S360，得到后续测站足迹内单元散射电流；
步骤S520：调用步骤S420存储的格林函数，重复步骤S430，得到后续测站测点磁场值。
应用本发明的技术方案，效果是：
1、本发明提供的足迹引导的高效航空电磁法数值模拟方法，摒弃了传统数值模拟方法要求将所有探测区域作为计算区域，选择将探测区域按单个测站足迹逐次计算，极大的缩小了计算区域。
2、本发明提供的足迹引导的高效航空电磁法数值模拟方法，使用均匀网格剖分，所有测站足迹由相同网格组成。在后续测站使用截断边界矢量有限元法时，可以重复使用格林函数第一测站存储的格林函数，很大程度上提高了后续测站计算效率。
具体请参考根据本发明的一种足迹引导的高效航空电磁法数值模拟方法提出的各种实施例的如下描述，将使得本发明的上述和其他方面显而易见。
图1是本实施例中足迹引导的高效航空电磁法数值模拟方法的示意图；
图2是测站装置形式及下方均匀地下介质均匀剖分示意图(使用均匀网格剖分航空电磁法探测区域)；
图3a是同轴装置形式地下均匀介质内电流实部分布图；
图3b是同轴装置形式地下均匀介质内电流虚部分布图；
图4a是共面装置形式地下均匀介质内电流实部分布图；
图4b是共面装置形式地下均匀介质内电流虚部分布图；
图5是地下足迹区域尺度评估过程图(实线框是第一次评估使用400×400×100m
3计算区域，虚线框是后续按照4:4:1扩展后计算区域)；
图6是不同方法计算大型水平延伸薄板异常体所需计算区域示意图(一号框是传统微分法所需计算区域，二号框是传统截断边界矢量有限元法所需计算区域，三号框是足迹引导的截断边界矢量有限元法所需计算区域)；
图7是后续测站足迹区域示意图(第一测站计算区域如实线黑框所示，第二测站计算区域如虚线黑框所示)；
图8是锯齿复杂三维模型示意图(一号框是本实施例第一测站计算区域，二号框是本实施例第45测站计算区域，三号框是传统截断边界矢量有限元法在任意测站都需要的计算区域)；
图9a是使用同轴装置时本实施例方法与传统截断边界矢量有限元法在锯齿复杂三维模型的磁场数值解虚部对比图(黑点是本实施例方法计算结果虚部，黑线是传统截断边界矢量有限元法计算结果虚部)；
图9b是使用同轴装置时本实施例方法与传统截断边界矢量有限元法在锯齿复杂三维模型的磁场数值解实部对比图(圆圈是本实施例方法计算结果实部，虚线是传统截断边界矢量有限元法计算结果实部)；
图10a是使用共面装置时本实施例方法与传统截断边界矢量有限元法在锯齿复杂三维模型的磁场数值解虚部对比图(黑点是本实施例方法计算结果虚部，黑线是传统截断边界矢量有限元法计算结果虚部)；
图10b是使用共面装置时本实施例方法与传统截断边界矢量有限元法在锯齿复杂三维模型的磁场数值解实部对比图(圆圈是本实施例方法计算结果实部，虚线是传统截断边界矢量有限元法计算结果实部)。
构成本申请的一部分的附图用来提供对本发明的进一步理解，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。
实施例：
见图1，本发明提供的一种足迹引导的高效航空电磁法数值模拟方法，包括以下步骤：
步骤S100：使用均匀网格剖分航空电磁法探测区域，优选的，此处使用的均匀网格为规则的六面体网格，网格x、y和z方向尺度相等，大小等于在5、10和25三个数之间选取，尺度越小精度要高。
步骤S200：计算所使用航空电磁法观测装置的足迹区域在x、y、z三个方向的延伸尺度，具体包括如下步骤：
第一步：计算单测站发射源在测点产生的均匀半空间解析解与全空间解析解，使用半空间解析解减去全空间解析解得到均匀地下介质内电流在测点产生的二次磁场解析解；
第二步：计算航空电磁装置发射源在其正下方均匀地下A区域(400×400×100m
3区域)均匀网格单元中心位置产生的激发电流；
第三步：计算所定义体积区域(400×400×100m
3区域)均匀网格单元内激发电流在测量装置处产生的二次场进行叠加，得到测点二次磁场数值解；
第四步：计算第三步得到的二次磁场数值解与第一步得到的二次磁场解析解的相对误差；
第五步：相对误差大小判断，若磁场实部与虚部相对误差均不超过5％，则步骤第二步中的A区域体积定义为该航空电磁装置足迹；
若磁场实部与虚部相对误差均超过5％，则按照4:4:1比例扩展(按照x、y与z方向)扩展均匀单元数量，重新定义第二步中A区域并计算区域内的激发电流和区域内激发电流在航空电磁测量装置产生的二次场数值解，返回步骤第三步。
步骤S300：使用截断边界矢量有限元法将第一个测站对应的足迹作为目标区域， 计算足迹内存在的散射电流，并将计算过程中产生的格林函数进行存储备用；
本实施例中截断边界矢量有限元法的使用包括以下步骤：
步骤S310：将计算区域定义为足迹及其一个单元厚度包裹层，电场定义在计算区域剖分单元棱边中心点，使用矢量有限元理论建立电场所满足的待求解方程组；
步骤S320：计算格林函数联系足迹内单元中心电流与计算区域边界电场，将格林函数排列为矩阵形式，行数为计算区域边界棱边数，列数为足迹区域内部单元棱边数；
步骤S330：将步骤S320关系方程组系数矩阵存储备用；
步骤S340：形成计算区域边界棱边中心电场由足迹内部单元棱边中心电场关系式方程组；
步骤S350：将步骤S310中计算区域边界电场使用步骤S320表达式替换，形成截断边界矢量有限元法控制方程组；
步骤S360：求解方程组，得到足迹内单元棱边中心电场值，使用线性插值得到单元中心散射电流。
步骤S400：使用测站足迹内离散单元散射电流计算测点电磁场响应，步骤如下：
步骤S410：计算联系第一测站测点与第一测站足迹内离散单元散射电流格林函数，并排列为行矩阵，矩阵列数为散射电流个数，即足迹离散单元个数；
步骤S420：将步骤S410形成的行矩阵存储备用；
步骤S430：将格林函数行矩阵乘以散射电流形成的列矩阵，得到测点磁场值。
步骤S500：对后续测站重复使用第一测站在S300与S400步骤所存储的格林函数，先按照S300计算测站足迹内散射电流，再按照S400计算站足迹内散射电流在测点电磁场。
优选的，具体S500包括以下步骤：
步骤S510：重复步骤S310，调用步骤S330存储的格林函数，重复步骤S340-S360，得到后续测站足迹内单元散射电流；
步骤S520：调用步骤S420存储的格林函数，重复步骤S430，得到后续测站测点磁场值。
应用本实施例的方案，详情如下：
1、使用均匀网格剖分探测区域
如图2所示，航空电磁法使用机载平台搭载观测装置，观测装置由发射线圈与接 收线圈组成。当发射线圈与接收线圈为水平线圈时称为共面装置形式，当发射线圈与接收线圈为垂直线圈时称为同轴装置形式。航空电磁法距离地表高度在40m之间，航线从右至左，发射线圈源位于接收线圈前方7.5m，发射线圈发射频率为10000hz。为了计算该航空电磁法装置在探测区域各个测站的磁场数据，本实施例使用均匀网格地下介质。
2、计算航空电磁观测装置足迹区域
航空电磁法发射装置所使用的垂直线圈时(同轴装置)可以当做水平磁偶极子源处理，发射装置所使用的水平线圈时(共面装置)可以当做垂直磁偶极子源处理。在使用时谐因子e
iwt时，设置直角坐标，地表为xy面，z坐标向下。
2.1、均匀地下初始电场求取，如下：
空气中发射源为单位水平磁偶极子源(同轴装置)在接收点(地下网格单元中心)产生的电场表达式满足表达式1)-3)：
E
z＝0 3)；
空气中单位垂直磁偶极子源(共面装置)在地下产生的电场表达式满足表达式4)-6)：
E
z＝0 6)；
表达式1)-6)中：i是单位虚数，w是与发射源频率对应的角速度，Δx与Δy分别是接收点(地下网格中心点)与航空电磁发射源的x与y坐标差，r与ρ分别是接收点(地下网格中心点)与偶极子源之间的距离与水平投影距离，λ是空间波数，J
0与J
1分别是第一类0阶与1阶贝塞尔函数，z与z'分别为接收点(地下网格中心点)与发射源点的纵坐标，μ是空气磁导率，
σ是地下均匀介质电导率。
2.2、均匀地下发射源产生电流分布，如下：
图3a、图3b、图4a和图4b分别展示了水平磁偶极子源(同轴装置)与垂直磁偶极子源(共面装置)电流幅值分布的水平切片图与竖直切片图。图3a与图4a第一列展示了地下z＝20m剖面，地下水平面电流分布，电流大于10
-9.5A/m主要位于400×400m
2水平区域内。图3b与图4b第二列展示了源水平偶极子方向竖直切面电流分布，电流大于10
-9.5A/m主要在地下100m深度内。图3a、图3b、图4a和图4b表明了航空电磁法发射源引起的地下电流主要集中在发射源正下方400×400×100m
3区域内，在足迹评估过程中，首先选定该区域计算测点二次磁场数值解。
2.3、测点二次磁场数值解与解析解求取，如下：
本实施例选择均匀半空间模型，通过比较两种不同方法计算地下介质在测点产生的磁场来定义航空电磁法装置足迹区域大小，具体是：首先，计算得到发射源在地下选定区域引起的电流分布之后；然后，计算选定区域网格内电流在航空电磁观测点引起的二次磁场叠加得到测点二次磁场数值解；最后，得到的二次磁场数值解与解析解相比较，判断所选电流区域是否已达到所要求规模。其中判断准则为以二次磁场解析解为标准，二次磁场数值解误差小于5％，解析解由测点磁场均匀半空间解析解减去全空间解析解得到。
在求取所选区域内网格二次磁场数值解时，网格单元内部电流被视为均匀分布，在计算其在测点产生的磁场时，网格单元被当做电偶极子源。在航空电磁法同轴装置时，单位电偶极子源在全空间产生的磁场为表达式7)-9)：
使用表达式1)-3)计算得到地下网格单元中心电场分布后，使用表达式7)-9)计算每个单元中心电流在测点产生的磁场，通过累加所选区域内单元在测点磁场贡献得到测点二次磁场数值解。
在航空电磁法共面装置时，单位电偶极子源在全空间产生的磁场为表达式10)-12)：
使用表达式4)-6)计算得到地下网格单元中心电场分布后，使用表达式10)-12)计算每个单元中心电流在测点产生的磁场，通过累加所选区域内单元在测点磁场贡献得到测点二次磁场数值解。
表达式7)-12)中：
i是单位虚数；σ
0是空气电导率；w是与发射源频率对应的角速度；μ是空气磁导率；H下标表示电偶极子方向，上标表示产生的磁场方向；Δx、Δy、Δz分别是地下单元网格中心点与航空装置测量点之间的x、y、z坐标差；r是地下单元网格中心点与航空装置接收点之间的距离。
在求取所选区域内网格二次磁场解析解时，在同轴装置下即发射装置使用x方向磁偶极子源，接收装置测量x方向磁场，单位磁偶极发射源在全空间条件下在测点产生的磁场表达式为表达式13)：
表达式13)中各参数与表达式7)-12)中各参数一致。
在均匀半空间条件下单位磁偶极发射源在测点产生的磁场表达式为
表达式14)中各参数与表达式1)-6)中各参数一致。使用表达式14)中
减去表达式15)中
可以得到航空电磁法同轴装置发射源在测量点产生的二次磁场解析解。
在共面装置下即发射装置使用z方向磁偶极子源，接收装置测量z方向磁场，单位磁偶极发射源在全空间条件下在测点产生的磁场为表达式15)：
表达式15)中各参数与表达式7)-12)中各参数一致。
在均匀半空间条件下单位磁偶极发射源在测点产生的磁场为表达式16)：
表达式16)中各参数与表达式1)-6)中各参数一致。使用表达式16)中
减去表达式15)中
可以得到航空电磁法共面装置发射源在测量点产生的二次磁场解析解。
2.4、自适应扩展网格区域求取足迹范围，如下：
图5展示了逐步扩大网格区域，使得测点二次磁场数值解逐渐接近解析解，求得航空电磁装置足迹区域的过程。第一次数值解计算使用400×400×100m
3区域内网格， 当区域网格内电流在测点二次磁场数值解与解析解差异大于5％时，按4:4:1比例扩展网格区域。
3、使用截断边界矢量有限元法计算第一测站足迹内散射电流，如下：
本实施例基于二次电场矢量有限元法，有限元离散网格单元棱边中心电场满足表达式17)：
其中：
是哈密顿算符，E
p是发射源在均匀地下产生的初始电场，可以由表达式1)-6)计算；E
s是地下介质电导率异常体在离散网格上产生的二次电场，σ
*是均匀地下介质电导率；σ是地下介质电导率，其他参数同表达式1)-6)中参数的含义。
根据Liu,R.,R.W.Guo,J.X.Liu,C.Y.Ma,and Z.W.Guo,“A hybrid solver based on the integral equation method and vector finite-element method for 3D controlled-source electromagnetic method modeling,”Geophysics,vol.83,no.5,pp.E319-E333,Sep.2018.矢量有限元离散网格中二次电场E
s满足表达式18)：
KE
s＝S+S
H 18)；
其中：K是整体刚度矩阵，S是与发射源在地下产生的初始电场相关项，S
H是电磁场在计算区域边界产生的坡印廷矢量项；当计算区域内部单元用I标示，计算区域边界单元用B标示，表达式18)可以分解为表达式19)：
其中：
是关于内部单元的电磁场积分项，根据电磁场在电导率边界的切向连续性，
内部会相互抵消。使用表达式19)的上半部分，计算区域内部二次电场满足表达式20)：
K
IIE
Is+K
IBE
Bs＝S
I 20)；
其中：K
II对称且稀疏；K
IB是不对称的。图6一号线框表示了矢量有限元法的计算区域，该计算区域边界远离异常体(水平黑色块状体)，使得E
Bs＝0。
如图6二号框所示，截断边界矢量有限元法将矢量有限元法计算区域截断为异常体及其单元厚度包裹层。表达式20)中的E
Bs与E
Is通过格林函数建立联系，截断边界矢量有限元法计算区域边界二次电场E
Bs可以由内部单元散射电流表示为表达式21)：
E
Bs＝g
eeJ
I 21)；
其中：g
ee是联系内部单元中心处散射电流J
I与计算区域边界电场E
Bs。通过将内部单 元棱边电场E
Is使用线性插值算子到单元中心，中心处散射电流J
I可表示为：
J
I＝V(σ-σ
*)N
e(E
Is+E
IP) 22)；
其中：E
IP是内部单元棱边初始电场，V是内部单元体积，N
e是线性插值算子，σ
*是均匀地下介质电导率；σ是地下介质电导率。然后，计算区域边界电场E
Bs可以由内部单元电场矩阵形式表达为表达式23)：
E
Bs＝G
ee(E
Is+E
IP) 23)；
其中：G
ee为表达式24)：
G
ee＝g
eeV(σ-σ
*)n
e 24)；
将表达式23)代入表达式20)中可得，计算区域内部单元二次电场控制方程组为表达式25)：
(K
II+K
IBG
ee)E
Is＝S
I-K
IBG
eeE
IP 25)；
通过求解方程组表达式25)，在得到E
Is后，使用表达式22)得到足迹内离散单元中心处电流强度J
I。
4、计算足迹区域内散射电流在第一测站测点的磁场，如下：
足迹内异常体离散单元中心电流J
I在测点产生的磁场可以视为，均匀地下电偶极源在空气中产生的磁场。使用共轴装置时，地下单位电偶极源在测点产生的二次磁场表达为为表达式26)-28)：
使用同面装置时，地下单位电偶极源在测点产生的二次磁场表达为表达式29)-31)：
其中：H下标表示电偶极子方向，上标表示产生的磁场方向。表达式26)-32)中各参数含义与表达式1)-6)一致。最终，测点二次磁场表达式可以写为矩阵形式表达式32)：
H
s＝g
emJ
I 32)；
其中g
em是联系均匀半空间地下电流与航空装置测点磁场的格林函数。
5、后续测站测点磁场计算，如下：
表达式21)与表达式32)中的矩阵g
ee与g
em计算包含第一类贝塞尔积分，此类积分计算需要耗费大量的时间。幸运的是，本发明足迹引导的截断边界矢量有限元法在对第一测站进行矩阵g
ee与g
em计算后，后续测站可以重复使用。如图7所示，第一测站计算区域为实线框，第二测站计算区域为虚线框，第二测站计算区域内部离散单元与计算区域边界之间相对位置关系，均可以从第一测站计算区域内部离散单元与计算区域边界之间相对位置关系中相应找到，鉴于此第二测站可以使用第一测站的g
ee。同样的，第二测站计算区域内离散单元与第二测站测点之间位置关系可以在第一测站计算区域内离散单元与第一测站测点之间位置关系中相应找到。鉴于此第二测站可以使用第一测站的g
em。以此类推，后续测站可以重复使用第一测站计算过程中产生的矩阵g
ee与g
em，由此很大程度上加快后续测站二次磁场响应计算。
以下结合具体实例对本发明的方法在不同发射装置及模型的应用进行详细说明。
本发明计算了一维层状模型，已验证本发明计算精度。航空电磁法使用图2所示装置，发射线圈与接收线圈相距7.5m，发射源接收装置距离地面40m。首先使用图5所示，首先使用均匀立方体网格剖分400×400×100m
3大小的计算区域，评估该区域内离散单元电流在测点产生的二次磁场与解析解误差，从而确定该区域是否满足足迹区域要求。计算结果表明在使用共轴装置时，400×400×100m
3体积中的激发电流在测点x方向二次磁场响应数值解实部误差0％，虚部误差2％；在使用同面装置时，400×400×100m
3体积中的激发电流在测点z方向二次磁场响应数值解实部误差3％，虚部误差2％。由此判断400×400×100m
3体积为该航空电磁法装置的足迹区域。图6虚线框为该航空电磁法装置足迹区域，其内部三号线框为足迹引导的截断边界矢量有限元法对层状介质模型的计算区域。表1给出了层状模型足迹引导的截断边界矢量有限元法(本发明方法)数值解与解析解在不同装置形式下的精度对比。由表1数据可知不同装置形式下，数值实部与虚部误差都保持在5％以内，符合航空电磁法精度要求。其中Hsx是x方向磁场计算结果，Hsz是z方向磁场计算结果。
表1 足迹引导的截断边界矢量有限元法数值解与解析解对比
| 同轴装置 | 实部(Hsx) | 虚部(Hsx) | 实部(Hsz) | 虚部(Hsz) |
| 截断边界矢量有限元法 | -3.899e-9 | 2.238e-9 | -8.032e-10 | 2.198e-10 |
| 解析解 | -3.940e-9 | 2.346e-9 | -8.112e-10 | 2.313e-10 |
| 相对误差(％) | -1.0 | -4.6 | -0.9 | -4.9 |
| 共面装置 | 实部(Hsx) | 虚部(Hsx) | 实部(Hsz) | 虚部(Hsz) |
| 截断边界矢量有限元法 | 8.117e-10 | -2.235e-10 | -7.977e-9 | 4.500e-9 |
| 解析解 | 8.119e-10 | -2.313e-10 | -7.928e-9 | 4.694e-9 |
| 相对误差(％) | 0.0 | -3.3 | 0.6 | -4.1 |
本发明计算了相对复杂三维模型，以验证本发明对复杂地质模型计算的高效性与准确性。复杂模型如图8所示，地下存在两个电导率不同的锯齿状低阻异常体。航空电磁法图2所示装置，第一测站位于x＝739.5m最后一个测站位于x＝211.5m，共45个测站，相邻测站相距12m。由于本次航空电磁法装置与层状模型一致，故装置的足迹区域同样为400×400×100m
3。本次航空电磁法第一测站使用足迹引导的截断边界矢量有限元法时，计算区域如图8一号框所示，区域大小为405×21×20m
3，选择3×3×4m
3网格剖分该区域。第二测站位于727.5m，使用足迹引导的截断边界矢量有限元法时，将第一测站剖分网格向左移动4个单元，得到第二测站计算区域。同理，后续测站陆续移动网格，在45号测站时，其计算区域如图8二号框所示。
在使用传统截断边界矢量有限元法计算该模型时，该模型的计算区域包括两个导电异常体和包裹层，计算区域如图8三号框所示，计算区域尺度为834×21×20m
3，表2对比了本发明方法与传统截断边界矢量有限元法在计算该模型单测站数据的计算效率。表2数据表明足迹引导的截断边界矢量有限元法形成的方程组自由度、需要的格林函数、计算时间明显小于传统截断边界矢量有限元法。两种方法在计算一次联系区域内部单元电场与边界节点电场之间的格林函数后，可以将其以矩阵形式存储。然而，由于足迹引导的截断边界矢量有限元法每次计算区域与接收装置相对位置不变，只需要计算一次联系计算区域内部单元电场与接收点磁场的格林函数，大大减少了计算时间。
图9a、图9b、图10a和图10b显示了基于相对复杂三维模型的截断边界矢量有限元法和足迹引导截断边界矢量有限元法(本发明)之间的数值比较。计算得到的测点二次磁场响应结果以发射源全空间磁场的百万分之一(ppm)绘制，并对这两种方法进行比较。对于同轴装置和共面装置，足迹引导的截断边界矢量有限元法计算得到的水平方向磁场(Hsx)与垂直方向磁场(Hsz)的实部与虚部结果与传统截断边界矢量有限元法的结果一致。该结果进一步验证了本发明计算结果正确，具有普遍适用性。
表2 本发明足迹引导的截断边界矢量有限元法与传统截断边界矢量有限元法计算效率对比
| 格林函数数量 | 方程组自由度 | 计算时间 | |
| 截断边界矢量有限元法 | 55823760 | 8280 | 147s |
| 足迹引导截断边界矢量有限元法 | 13206900 | 3990 | 27s |
本发明算法计算层状模型电磁响应具体实施如下：
步骤一：以边长为4m的相同大小立方体作为单元体，以发射源地下投影为中心分别在x与y轴方向将剖分网格50个，往z轴方向剖分网格25个，形成100×100×25规模网格，计算区域大小为400×400×100m
3。
步骤二：计算400×400×100m
3内，离散剖分网格中心发射源激发电流大小；计算各离散剖分网格中心发射源激发电流在接收点产生的二次磁场，并将各网格内电流在测点产生的磁场叠加，形成测点磁场数值解。
步骤三：计算均匀半空间地下介质在测点产生的磁场解析解，以解析解为标准，计算数值解误差。得到误差在5％以内，确定航空电磁法装置足迹区域大小为400×400×100m
3。
步骤四：将无限延伸低阻板状模型，截断为400×400×4m
3的有限体积板状体。足迹引导的截断边界矢量有限元法将计算区域定为400×400×12m
3，使用4×4×4m
3网格剖分，
步骤五：使用截断边界矢量有限元法计算400×400×12m
3区域内离散单元中心散射电流。
步骤六：计算400×400×12m
3区域内离散单元中心散射电流在测点产生的磁场，将每个单元产生的磁场叠加得到层状模型足迹引导的截断边界矢量有限元法数值解。
本领域技术人员将清楚本发明的范围不限制于以上讨论的示例，有可能对其进行若干改变和修改，而不脱离所附权利要求书限定的本发明的范围。尽管己经在附图和说明书中详细图示和描述了本发明，但这样的说明和描述仅是说明或示意性的，而非限制性的。本发明并不限于所公开的实施例。
通过对附图，说明书和权利要求书的研究，在实施本发明时本领域技术人员可以理解和实现所公开的实施例的变形。在权利要求书中，术语“包括”不排除其他步骤或元素，而不定冠词“一个”或“一种”不排除多个。在彼此不同的从属权利要求中引用的某些措施的事实不意味着这些措施的组合不能被有利地使用。权利要求书中的任何参考标记不构成对本发明的范围的限制。
Claims (8)
- 一种足迹引导的高效航空电磁法数值模拟方法，其特征在于，包括以下步骤：步骤S100：使用均匀网格剖分航空电磁法探测区域；步骤S200：计算所使用航空电磁法观测装置的足迹区域在x、y、z三个方向的延伸尺度；所述足迹区域在x与y方向网格数量相等，z方向网格数量在x或y网格数量的1/4至1/2之间，网格具体数量由均匀半空间模型测点电磁场解析解与数值解比率决定；步骤S300：使用截断边界矢量有限元法将第一个测站对应的足迹作为异常体，计算足迹内存在的散射电流，并将计算过程中产生的第一格林函数进行存储备用；所述第一格林函数在使用截断边界矢量有限元时产生，是用来计算足迹内离散单元散射电流在截断边界矢量有限元法计算边界产生的电磁场；步骤S400：使用第一测站足迹内离散单元散射电流计算测点电磁场响应，并对计算过程中产生的第二格林函数进行存储备用；所述第二格林函数用来计算足迹内离散单元散射电流在测站观测点产生的电磁场；步骤S500：对后续测站重复使用第一测站在S300与S400步骤所存储的格林函数，先按照S300计算测站足迹内存在的散射电流，再按照S400计算站足迹内离散单元散射电流在测点电磁场。
- 根据权利要求1所述的一种足迹引导的高效航空电磁法数值模拟方法，其特征在于，所述步骤S100中网格为规则六面体网格。
- 根据权利要求1所述的一种足迹引导的高效航空电磁法数值模拟方法，其特征在于，步骤S200中：要求足迹区域以发射源在地表投影为中心，按4:4:1比例逐步增加x、y、z方向的网格数量，直到足迹内散射电流在测点产生的电磁场与测点解析解误差不超过5％。
- 根据权利要求3所述的一种足迹引导的高效航空电磁法数值模拟方法，其特征在于，所述步骤S200具体是：步骤S210：计算航空电磁装置发射源在接收点产生的二次电磁场解析解；计算航空电磁装置发射源正下方均匀地下A区域的激发电流，并计算A区域内激发电流在航空电磁测量装置产生的二次电磁场数值解；此处A区域为400×400×100m 3的区域；步骤S220：计算得到的二次电磁场数值解与得到的二次电磁场解析解的相对误差；步骤S230：相对误差大小判断，若磁场实部与虚部相对误差均不超过5％，则步骤S210中的A区域体积定义为该航空电磁装置足迹；若磁场实部与虚部相对误差均超过5％，则按照4:4:1比例扩展重新定义S210中A区域并重新计算A区域内的激发电流和A区域内激发电流在航空电磁测量装置产生的二次电磁场数值解，返回步骤S220。
- 根据权利要求4所述的一种足迹引导的高效航空电磁法数值模拟方法，其特征在于，步骤S300中不论测站下方异常电导率规模，只需要将足迹区域作为计算目标区域，使用截断边界矢量有限元法计算足迹内存在的散射电流。
- 根据权利要求4所述的一种足迹引导的高效航空电磁法数值模拟方法，其特征在于，所述步骤S300中截断边界矢量有限元法的使用包括以下步骤：步骤S310：将计算区域定义为足迹及其一个单元厚度包裹层，电场定义在计算区域剖分单元棱边中心点，使用矢量有限元理论建立电场所满足的待求解方程组；步骤S320：计算第一格林函数，联系足迹内单元中心电流与计算区域边界电场，并将以第一格林函数排列为矩阵形式，行数为计算区域边界棱边数，列数为足迹区域内部单元棱边数；步骤S330：将步骤S320得到的第一格林函数矩阵存储备用；步骤S340：形成计算区域边界棱边中心电场与足迹内部单元棱边中心电场关系式方程组；步骤S350：将步骤S310中计算区域边界电场使用步骤S340得到的关系方程组表达式替换，形成截断边界矢量有限元法控制方程组；步骤S360：求解方程组，得到足迹内单元棱边中心电场值，使用线性插值得到单元中心散射电流。
- 根据权利要求6所述的一种足迹引导的高效航空电磁法数值模拟方法，其特征在于，步骤S410：计算联系第一测站测点与第一测站足迹内离散单元散射电流第二格林函数，并排列为行矩阵，矩阵列数为散射电流个数，即足迹离散单元个数；步骤S420：将步骤S410形成的第二格林函数行矩阵存储备用；步骤S430：将第二格林函数行矩阵乘以散射电流形成的列矩阵，得到测点磁场值。
- 根据权利要求7所述的一种足迹引导的高效航空电磁法数值模拟方法，其特征在于，所述步骤S500包括以下步骤：步骤S510：重复步骤S310，调用步骤S330存储的第一格林函数矩阵，重复步骤S340-S360，得到后续测站足迹内单元散射电流；步骤S520：调用步骤S420存储的第二格林函数行矩阵，重复步骤S430，得到后续测站测点磁场值。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| ZA2021/08251A ZA202108251B (en) | 2019-10-12 | 2021-10-26 | Numerical simulation method for footprint-guided high-efficiency airborne electromagnetic survey |
Applications Claiming Priority (2)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910966894.2A CN110598367A (zh) | 2019-10-12 | 2019-10-12 | 一种足迹引导的高效航空电磁法数值模拟方法 |
| CN201910966894.2 | 2019-10-12 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| WO2021068527A1 true WO2021068527A1 (zh) | 2021-04-15 |
Family
ID=68866887
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| PCT/CN2020/093495 WO2021068527A1 (zh) | 2019-10-12 | 2020-05-29 | 一种足迹引导的高效航空电磁法数值模拟方法 |
Country Status (3)
| Country | Link |
|---|---|
| CN (1) | CN110598367A (zh) |
| WO (1) | WO2021068527A1 (zh) |
| ZA (1) | ZA202108251B (zh) |
Cited By (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113569447A (zh) * | 2021-07-06 | 2021-10-29 | 武汉市市政建设集团有限公司 | 一种基于舒尔补方法的瞬变电磁三维快速正演方法 |
| CN113886950A (zh) * | 2021-09-19 | 2022-01-04 | 中国航空工业集团公司西安飞机设计研究所 | 一种机载设备质量特性仿真方法 |
| CN113887106A (zh) * | 2021-10-11 | 2022-01-04 | 吉林大学 | 一种基于Chikazumi模型的感应-磁化效应三维数值模拟方法 |
| CN114065585A (zh) * | 2021-11-22 | 2022-02-18 | 中南大学 | 一种基于库伦规范的三维电性源数值模拟方法 |
| CN114065586A (zh) * | 2021-11-22 | 2022-02-18 | 中南大学 | 一种三维大地电磁空间-波数域有限元数值模拟方法 |
| CN115906559A (zh) * | 2022-10-31 | 2023-04-04 | 重庆大学 | 一种基于混合网格的大地电磁自适应有限元正演方法 |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110598367A (zh) * | 2019-10-12 | 2019-12-20 | 中南大学 | 一种足迹引导的高效航空电磁法数值模拟方法 |
Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102565862A (zh) * | 2011-12-16 | 2012-07-11 | 朱德兵 | 一种瞬变电磁响应信号梯度测量方法及观测装置 |
| CN106199742A (zh) * | 2016-06-29 | 2016-12-07 | 吉林大学 | 一种频率域航空电磁法2.5维带地形反演方法 |
| CN106980736A (zh) * | 2017-04-11 | 2017-07-25 | 吉林大学 | 一种各向异性介质的海洋可控源电磁法有限元正演方法 |
| CN107121706A (zh) * | 2017-05-08 | 2017-09-01 | 厦门大学 | 基于波恩迭代法的航空瞬变电磁电导率三维反演方法 |
| US20190195067A1 (en) * | 2017-12-26 | 2019-06-27 | Saudi Arabian Oil Company | Determining sand-dune velocity variations |
| CN110598367A (zh) * | 2019-10-12 | 2019-12-20 | 中南大学 | 一种足迹引导的高效航空电磁法数值模拟方法 |
Family Cites Families (13)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6997097B2 (en) * | 2001-05-02 | 2006-02-14 | Lockheed Martin Corporation | Autonomous mission profile planning |
| JP2003043157A (ja) * | 2001-07-31 | 2003-02-13 | Hisatoshi Konishi | 空中電磁探査法のドリフト補正方法 |
| US6870370B2 (en) * | 2002-07-10 | 2005-03-22 | Agri Dynamics, Inc. | Electromagnetic induction detection system |
| CN101710187B (zh) * | 2009-12-17 | 2013-01-09 | 成都理工大学 | 一种时间域航空电磁高度校正方法 |
| CN101915943B (zh) * | 2010-08-10 | 2012-11-07 | 中南大学 | 均匀背景介质的介电常数和隐蔽目标参数的联合反演方法 |
| CN102043759A (zh) * | 2010-12-31 | 2011-05-04 | 中国航空工业集团公司第六三一研究所 | 一种数学模型数值计算程序的验证方法 |
| US8854255B1 (en) * | 2011-03-28 | 2014-10-07 | Lockheed Martin Corporation | Ground moving target indicating radar |
| US9542359B2 (en) * | 2011-12-29 | 2017-01-10 | Technoimaging, Llc | Method of subsurface imaging using superposition of sensor sensitivities from geophysical data acquisition systems |
| CN106199697A (zh) * | 2016-06-29 | 2016-12-07 | 中国石油化工股份有限公司 | 模拟微地震的弹性波正演方法 |
| CN108509693B (zh) * | 2018-03-13 | 2019-08-06 | 中南大学 | 三维频率域可控源数值模拟方法 |
| CN108984818A (zh) * | 2018-05-22 | 2018-12-11 | 吉林大学 | 固定翼时间域航空电磁数据拟三维空间约束整体反演方法 |
| CN110068873B (zh) * | 2019-05-10 | 2020-09-25 | 成都理工大学 | 一种基于球坐标系的大地电磁三维正演方法 |
| CN110210129B (zh) * | 2019-06-03 | 2021-05-11 | 中南大学 | 自适应有限元gpr频率域正演方法 |
-
2019
- 2019-10-12 CN CN201910966894.2A patent/CN110598367A/zh active Pending
-
2020
- 2020-05-29 WO PCT/CN2020/093495 patent/WO2021068527A1/zh active Application Filing
-
2021
- 2021-10-26 ZA ZA2021/08251A patent/ZA202108251B/en unknown
Patent Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102565862A (zh) * | 2011-12-16 | 2012-07-11 | 朱德兵 | 一种瞬变电磁响应信号梯度测量方法及观测装置 |
| CN106199742A (zh) * | 2016-06-29 | 2016-12-07 | 吉林大学 | 一种频率域航空电磁法2.5维带地形反演方法 |
| CN106980736A (zh) * | 2017-04-11 | 2017-07-25 | 吉林大学 | 一种各向异性介质的海洋可控源电磁法有限元正演方法 |
| CN107121706A (zh) * | 2017-05-08 | 2017-09-01 | 厦门大学 | 基于波恩迭代法的航空瞬变电磁电导率三维反演方法 |
| US20190195067A1 (en) * | 2017-12-26 | 2019-06-27 | Saudi Arabian Oil Company | Determining sand-dune velocity variations |
| CN110598367A (zh) * | 2019-10-12 | 2019-12-20 | 中南大学 | 一种足迹引导的高效航空电磁法数值模拟方法 |
Cited By (11)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113569447A (zh) * | 2021-07-06 | 2021-10-29 | 武汉市市政建设集团有限公司 | 一种基于舒尔补方法的瞬变电磁三维快速正演方法 |
| CN113569447B (zh) * | 2021-07-06 | 2023-11-03 | 武汉市市政建设集团有限公司 | 一种基于舒尔补方法的瞬变电磁三维快速正演方法 |
| CN113886950A (zh) * | 2021-09-19 | 2022-01-04 | 中国航空工业集团公司西安飞机设计研究所 | 一种机载设备质量特性仿真方法 |
| CN113886950B (zh) * | 2021-09-19 | 2022-09-06 | 中国航空工业集团公司西安飞机设计研究所 | 一种机载设备质量特性仿真方法 |
| CN113887106A (zh) * | 2021-10-11 | 2022-01-04 | 吉林大学 | 一种基于Chikazumi模型的感应-磁化效应三维数值模拟方法 |
| CN113887106B (zh) * | 2021-10-11 | 2024-04-12 | 吉林大学 | 一种基于Chikazumi模型的感应-磁化效应三维数值模拟方法 |
| CN114065585A (zh) * | 2021-11-22 | 2022-02-18 | 中南大学 | 一种基于库伦规范的三维电性源数值模拟方法 |
| CN114065586A (zh) * | 2021-11-22 | 2022-02-18 | 中南大学 | 一种三维大地电磁空间-波数域有限元数值模拟方法 |
| CN114065585B (zh) * | 2021-11-22 | 2024-05-10 | 中南大学 | 一种基于库伦规范的三维电性源数值模拟方法 |
| CN115906559A (zh) * | 2022-10-31 | 2023-04-04 | 重庆大学 | 一种基于混合网格的大地电磁自适应有限元正演方法 |
| CN115906559B (zh) * | 2022-10-31 | 2023-08-15 | 重庆大学 | 一种基于混合网格的大地电磁自适应有限元正演方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN110598367A (zh) | 2019-12-20 |
| ZA202108251B (en) | 2022-01-26 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| WO2021068527A1 (zh) | 一种足迹引导的高效航空电磁法数值模拟方法 | |
| CN106199742B (zh) | 一种频率域航空电磁法2.5维带地形反演方法 | |
| CN105223480B (zh) | 天线阵列时差法定位变电站局部放电源的定位误差仿真方法 | |
| CN110058315B (zh) | 一种三维各向异性射频大地电磁自适应有限元正演方法 | |
| CN113553748B (zh) | 一种三维大地电磁正演数值模拟方法 | |
| US9959670B2 (en) | Method for rendering terrain | |
| CN112685928B (zh) | 一种基于三相电抗器声源模型的噪声预测方法及*** | |
| CN109872394B (zh) | 基于最小二乘支持向量机的狭长三角形网格优化方法 | |
| CN112363236B (zh) | 一种基于pde的重力场数据等效源延拓与数据类型转换方法 | |
| CN107038308B (zh) | 一种基于线性内插的规则格网地形建模方法 | |
| CN106294894B (zh) | 快速分析非均匀目标电磁散射特性的有限元边界积分方法 | |
| CN111932669A (zh) | 一种基于边坡岩体特征对象的变形监测方法 | |
| WO2016074202A1 (zh) | 粒子刻蚀或沉积演化仿真方法、装置和计算机可读介质 | |
| CN103914879A (zh) | 一种在抛物线方程中由三角面元数据生成立方网格数据的方法 | |
| CN108764741B (zh) | 用于确定工厂设定区域内的生产设备布局的方法及装置 | |
| CN104778151A (zh) | 基于矩量法和抛物线方程的含腔目标电磁散射分析方法 | |
| Barakou et al. | Fractal geometry for distribution grid topologies | |
| CN112966404B (zh) | 一种三维雷电先导发展路径的生成方法 | |
| CN104778286B (zh) | 掠海飞行器电磁散射特性快速仿真方法 | |
| CN112733364A (zh) | 一种基于阻抗矩阵分块的箔条云散射快速计算方法 | |
| CN104155703A (zh) | 评价三维观测***的方法和设备 | |
| CN111339688A (zh) | 基于大数据并行算法求解火箭仿真模型时域方程的方法 | |
| CN115495981A (zh) | 一种基于遗传与粒子群算法的雷电定位优化方法 | |
| CN113868919B (zh) | 一种随钻电磁波测井3d模拟简化方法 | |
| CN107239629B (zh) | 一种岩石结构面实验室合理尺寸确定的分形维数分析方法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| 121 | Ep: the epo has been informed by wipo that ep was designated in this application |
Ref document number: 20874562 Country of ref document: EP Kind code of ref document: A1 |
|
| NENP | Non-entry into the national phase |
Ref country code: DE |
|
| 122 | Ep: pct application non-entry in european phase |
Ref document number: 20874562 Country of ref document: EP Kind code of ref document: A1 |
|
| 122 | Ep: pct application non-entry in european phase |
Ref document number: 20874562 Country of ref document: EP Kind code of ref document: A1 |
|
| 32PN | Ep: public notification in the ep bulletin as address of the adressee cannot be established |
Free format text: NOTING OF LOSS OF RIGHTS PURSUANT TO RULE 112(1) EPC (EPO FORM 1025A DATED 02/11/2022) |
|
| 122 | Ep: pct application non-entry in european phase |
Ref document number: 20874562 Country of ref document: EP Kind code of ref document: A1 |