摘要
直升机瞬变电磁方法是为了适应地表单元复杂多变、起伏剧烈和地质结构复杂等特殊环境而产生的一种探测技术。直升机电磁法自诞生以来，已被大量用于金属矿勘查、煤炭、油气勘查、海水入侵、水文地质调查和环境监测等领域，并以快速、经济而著称。我国矿产资源分布形态多样，多呈透镜状、扁豆状、脉状等分布，为了减小层状大地模型近似实际矿体模型而产生的误差，应采用三维异常体近似透镜状等矿体分布。已有的三维有限差分方法计算直升机瞬变电磁响应存在计算时间长、无吸收边界条件等问题。
本文基于以上问题，在国家863重大科技攻关项目“吊舱式时间域直升机航空电磁勘查系统”和国家自然基金项目“基于全波形航空时域电磁探测的三维地质体识别关键技术研究”的支撑下，做了以下几方面的研究。在小宗量近似条件下，推导出了回线源在空中发射、地下接收的初始场计算式，并且实现了初始场的并行计算；以三次样条插值为基础，研究了非扭结边界条件下的双三次样条插值并行计算方法；将MAXWELL的旋度方程差分离散，推导出了电磁场在计算区域内三个分量上的计算式，采用三维Mur吸收边界条件，推导出了截断边界上的电场分量计算式，并实现了迭代场的并行计算，缩短了计算时间；分析了阶跃波直升机瞬变电磁响应特征，并研究了斜阶跃波的直升机瞬变电磁响应计算方法。通过上述研究提高了三维有限差分方法计算直升机瞬变电磁响应的速度，在相同计算区域内，提高了响应的计算精度，为三维反演提供了基础。
With the rapid development of China's economy, the demand for mineral resourcesalso increases dramatically. The contradiction between supply and demand has becomeincreasingly intensified, where the demand of the bulk minerals such as iron, manganese,chromium, copper, aluminum, potassium salt is particularly prominent. Helicoptertransient electromagnetic method is an exploration technology which is developed becauseof the special environment, for example, forest-swamp landscape area and desert-Gobilandscape area. In China, the shape of the mineral resource's distribution of is diverse.Most of them present as lenticular, lentiform and vein type. However, the actual situationis just the opposite. The theoretical investigation of helicopter transient electromagneticmethod is mainly based on the uniform half-space model or layered earth model. On thebasis of that, one-dimensional or two-dimensional approximate inversion surely has awide gap comparing with the actual structure of the ore body. In order to reduce the error,we use a three-dimensional anomalous body model to approximate lenticular, lentiform etal mineral distribution. The following problems exist in3D finite difference numericalcalculation method. Firstly, there are no absorbing boundary conditions in the truncatedborders of computational domain. Artificially, we let the value of truncated borders beequal to zero. Do like this, if computational domain is not large enough, the field will bedivergence in the late stage. Secondly, the calculation time is long. So it is unable to meetthe rapid demand of the3D electromagnetic inversion.
Based on the problems above, the thesis is in the support of863major tackle keyproblems in science and technology project "the bird time domain helicopter-borneelectromagnetic exploration system", and national natural science foundation project "research about the3D geological identify key technology based on the full waveformaviation time-domain electromagnetic detection". In the basis of serial computing of3Dfinite-difference transient electromagnetic response which doesn't have the absorbingboundary conditions, we use heterogeneous platforms based on CPU and GPU to researchthe parallel computing of3D finite difference helicopter transient electromagneticresponse with Mur absorbing boundary conditions. Major researches and innovations areas follows:
1．According to the expressing which the loop is on the surface and receiving isunderground, we derive time domain initial field formula of step incentives which the loopis in the air and receiving is underground. By means of the small amount of approximationmethod, making the formula only have one Bessel function. Applying the MAXWELLfirst and second equation, through differential discrete method, I derive pulse initialelectromagnetic field which is needed in the iteration field computation. At last, realizeinitial field parallel computation. When the grid is101×101×50, acceleration ratio is21.04.When the grid is201×201×150, acceleration ratio is13.2. The denser the grid is, the moretime can be saved by parallel calculation.
2． Research the calculation method of bicubic spline interpolation under thenot-a-knot boundary condition: decompose the bicubic spline interpolation into two cubicspline interpolations, interpolating along the x direction first, and then y direction. At last,the row or column of two-dimensional array to be interpolated is assigned to Thread in theGPU to realize parallel calculation of interpolation. Ultimately, the parallel computingresult and interpolation function in the matlab are compared; the error is controlledbetween10-4%~10-8%.
3．Using3D finite differential time domain to discrete curl equations in theMAXWELL equation，I derived x, y, z components of electric and magnetic field in thecomputational domain. Adopting3D Mur absorbing boundary conditions, making use offirst order and second order approximation derive electric field formula on the truncatedboundary. Solving the problem, when the rock resistivity is larger, using no absorbingboundary conditions to calculate3D electromagnetic response will produce the divergence in the late time. Effectively extend the electromagnetic response calculation time andimprove accuracy. Applying upward-continuation theory solves the problem aboutmagnetic field component calculation in the air. At last, achieve parallel computing ofiteration field, acceleration ratio is12.03and7.32respectively, and supply the foundationfor3D inversion of anomaly.
4. Transmitter-receiver loop applies the central loop device. I studies helicoptertransient electromagnetic response characteristics of layered earth model, plate-like bodywith or without covering layer model and3D abnormal body model. Supply forwardsample data for inverse explain of3D target. According to the abnormal morphology ofprofile curve of plate-like body in the homogeneous half-space model, the tendency of theplate-like can be roughly determined. Single peak is a horizontal plate-like body,equivalent double peak is an upright plate-like body, and the non-equivalent double peakis an inclined plate-like body. The meters changing which the aircraft deviates from thenormal flight altitude will cause false anomaly, and therefore it need to be corrected. Bythe thinking of distinguished time window, I research the detection capability of helicoptertransient electromagnetic detection system. Only changing the conductivity of theanomalous body, the other parameters are constant, then the greater the conductance ratiobetween the anomaly and surrounding rock, the more easily detected the anomaly. Whenonly change the depth of anomaly, the shallower the anomaly is, the more easily detectedthe anomaly.
5. Combining the G-S inverse Laplace transform with F.N.Kong241's Hankeltransform calculates helicopter transient electromagnetic response of the layered earthwhen transmitted current is a ramp function and is completely switched off. By theconvolution between the first derivative of trapezoidal wave and the3D anomaly responsewhen transmitted current is a step function, I calculate3D anomaly helicopter transientelectromagnetic response when the transmitted current is a trapezoidal form.
The innovation of the paper are as follows：
1. Put forward the bicubic spline interpolation computing method under not-a-knotboundary condition. Solve the problem that we can't calculate the first and second order derivative of magnetic field value at two end points of the computing area and solveinterpolation calculation problem in the upward-continuation, to lay the groundwork forthe GPU-based iterative field of parallel computing.
2. Put forward applying the truncated boundary’s electric field formula to computeiteration field of3D finite difference helicopter transient electromagnetic response.Solving the problem, when the rock resistivity is larger, using no absorbing boundaryconditions to calculate3D electromagnetic response will produce the divergence problemin the late time. Effectively improve the accuracy of numerical computation ofelectromagnetic lately.
3. By analyzing numerical calculation features of initial field, electromagnetic fieldwithin the region, truncated boundary’s electric field and magnetic field in the air, putforward applying parallel computing method based on GPU to initial field and iterationfield of electromagnetic response. Complete the parallel calculation of helicopter transientelectromagnetic response of3D finite differential method based on the GPU, savingcomputing time.
引文
[1]余丽秀，孙亚光.矿产资源综合利用与地质调查关系探讨[J].矿产保护与利用，2008（5）：10-13.
[2]高乐，王琨，李红中，等.地球物理方法在矿产资源探查中的运用[J].矿物学报，2011，778.
[3]张潇，崔斌，岑况，等.利用“两个市场”开发“两种资源”推动我国矿业发展[J].矿产综合利用，2006（3）：47-49.
[4]陈从喜.试论矿产资源综合利用与地质找矿[J].国土资源情报，2009（7）：13-19.
[5]杨沈生.我国矿产资源现状与应对措施[N].科技创新导报，2010（12）：70.
[6]滕吉文.中国地球物理学研究面临的机遇、发展空间和时代的挑战[J].地球物理学进展，2007，22（4）：1101-1112.
[7]杨少平，焦保权，孙忠军，等.森林沼泽景观区区域化探异常追踪方法技术[J].物探与化探，2008，32（5）：480-487.
[8]金浚，陈伟民，李占龙，等.中国东北地区森林沼泽景观地球化学特征与勘查方法[J].矿产勘查，2010，1（1）：60-67.
[9]陈红恩，倪新辉，韩爱英，等.高精度三维地震勘探技术在青藏高原沙漠戈壁区煤田的应用[J].中国煤田地质，2007，19:117-119.
[10]沈骥千.沙漠戈壁地区地震资料层析静校正研究与应用[J].中国煤炭地质，2012，24（6）：57-59.
[11]Richard S. Smith, A Peter Annan, Patrick D McGowan.A comparison of data from airborne,semi-airborne, and ground electromagnetic systems[J].Geophysics,2001,66(5):1379–1385.
[12]薛国强，李貅，底青云.瞬变电磁法正反演问题研究进展[J].地球物理学进展，2008，23（4）:1165-1172.
[13]欧阳立胜.地面瞬变电磁法(TEM)研究进展[J].四川地震，2011，4：40-45.
[14]杨光，陈玉玖，姜志海.小回线源瞬变电磁法在煤矿积水[J].工程地球物理虚报，2011，8（4）：309-402.
[15]杨云见，何展翔，赵晓明.接地长导线源瞬变电磁法全区视电阻率定义探讨[J].物探装备，2010，20（2）：117-120.
[16]朱凯光，林君，韩悦慧，等.基于神经网络的时间域直升机电磁数据电导率深度成像[J].地球物理学报，2010，53（3）：743-750.
[17]毛立峰，陈斌，吕东伟.测高数据不准时的直升机航空瞬变电磁一维反演方法理论研究[J].物探与化探，2011，35(3)：402-405.
[18]雷栋，胡祥云，张素芳.航空电磁法的发展现状[J].地质找矿论丛，2006，21(1)：40-53.
[19]张昌达.航空时间域电磁法测量系统：回顾与前瞻[J].工程地球物理学报，2006，3(4)：265-273.
[20]Fountain D.60Years of airborne EM-focus on the last decade,2008[C]AEM2008-5th InternationalConference on Airborne Electromagnetics, Finland.
[21]David Fountain. Airborne electromagnetic system-50years of development [J]. ExplorationGeophysics,1998,29：1-11.
[22]朴化荣.电磁测深法原理[M].北京：地址出版社，1990.
[23]J H Knight, A P Raiche.Tansient electromagnetic calculations using the Gaver-Stehfest inverseLaplace transform method[J].Geophysics,1982,47(1),47-50.
[24]罗延钟，昌彦君.G-S变换的快速算法[J].地球物理学报，2000，43(5)：684-690.
[25]D. Guptasarma.Computation of the time-domain response of a polarizable ground[J]. Geophysics,1982,47(11):1574-1576.
[26]昌彦君，张桂青.电磁场从频率域转换到时间域的几种算法[J].比较物探化探计算技术，1995，17(3)：25-29.
[27]张榴晨，徐松.有限元法在电磁计算中的应用[M].北京：中国铁道出版社出版发行，1996.
[28]曾余庚，徐国华，宋国乡.电磁场有限单元法[M].北京：科学出版社，1982.
[29]葛德彪，闫玉波.电磁波时域有限差分方法（第三版）[M].西安：西安电子科技大学出版社，2011.
[30]Wannamaker Philip E, Hohmann G W, SanFilipo W A. Electromagnetic modeling of three-dimensional bodies in layered earths using integral equations[J]. Geophysics,1984,49(1):60~74.
[31]Sanfilipo William, Hohmann G W. Integral equation solution for the transient electromagneticresponse of a three-dimensional body in a conductive half-space[J]. Geophysics,1985,50(5):798~809.
[32]Wannamaker P E. Advances in three-dimensional magneto telluric modeling using integralequations [J]. Geophysics,1991,56(11):1716~1728.
[33]Dmitriev V I, Nesmeyanova N I. Integral equation method in three-dimensional problem s of low-frequency electrodynamics[J]. Computation al Mathematics and Modeling,1992,3(3):313~317.
[34]Juhl T lb ll R, B ie Christensen N. Robust1D inversion and analysis of helicopterelectromagnetic (HEM) data[J].GEOPHYSICS,2006,71(2),53–62.
[35]Christiansen, Anders Vest,Auken, Esben,et al. Quantification of modeling errors in airborne TEMcaused by inaccurate system description [J]. GEOPHYSICS,2011,76(1),43-52.
[36]Burschil Thomas, Wiederhold Helga, Auken Esben. Seismic results as a-priori knowledge forairborne TEM data inversion-A case study[J]. APPLIED GEOPHYSICS,2012,80,121-128.
[37]Smiarowski Adam, Macnae James, Bailey R C. Predictive filter calculation of primary fields in afixed-wing time-domain AEM system[J].Geophysics,2010,75(4),97-106.
[38]Sattel Daniel. Inverting airborne electromagnetic (AEM) data with Zohdy's method[J].GEOPHYSICS,2005,70(4):77-85.
[39]Reid J, Pfaffling A, Vrbancich J. Airborne electromagnetic footprints in1D earths[J].GEOPHYSICS,2006,71(2):63–72.
[40]Haoping Huang, Douglas C Fraser. Inversion of helicopter electromagnetic data to a magneticconductive layered earth[J].GEOPHYSICS,2003,68(4),1211–1223.
[41]P H Nelson, D B Morris. Theoretical response of a time-domain, airborne, electromagneticsystem[J].Geophysics,1969,34(5):729-738.
[42]Richard S Smith, Pierre B Keating. The usefulness of multicomponent, time-domain airborneelectromagnetic measurements[J].Geophysics,1996,61(1):74–81.
[43] YIN C C, SMITH R S, HODGES G, et al. Modeling Results of on-and off-time B and dB/dt forTime-Domain Airborne EM Systems:70th EAGE Conference&Exhibition-Rome,Italy,9-12June2008[C]. Roma: EAGE,2008.
[44]黄皓平.航空电磁法应用的新领域—盐渍土及次生盐渍化调查[J].外地质勘探技术，1992，2(2)：16-18.
[45]丁志强，程志平，董浩，等.航空电磁法筛选金属矿异常技术研究[J].地质与勘探，2012，48(3)：601-610.
[46]王卫平，王越胜.航空电磁法在黄河口地区寻找浅层淡水的地质效果[J].物探与化探，1999，23(2)：115-121.
[47]周凤桐，陈本池.航空电磁法的应用现状与前景[J].国外地质勘探技术，1997，2：1-5.
[48]刘慧勤.山东省昌邑—平度地区航空电磁法测量及应用效果[J].山东地质，1998，14(1)：53-58.
[49]罗延钟，张胜业，王卫平.时间域航空电磁法一维正演研究[J].地球物理学报，2003，46(5)：719-724.
[50]黄皓平，王堆中.时间域航空电磁数据的反演[J].地球物理学报，1990，33(1)：87-97.
[51]嵇艳鞠，林君，关珊珊等.直升机航空TEM中心回线线圈姿态校正的理论研究[J].地球物理学报，2010，53(1)：171-176.
[52]嵇艳鞠，栾卉，李肃义.全波形时间域航空电磁探测分辨率[J].吉林大学学报(地球科学版)，2011，4l(3)：885-891.
[53]闵刚，王绪本，毛立峰，等.磁偶极子源航空瞬变电磁对飞行高度的响应特征[J].物探与化探，2012，36(3)：1-4.
[54]许洋铖，林君，李肃义，等.全波形时间域航空电磁响应三维有限差分数值计算[J].地球物理学报，2012，55(6)：2015-2114.
[55]朱凯光，林君，韩悦慧，等.基于神经网络的时间域直升机电磁数据电导率深度成像[J].地球物理学报，2010，53(3)：743-750.
[56]强建科，罗延钟，汤井田，等.航空瞬变电磁法的全时域视电阻率计算方法[J].地球物理学进展，2010，25(5)，1657-1661.
[57]毛立峰，陈斌，吕东伟.测高数据不准时的直升机航空瞬变电磁一维反演方法理论研究[J].物探与化探，2011，35(3)：402-405.
[58]Raiche A P. An Integral equation approach to three-dimensional modeling[J]. Geophys,1974,36:363~376.
[59]Best M E,Duncan P, Jacobs F J,et al.Numerical modeling of the electromagnetic response ofthree-dimensional conductors in a layered earth[J].Geophysics,1985,50:665-676.
[60]Jopie I Adhidjaja,Gerald W Hohmann. A finite difference algorithm for the transientelectromagnetic response of a three-dimensional body[J].Geophys,1989,98:233-242.
[61]Zivanovic S S,Yee K S, Mei K K. A subgridding method for the time-domain finite-differencemethod to solve MAXWELL's equation[C].IEEE trans.Microwave Theory and Tech,1991,471-479.
[62]Tsili Wang, Gerald W Hohmanni. A finite-difference, time-domain solution for three-dimensionalelectromagnetic modeling[J].GEOPHYSICS,1993,58:797-809.
[63]邓正栋，关洪军，聂永平，等.稳定地电场三维有限差分正演模拟[J].石油物探，2001，40(1)：107-114.
[64]王志刚，何展翔，魏文博.井中垂直双极源体积分方程法三维模拟研究地[J].球物理学进展2007，2(6)：1802-1808.
[65]蒋大青，付志红，侯兴哲，等.基于MAXWELL3D瞬变电磁法三维正演研究[J].电测与仪表，2012，49(588)：29-32.
[66]印红军.时间域航空电磁法二维正演研究[D].北京：中国地质大学硕士论文，2012.
[67]G A Newman, D L Alumbaugh. Three-dimensional massively parallel electromagneticinversion-I[J]. Theory Geophys,1997,128：345-354.
[68]Fabio I Zyserman, Juan E Santos. Parallel finite element algorithm with domain decompositionforthree-dimensional magnetotelluric modeling[J]. Applied Geophysics,2000,44:337–351.
[69]匡斌，李心友，王华忠，等.三维有限差分深度偏移并行算法的设计和实现[J].同济大学学报，2000，28（2）：183：188.
[70]Michael Commer,Gregory Newman.A parallel finite-difference approach for3D transientelectromagnetic modeling with galvanic sources[J]. Geophysics,2004,11:1192-1202.
[71]Tan Han dong,ong Tu, Lin Chang hong.The Parallel3D magnetotlluric forward modelingalgorithm[J].APPLIED GEOPHYSICS,2006,3(4):197-202.
[72]终拓，谭捍东，林昌洪，等.基于MPI的频率域航空电磁法二维有限元数值模拟并行算法研究[C].第十届中国国际地球电磁学术讨论会，2011，73-76.
[73]张帆.基于MPI和GPU直流电法和大地电磁法的三维正演并行算法[D].北京：中国地质大学，2011.
[74]史维，严良俊，谢兴兵.MPI环境的并行算法在瞬变电磁[J].工程地球物理学报，2012，9（1）：75-80.
[75]牛之琏.时间域电磁法原理[M].长沙：中南大学出版社，2007.
[76] Misac N.Nabighian.勘查地球物理-电磁法[M].北京：地质出版社，1992.
[77]李貅.瞬变电磁测深的理论与应用[M].西安：陕西科学技术出版社，2002.
[78]朴化荣.电磁测深法原理[M].北京：地址出版社，1990.
[79]蒋邦远.实用近区磁源瞬变电磁法勘探[M].北京：地址出版社，1998.
[80]傅良魁.电法勘探教程[M].北京：地址出版社，1983.
[81]Report on a combined TURAIR electromagnetic and magnetic survey,1975, Mining SpectersOffice Whitehouse, Y.T.
[82]Toru Mogi, Ken’ichirou Kusunoki, Hideshi Kaieda, et al. Grounded electrical-source airbornetransient electromagnetic (GREATEM) survey of Mount Bandai,north-eastern Japan[J].Exploration Geophysics,2009,40:1-7.
[83]Peter Elliott. The Principles and Practice of FLAIRTEM[J]. Exploration Geophysics,1998,29:58-60.
[84]陈国良.并行计算-结构·算法·编程[M].北京：高等教育出版社，2003.
[85]张武生，薛巍，李建江，等.MPI并行程序设计实例教程[M].北京：清华大学出版社，2009.
[86]张舒，褚艳利，赵开勇，等.GPU高性能运算之CUDA[M].北京：中国水利水电出版社，2009.
[87]翟艳堂，涂强，郎显宇.基于CUDA的蛋白质翻译后修饰鉴定MS-Alignment算法加速研究[J].计算机应用研究，2010，27(9)：3409-3414.
[88]王占广，罗忠涛，李军.基于CUDA的机载MIMO雷达杂波建模[J].火控雷达技术，2011，40(4)：35-39.
[89]王英俊，王启富，王钢.CUDA架构下的三维弹性静力学边界元并行计算[J].计算机辅助设计与图形学学报，2012，24(1)：112-119.
[90]夏春兰，石丹，刘东权.基于CUDA的超声B模式成像[J].计算机应用研究，2011，28(6)：2011-2015.
[91]吴连贵，易瑜，李肯立.基于CUDA的地震数据相干体并行算法[J].计算机应用，2009，29(3)：912-914.
[92]黄兴，宋建新.基于GPU的视频转码技术研究[J].视技术，2012，36(1)：26-29.
[93]Zhongliangv译. NVIDIA_Fermi架构_白皮书_中文详细版[EB/OL].[2010-09-28].http://wenku.baidu.com/view/6b4e8ce1524de518964b7d9d.html.
[94]Jason sanders,Edward kandrot.CUDA范例精解-通用GPU编程（影印版）[M].北京：清华大学出版社，2010.
[95]NVIDIA Corporation. NVIDIA CUDA C Programming Guide Version4.2[EB/OL].[2012-04-05].http://developer.nvidia.com/cuda-downloads.
[96]风辰.风辰的CUDA入门教程[EB/OL].[2010-07-24]. http://wenku.baidu.com/view/bdcac5260722192e4536f6da.html.
[97]李庆扬，王能超，易大义.数值分析[M].北京：清华大学出版社，2001.
[98]王新民，术洪亮.计算方法[M].北京：高等教育出版社，2005.
[99]NVIDIA Corporation.CUDA CUFFT Library[EB/OL].[2010-02]. http://wenku.baidu.com/view/5814587c31b765ce0508140f.html.
[100]Raiche A P.The Effect of Ramp Function Turn-off on the TEM Response of Layered Earth[J].Exploration Geophysics,1984,15:37-41.