时间域航空电磁法正演模拟

时间域航空电磁法正演模拟指导教师：谭捍东教授答辩人：印红军2012年5月19日2012年中国地质大学（北京）硕士学位论文终期答辩方法简述时间域航空电磁法（ATEM），是航空电磁方法的一个重要分支。该方法以飞机作为载体，基于不同探测目标体的电学和磁学性质差异，通过观测及研究断电后二次涡流场随时间的变化规律来探测介质的电性特征。该方法被广泛应用于矿产资源勘察、地质调查和环境监测、地下水勘察等领域。种类繁多；成本低；效率高；分辨率高；操作灵活；对屏蔽层下伏地层或构造有较好的分辨能力。A.固定翼系统B.直升机系统C.半航空系统国内外研究现状（1）1948年，加拿大试飞第一套航空电磁系统。1950年INCO航空电磁系统研制成功。60－70年代，首次在矿产开发上得到应用。70－80年代，航空电磁系统进一步优化。80年代后，几经曲折，航空电磁法仪器已经相对稳定。时间域航空电磁系统（INPU-T）,在1980年夏首次试飞，随后相继推出多种改进的型号。MEGATEMIITDEM系统INPUT系统国内外研究现状（2）我国航空电磁法研究起步稍晚。七十年代末，物化探研究所研制的双频航空电磁系统取得成功并投入使用。航空物探遥感中心于七十年代末引进的加拿大TRIDEM三频航空电磁系统也进行了大量的生产飞行，现在该系统已报废。目前，使用中的有“松鼠”型直升机吊挂型航空电磁测量系统。Y12F系统SKYTEMSURVING系统SKYTEMSURVING系统本文需要完成任务利用有限元法实现了时间域航空电磁法二维正演数值模拟。通过模型试算对二维正演算法进行精度验证。完成高低阻、组合体、倾斜体以及高低阻屏蔽层模型的数值模拟。研究加入cole-cole模型参数，对正演模拟的影响研究视电阻率求取方法边值问题基本微分方程推导由于外加源Mp和Jp很难通过离散公式表达，因此求二次场相对容易。二次场Es和Hs满足：设电磁场为按eiwt变化的谐变场，并忽略位移电流的影响，则总电场和磁场可以通过麦克斯韦方程表述0ppEiHMHEJ0sssspaEiHHEEHs不仅和Es有关，也与一次电场有关102()()()()()()()()ssssyyyysyssssyyyysyHHEEppqqiHSxxzzxzzxHHEEqqrrESxzzxxxzz12()()()()ppaxazpppayaxazSpEpEzxSEqEqExz降维处理得到：其中：220222220yyyyrkkipkkikqkkki，，，右端项是和一次场有关的量，本文一次场取无限均匀空间，对于垂直磁偶极子有：22(1)4(1)40ikrxikryzimyEikrerrimxEikrerrE边值问题对应的泛函问题：对应泛函极值问题：1,11,16112,12,162216,116,161616KKWSKKWSKKWS单元网格2222220rr(u,v){(()())(()())()221()}2uuvvuvuvJqxzxzxzzxuivdxdz,ssyyuEvH(u,v)=0u(u,v)=0vu(,)g(x,z)v(,)(,)JJxzxzgxz当边界网格取足够大时，边界节点上的电场u，v可以近似等于0.等参单元法求系数矩阵1234252627281(1)(1)(1)41(1)(1)(1)41(1)(1)(1)41(1)(1)(1)41(1)(1)21(1)(1)21(1)(1)21(1)(1)2NNNNNNNN112211112211(,)(,)(,)(,)rv((()())dx())22(,)(,)(,)(,)r((()())dx())22()(ixjxizjzTixjxizjzTixTFFFFvrdzvddvxzJFFFFuurdzudduxzJFFuvuvqdxdzquxzzx11111122001111)1(())()()222jzizjxTTTTFFddvJiwuivdxdzuNNJdduvNNJddv双二次线性插值：11()(,)11()(,)iiiixiiiiyNNNyyFxJJNNNxxFzJJ其中：xyJxx泛函可以表示为：11111111111110111112(,)(,)(,)(,)(,)(,)(,)(,)(u,v)())())22()()()22(ixjxizjzixjxizjzTTixjzizjxTTTTTFFFFFFFFrrJudduvddvJJFFFFiwquddvuNNJdduvNNJddvJS112-1-1))uSvdxdz上述公式中所示积分，采用高斯积分法求取。至此系数矩阵可以完整的表达出来。对矩阵方程的求取，目前方法比较多，本文采用波前法求解。由于篇幅所限，波前法原理在此不再赘述。1111111111211111111(,)(,)(,)(,)(u,v)())u()u2()u02(,)(,)(,)(,)(u,v)())(v2ixjxizjzixjzizjxTTTixjxizjzixTFFFFFFFFJrddqddvJJNNJddSddFFFFFJrddvqJ111111011111)u()02jzizjxTTFFFddJiwNNJddvSdd求偏导一次场计算222222232223222320(1)4(1)4[(33)(1)]4(33)4(33)4xikryikrzikrxikryikrzEimzEikrerrimyEikrerrmxHekrikrkrikrrrmxyHekrikrrrmxzHekrikrrr水平磁偶极子22223222322222232(1)4(1)40(33)4(33)4[(33)(1)]4ikrxikryzikrxikryikrzimyEikrerrimxEikrerrEmxzHekrikrrrmzyHekrikrrrmzHekrikrkrikrrr垂直磁偶极子1,11,16112,12,162216,116,161616KKWSKKWSKKWS解方程,ssyyuEvH2i-12i=u,=viiWW12()()()()ppaxazpppayaxazSpEpEzxSEqEqExz右端项令，在一个单元格中，矩阵方程如上式，将各单元的积分相加扩展到NPN*NPN全体系数矩阵，其中NPN=（2*NE-1）*(2*NZ-1)-(NE-1)*(NZ-1)，NE为横向剖分节点数，NZ为纵向剖分节点数。2i-12i=u,=viiWW时间域转换正余弦滤波算法：该算法是由汉克尔演变而来，用来解决形如（a）式形式的数值积分。1t2iwtiwtFwftedtfFwedw（a）采用傅里叶变换，从频率域转换为时间域*12*12()1[]()()sin()22()1[]()()cos()22nnnnnnnnnnnneFeetftHnFcntttteteFeetftHnFcnttttet本文采用300点滤波系数，经验算，该滤波方法在所求解数值范围内能最大限度地减小由于算法引起的误差。斜阶跃波激励响应'1()()BfBfi理想脉冲函数响应阶跃波函数响应斜阶跃函数响应''''''()1()()()oftttofdItBtBtsdsBrdrdtt阶跃波函数响应由等参单元法直接得到的脉冲激励响应，需经过变换才能得到实际应用中的斜阶跃波激励响应。斜阶跃波激励响应与发射关闭关断时间成反比，其实质是对理想阶跃关断函数响应值在时间内求平均。视电阻率求取目前方法众多，但基本都是先分别求早晚期视电阻率，再经过插值求得全时区视电阻率。本文采用，强建科在“航空瞬变电磁法的全时域视电阻率计算方法”一文中所述方法。当实际感应电动势或磁感应强度值，等于在相同参数下均匀大地正演计算值，该均匀的大地电阻率值称为“视电阻率”。因此，在某一时刻，只要找到对应采集到的感应电动势或磁感应强度值的均匀大地电阻率，便可以得到该时刻的视电阻率值。否是给定适当视电阻率初值ρ0开始δ=0.0001计算电阻率00的均匀大地上，相同观测条件下，t时刻的磁感应强度（第一次迭代0）比较磁感应强度实测值和理论值，计算其相对拟合差：()()()lBtBtBt结束计算电阻率修改量：'()()ltt0a视电阻率求取图（a）电阻率为10Ω.m的低阻异常体，埋深为40m，尺寸为40m*60m，围岩电阻率为100Ω.m。视电阻率等值线图中异常范围能大体反映出来，但是埋深偏深。图a图b图（b）埋深分别为20m和100m，电阻率分别为1000Ω.m和10Ω.m，围岩电阻率为100Ω.m，异常体尺寸都是80m*40m。算法正确性验证五层大地模型，如右图所示，飞机飞行高度保持在50m不变。选用水平共面装置，接收与发射间隔8.1m。供电电流采用幅值为10A、时长为2ms的正向和反向交替供电的阶跃波形，采集时长14ms。解析解采用时间域航空电磁一维正演算法。解析解和等参单元法解的衰减曲线和误差曲线如图所示。等参单元法误差保持在2.5%以内，说明该算法精确度能满足要求。模型分析（1）电阻率为10Ω.m（图a）和电阻率为100Ω.m（图a）的异常体，埋深为40m，尺寸为40m*60m，围岩电阻率为100Ω.m。该方法对低阻异常体反应灵敏；对高阻体反应微弱。图a图b模型分析（2）围岩电阻率为100Ω.m，异常体电阻率分别为1000Ω.m和10Ω.m，异常体间隔80m，且尺寸都为40m*180m，埋深为20m。对纵向排列组合体分辨率较横向排列组合体高。图a图b模型分析（3）围岩电阻率为100Ω.m，分别为高阻屏蔽层下伏低阻异常体和低阻屏蔽层高阻异常体模型。屏蔽层埋深20m，厚度10m，异常体埋深50m，尺寸