文章编号：1004-0609(2015)-11-3196-06
频率域线源CSAET二维反演数值计算
曹创华1, 2，童孝忠1, 2，柳建新1, 2
(1. 中南大学 地球科学与信息物理学院，长沙 410083；
2. 中南大学 有色资源与地质灾害探查湖南省重点实验室，长沙 410083)
摘 要：由于实测数据包含噪声和其固有的不完备性，决定了频率域电磁测深数据在数学计算上是一个非唯一解的、不稳定的问题，这就意味着不同的地电模型可以拟合在同一剖面上的视参数数据。面对这种问题可以通过给定加权的正则化方法使得病态逆问题变为稳定。本研究实现了频率域线源可控源音频大地电流法(CSAET法)剖面的二维反演计算。为了使反演过程稳定，采用了模型参数的平滑度约束的方法，并使用自适应正则化参数技术客观上调整了反演参数。并且在计算过程中利用了伴随方程近似计算方法，减少了计算灵敏度的耗时。数值计算表明，利用合成的数据可以得到精确的反演结果。
关键词：CSAET法；线源；二维反演；平滑限制
中图分类号：P631.3 文献标志码：A
Inversion numerical calculation of 2D CSAET generated by line current source
CAO Chuang-hua1, 2, TONG Xiao-zhong1, 2, LIU Jian-xin1, 2
(1. School of Geosciences and Info-Physics, Central South University, Changsha 410083, China;
2. Key Laboratory of Non-ferrous Resources and Geological Detection,
Ministry of Hunan province, Changsha 410083, China)
Abstract: The frequency domain electromagnetic sounding calculation is a mathematically non-unique problem with instability because the observed data contain noise and its inherent incompleteness. Thus different models can fit the data from the same cross-section well. This problem can be solved by using regularization method. This study implements the two-dimensional inversion of frequency domain line source CSAET. The method of model parameters smoothness constraints and adaptive regularization technique was used to stabilize the inversion process. Besides, the adoption of adjoin equation approximate calculation method is more time efficient in the calculation of sensitivity. The numerical calculation results show that the use of synthetic data can contribute to the accurate inversion results.
Key words: CSAET method; line current source; 2D inversion; smoothness-constrained
频率域人工源电磁测深法勘探自从推广以来已经得到广泛的应用[1]。近年来不论是正演模拟，还是反演计算，地球物理工作者都尝试了很多，也相应地取得了一定的成功[2]。尽管三维建模和反演算法现已取得很多进步[3]，野外试验也逐渐展开，由于它们计算耗时和操作麻烦，成本较大，解释技术还不成熟，还没有广泛应用。
目前，二维反演分析仍然在实际应用中扮演主要角色，但在国内外常用的还是以二维的TM极化模式进行反演，成功与失败案例平分秋色。本文作者在科研和野外工作中发现，可以尝试线源可控源音频大地电流法(以下皆称作CSAET法)剖面的二维反演。国内几乎没有单位线源二维CSAET法反演进行野外试验尝试，数值模拟和反演计算更鲜有文献论述。因此，尝试优化的线源二维CSAET反演算法显得尤为必要。
众所周知，有限单元法(FEM)和有限差分法(FDM)是对频率域大地电磁法响应进行数值模拟的主要工 具[4-5]，三维频率域电磁法正演响应也经常用积分方程法[6]。对于一维层状介质，其一般都有解析解，较为简单；但当面对二维复杂地形和复杂介质模拟时，电磁场的正演响应没有解析解，且有限差分法在模拟起伏地形条件下的频率域电磁响应时没有有限单元法精度高[7]，因此，本文作者利用有限单元法对线源二维CSAET法进行数值模拟计算。
近年来，国内外学者在电磁数据反演方面也取得了新的进展，比如TORRES-VERDIN等[8]实现了起伏地形下的Born近似层析电磁数据反演成像问题；SAKASHITA等[9]利用这一原理从钻孔电磁数据中提取了钻孔的电阻率和磁化率特征曲线；ZHDANOV等[10]发展了一种把非线性方程组线性化的拟线性策略。本文作者在其基础上，提出了一种平滑约束技术，通过反演模型参数转换，来更好的解决这个问题；通过这种方法，反演迭代中的电导率更好的得到了限制，提高了反演的准确度。
本研究中，利用平滑度约束反演方法反演了由线电流源产生的二维频域电磁响应数据。首先，采用有限单元法进行了线源二维CSAET数值模拟得到了正演响应数据。为了检查算法的有效性，在文中把解析解和数值结果做了对比分析。其次，开发了一个平滑约束最小二乘反演方法，并利用伴随方程法计算了反演过程中的灵敏度矩阵。最后，通过合成模型验证了本文算法的正确性。
1 正演计算
1.1 变分方程
本研究中假设的地电模型如图1所示， y方向为走向方向，沿y方向电阻率不变，电阻率的变化在xz平面，假设电磁波按照时间域谐波e-iωt变化，由电流源产生的电磁波在二维情况下，图示测量模式为TE模式[11]，如下式所示：
(1)
式中：Ey表示y方向的电场强度，I为电流源的电流强度，δ为狄拉克函数，ω是角频率，σ为模拟空间的电导率，μ(=4p×10-7 H/m)为磁导率。
图1 二维正演模型三维坐标系示意图
Fig. 1 3D Coordinate system for 2D forward modeling
将水平电偶极沿地层走向延伸到无限远，电偶极源变成了线源，那么对于二维地电模型来说，场也就是二维，对于外边界条件，采用无穷远边界条件，即认为电场在外边界已经衰减到非常小，接近于零，对结果的影响不大，可以忽略。式(2)变分式与边值式(1)等价，如下：
(2)
式中：Ω为计算空间，即为图1中的Site1~Site3之间的地表以下垂直于异常体的断面。
1.2 有限单元法
图2 有限元剖分模型示意图
Fig. 2 Schematic model of finite element subdivision
为了解式(2)微分方程的解，本研究利用有限单元法，把计算域剖分成能够进行模拟的四边形，如图2所示。x方向上的网格在源附近划分最密，在测区部位相对稀疏一些，其他部分划分稀疏；y方向上以地表为界，空气层划分稀疏，地层下网格剖面在探测区域划分较为密集，深部相对稀疏，本研究中采用矩形单元将整个区域剖分，矩形单元采用双二次插值，在每个单元上取8个点(4个端点和4条边的中点)。每个矩形网格内的电场场值由若干个节点值相加而成，用
表示，如下式：
(3)
式中：
表示单元节点的形函数，
表示每个节点的电场值，在局部坐标系ξ(横坐标)和η(纵坐标)下的插值函数可以表示为[12]：
(4)
利用上述方法及式(4)插值公式对式(2)进行了离散插值，由方程式(2)得到：
(5)
(6)
(7)
式中：P为与电流源相关的参数。假设每个单元的计算方程在满足边界条件下都可以得到解，然后将它们扩展矩阵相加，得到
，其中K为总体扩展矩阵，Eye代表单元各节点y方向上的电位场值的列 向量。
1.3 正演算例验证
为了验证有限单元法正演模拟的正确性，设置地电模型为均匀半空间，把利用有限单元法求出的数值解和均匀半空间条件下的解析解进行了对比。均匀半空间的解析解的表达式已经由WARD等[11]给出，如式(8)所示：
(8)
式中：K1为第二类1阶的修正贝塞尔函数；k0和k1都为电磁波在模拟介质中的波数。
假设均匀半空间的电阻率为100 Ω·m，电流源的方向沿y轴中心坐标在地表的x=0 m处。在频率100 Hz沿着x轴的有限单元数值解的场值计算结果如图3所示，其结果用红色小圆圈表示；细蓝色实线表示同样条件下由式(8)计算出的解析解，发现两者对应较好。通过此例说明正演算法是有效正确的。
图3 电阻率100 Ω·m的均匀半空间解析解和有限单元法数值解对比图
Fig. 3 Comparison of finite element method response and analytical (solid line) solution on surface of 100 Ω·m half-space
2 二维反演
2.1 平滑约束最小二乘反演
频率域电磁数据的反演问题是病态的、非线性的，可以通过正则化方法对其进行求解[13]。如下式所表示：
, (9)
式中：
为反演目标函数，
为稳定因子函数，
为正则化参数，正则化参数在求解
最小值时起权衡作用。其中和可用下两式表示：
(10)
(11)
式中：F(m)为正演响应函数，m表示模型空间向量，d表示数据空间向量，C表示模型参数的权重矩阵。
目前解此类方程有很多种数学方法[14-16]，此处利用平滑约束最小二乘反演方法求解正则化反演问题。使得式(9)线性化，如下式表示：
(12)
式中：
表示观测的视参数和正演响应数据的误差向量，
表示每次反演计算迭代过程中的模型误差向量，J为反演过程中的灵敏度矩阵或者由正演响应函数F产生的雅可比矩阵，L为拉普拉斯(二阶)平滑因子。
2.2 灵敏度计算
理论上，灵敏度的计算有如下3种方式：1) 直接计算方法；2) 灵敏度方程法；3) 伴随方程法[17]。这3种方法的计算时间分别与
、和
的正演计算相关(N是正演参数的数量，Mf为参与计算的频率的数量的数量，Mo为测量点的数量)。因为在反演计算过程中，模型参数的数量远远大于观测数据的量，所以伴随方程法计算灵敏度矩阵最为有效。
利用伴随方程法计算了灵敏度矩阵，用此方法计算了电磁数据的反演问题。为了达到反演的目的，电导率被表示为基函数的有限线性组合
[18]：
(13)
式中：
表示基函数，
表示权系数；r为搜索半径。
这样灵敏度矩阵的计算就可以由下式得出：
(14)
式中：
表示在电磁场空间
里面的辅助场场值。
2.3 正则化参数的选择
正则化参数平衡着最小的数据误差和模型粗糙度之间的关系，是评判反演过程稳定性和分辨率的重要参数，使得它们两者达到平衡。在此过程中利用自适应正则化参数计算方法，计算式如下：
(15)
式中：k表示第k次反演拟合迭代的次数，正则化参数选择与第k次迭代关系如表达式：
(16)
式中：0.01≤c≤0.5，βk为正则化因子最大值；
由式(15)自适应算法的极小值得到。
3 算例
3.1 模型一
如图4(a)所示，设计的模型在x和z轴组成的平面上，设计剖面的长度为6000 m，深度为3000 m，其中x轴上的坐标范围为(-3000 m，3000 m)，z轴上的坐标范围为(-3000 m，0 m)，围岩电阻率设置为500 Ω·m；异常体设置2个，位置在x轴(-1000 m，-500 m)、z轴(-1000 m，-500 m)的正方形异常设置为1000 Ω·m，位置在x轴(500 m，1000 m)、z轴(-1000 m，-500 m)的正方形异常设置为250 Ω·m。
反演的初始网格由52列和26行组成，共计1352个剖分单元，选择的背景场电阻率值为500 Ω·m，测量点为30个，点距200 m，选择的频点为分别为0.125、0.25、0.5、1、2、4、8、16、32、64和128 Hz，共计11个频点。图4(b)所示即为反演结果，从图4(b)中可以看出：反演后最大的电阻率变为660 Ω·m，最小电阻率为320 Ω·m，从反演结果的模型电阻率等值线断面可以清楚的识别低阻区域和高阻区域。这从侧面说明平滑度约束反演算法处理CSAET频率域电磁数据时能很好地反演简单的组合模型，但高阻异常与初始模型相比，电阻率变低；低阻异常与初始模型相比，电阻率变高。
图4 模型一正演初始模型和反演结果
Fig. 4 Performance test of 2D smoothness-constrained inversion for simple model
3.2 模型二
图5(a)所示为SASAKI提出的电磁法经典模型[19]。剖面线总长为17500 m，深度4000 m，其背景电阻率值为50 Ω·m，剖面深部镶嵌两个5 Ω·m的低阻异常，在地表附近设置了一个100 Ω·m的相对高阻异常和一个10 Ω·m的相对中等电阻率异常。此剖面设置为35个测点，点距为500 m，选择的频点仍分别为0.125、0.25、0.5、1、2、4、8、16、32、64和128 Hz共11个频点。
此模型的正演响应仍然利用上文中提到的有限单元法，其初始网格由40列、和28行组成，共计1120个剖分单元，利用本研究中提到的方法进行反演。图5(b)所示为此模型的反演结果，可见反演结果分层较为明显，电阻率范围与模型一类似，具有高于背景电阻率的异常反演后电阻率变低而低于背景电阻率的异常反演后电阻率值变高的现象。
通过上述两个模型的计算，表明二维平滑度约束方法可以有效地对CSAET频率域电磁数据进行反演。
图5 模型二正演初始模型(a)和反演结果(b)
Fig. 5 Performance test of 2D smoothness-constrained inversion for model as same as one tested by SASAKI[19]
4 结论
1) 推导了垂直频率域线源频率域CSAET法二维电场微分方程，并推导了其变分方程，利用有限单元法进行程序编制，实现了其正演响应的求解。
2) 编制了二维平滑度约束方法反演算法，实现了简单模型和SASAKI复杂模型的反演计算，对设置地层和异常体反演结果较好。
3) 反演过程中采用自适应正则化参数技术实现了自动调节反演过程中的正则化因子。
4) 从理论上分析了伴随方程近似法计算灵敏度的耗时较小，并成功的进行了反演计算。
REFERENCES
[1] 刘春明, 佟铁钢, 何继善. 多种电磁法在某金矿的野外勘探应用[J]. 中国有色金属学报, 2013, 23(9): 2422-2429.
LIU Chun-ming, TONG Tie-gang, HE Ji-shan. Exploration of various electromagnetic methods in some gold mine[J]. The Chinese Journal of Nonferrous Metals, 2013, 23(9): 2422-2429.
[2] UNSWORTH M J, TRAVIS B J, CHAVE A D. Electromagnetic induction by a finite electric dipole source over a 2-D earth[J]. Geophysics, 1993, 58, 198-214.
[3] COMMER M, NEWMAN G A. New advances in three-dimensional controlled-source electromagnetic inversion[J]. Geophysical Journal International, 2008, 172: 513-535.
[4] COGGON J. Electromagnetic and electrical modeling by finite element method[J]. Geophysics, 1971, 36: 132-155.
[5] 阎 述, 陈明生. 线源频率电磁测深二维正演(一)[J]. 煤田地质与勘探, 1999, 27(5): 60-62.
YAN Shu, CHEN Ming-sheng. Two-dimensional forward modeling of line source frequency electromagnetic sounding (Ⅰ)[J]. Coal Geology & Exploration, 1999, 27(5): 60-62.
[6] 李帝铨, 谢 维, 程党性. E-Ex 广域电磁法三维数值模拟[J]. 中国有色金属学报, 2013, 23(9): 2459-2470.
LI Di-quan, XIE Wei, CHEN Dang-xing. Three-dimensional modeling for E-Ex wide field electromagnetic methods[J]. The Chinese Journal of Nonferrous Metals, 2013, 23(9): 2459-2470.
[7] 柳建新, 童孝忠, 郭荣文, 李爱勇, 杨 生. 大地电磁测深法勘探—资料处理、反演与解释[M]. 北京: 科学出版社, 2012: 58-98.
LIU Jian-xin, TONG Xiao-zhong, GUO Rong-wen, LI Ai-yong, YANG Sheng. The magnetotelluric sounding prospecting—Data processing, inversion and interpretation[M]. Beijing: Science Press, 2012: 58-98.
[8] TORRES-VERDIN C, HABASHY T M. Rapid 2.5-D forward modeling and inversion via a new nonlinear scattering approximation[J]. Radio Science, 1994, 29: 1051-1079.
[9] SAKASHITA S, FUKUOKA K. An imaging technique for magnetic susceptibility and resistivity by electromagnetic tomography-Numerical experiments: Buturi-Tansa[J]. Geophysics Exploration, 1998, 51: 493-506.
[10] ZHDANOV M S, FANG S. Three-dimensional quasi-linear electromagnetic inversion[J]. Radio Science, 2012, 31: 741-754.
[11] WARD S H, HOHMANN G W, NABIGHIAN M N. ElectromagneticMethods in applied geophysics[C]// ElectromagneticTheoryforGeophysicalApplications. Tulsa, USA: SocietyforExploration in Geophysics, 1988: 131-312.
[12] 徐世浙. 地球物理中的有限单元法[M]. 北京: 科学出版社, 1994: 5-28.
XU Shi-zhe. Finite element method in geophysics[M]. Beijing: Science Press, 1994: 5-28.
[13] TIKHONOV A N, ARSENIN V Y. Solution of Ill-Posed Problems[M]. New York: Wiley, 1977: 3-22.
[14] ZHDANOV M S. Geophysical inverse theory and regularization problems[M]. Amsterdam: Elsevier, 2002: 5-29.
[15] MITSUHATA Y, UCHIDA T, AMANO H. 2.5-D inversion of frequency-domain electromagnetic data generated by a grounded-wire source[J]. Geophysics, 2002, 74: 1753-1768.
[16] CANDANSAYAR M E. Two-dimensional inversion of magnetotelluric data with consecutive use of conjugate gradient and least-squares solution with singular value decomposition algorithms[J]. Geophysical Prospecting, 2008, 56: 141-157.
[17] MCGILLIVRAY P R, OLDENBURG D W. Methods for calculating Fréchet derivatives and sensitivities for the non-linear inverse problem: A comparative study[J]. Geophys Prospect, 1990, 38: 499-524.
[18] FARQUHARSON C G, OLDENBURG D W. Approximate sensitivities for the electromagnetic inverse problem[J]. Geophysical Journal International, 1996, 126: 235-252.
[19] SASAKI Y. Two-dimensional joint inversion of magnetotelluric and dipole-dipole resistivity data[J]. Geophysics, 1989, 54: 254-262.
(编辑 何学锋)
基金项目：国家高技术研究发展计划资助项目(2014AA06A615)；国家自然科学基金资助项目(41174103)；中南大学中央高校基本科研业务费专项资金博士生自由探索项目(2013zzts054)；国家教育部博士点基金资助项目(20110162120064)
收稿日期：2015-03-10；修订日期：2015-05-15
通信作者：童孝忠，讲师，博士；电话：13469439855；E-mail：csumaysnow@163.com
