0 引言
积分方程法在计算总场时形成的矩阵是稠密矩阵，对于小规模三维数值模拟，积分方程法占据着精度高和速度快的优势；但随着异常区域规模的增大，积分方程法对内存的需求和计算时间都会急剧增加，限制了该方法在三维反演中的应用。为克服积分方程法这个缺点，很多学者将目光放在了总场的近似求解方法上。
Born近似方法[1]最早被引入地球物理数值模拟及反演研究当中。它假设散射场为零、总场为一次入射场，避开了求解方程组，大大缩短了正演模拟的时间。但Born近似方法只能解决弱反射问题，在电磁法中只适用于电导率差异较小、异常区域较小的低频电磁场计算；电导率差异增大或者异常区域增大都会增加电磁散射强度，此时利用Born近似方法计算的结果已经失去了准确性[2]。在Born近似方法的基础上，Torres-Verdin等[3-4]提出了局部非线性散射近似方法，扩展了Born近似方法的适用范围且计算速度较快；但该方法仅适用于2维或2.5维情况，且异常电导率在背景电导率50倍以内[3]。Zhdanov等[5-6]首先提出拟线性近似方法，包括对角拟线性近似法与标量拟线性近似法。拟线性近似方法通过引入反射系数λ，并假设λ在异常区域内部是慢变的、异常体的散射场等于λ与入射场的乘积，将非线性问题转化为求解λ的拟线性问题。随后，Zhdanov等[7]又提出标量反射系数形式的拟解析近似方法，进一步简化了总场计算方程。受Torres-Verdin等[3-4]思路的影响，Zhdanov等[8]提出了适用于多源电磁模拟的局部拟线性近似方法。该方法在拟线性近似方法的基础上，假设背景场在异常区域内缓慢变化，继而反射系数的求取与背景场无关，只需要计算一次反射系数。Ueda等[9]提出了多重网格拟线性近似方法，通过计算大网格单元的反射系数再插值出小网格单元的反射系数来加快反射系数矩阵的生成速度。孙建国[10]对积分方程中的线性近似及拟线性近似方法进行了综合性评述。
Zhdanov等提出的拟线性系列近似方法相对于Born近似等其他近似方法在计算电磁场的精度上有了较大的提高，适用于电导率差异变化更大的模型，但这一系列的拟线性近似方法各自的适用范围至今仍然比较模糊，导致实际应用的计算速度和计算精度还没有一个清晰的概念。本文首先系统地推导了对角拟线性近似、标量拟线性近似、拟解析近似以及局部拟线性近似的计算公式并实现了其对应算法。接着，开展了不同近似方法的三维电磁场数值模拟。通过对模拟结果的对比分析，总结了各种近似方法适用的电导率变化范围，特别分析了对角拟线性近似方法精度高和适用性强的原因。最后，对比了各种近似方法的计算速度。
1 拟线性近似理论均匀大地中的三维异常体如图 1所示。利用积分方程法对其进行模拟，将异常体剖分为M个均匀小单元，则积分方程离散，转化为M个未知量的线性代数方程组。通过推导Maxwell方程组，我们直接给出体积分方程响应公式：
|
| 图 1 三维异常体模型 Figure 1 3D abnormal body |
|
|
(1) 式中：E(rj)为总电场；Eb(rj)为一次电场；Ω为异常区域；j为当前单元编号；rj为当前单元；r为任意单元；
(2) 式中，Ea为二次电场。
拟线性近似理论[5-6]基于这样的假设：异常区域内的异常场与背景场呈线性关系，通过反射系数λ可将异常场表示为
(3) 将方程(3) 代入方程(2) 会得到拟线性近似解EQLa：
(4) 式中, 
λ可以通过求解下式极小值问题得到：
(5) 式中，L2(D)为二范数；
下面给出具体的求解方法。
为了便于数值计算，将式(4) 写成
(6) 式中：λβγ为反射系数张量的笛卡尔坐标系分量；δαβ和δβγ为delta函数；Eαa(rj)、Eγb(rj)为电场矢量的笛卡尔坐标系分量；
在异常区域内反射系数是缓慢变化的[8]，因此，可近似地将积分中的反射系数提到积分之外，则式(6) 写为
(7) 然后可以得到
(8) 式中，EαB(rj)为Born近似解[1]。
通过最小二乘方法可确定满足式(5) 的λ：
(9) 式中：J为剖分块数；*代表复数的共轭；ΔEαβγ为
(10) 将式(9) 中的χ(λ)泰勒展开并取其一阶变化函数：
(11) 对于任意的δλ，令δχ(λ, δλ)=0，则
(12) 定义对角形式的反射系数为
(13) 二次电场可通过下式求解：
(14-1) 
(14-2) 
(14-3) 在异常区域内，反射系数是缓慢变化的，近似地，可将积分中的反射系数提到积分之外并整理可得
(15) 将方程组(15) 写成矩阵相乘的形式：
(16) 其中：
进一步对反射系数进行标量化近似处理。令式(13) 中对角线上的反射系数与标量反射系数相等：
(17) 则标量反射系数可由下式解出：
(18) 如果反射系数λ采用标量形式，那么式(4) 可写为
(19) 对式(19) 两端乘以Eb*(rj)，有
(20) 如果Eb(rj)Eb*(rj)≠0，将式(20) 两端除以Eb(rj)Eb*(rj)，有
(21) 式中，
(22) 如果Eb(rj)Eb*(rj)=0，取λ(rj)=－1。
将式(21) 代入式(4) 得到
(23) 那么二次电场的求解可由下式给出：
(24) 式中，EQAa(rj)为拟解析近似解。
1.4 局部拟线性(LQL)近似方法由于异常区域内部背景场是缓慢变化的，那么
(25) 由此，反射系数可通过下式取极小得到
(26) 因此，反射系数的计算与背景场无关，在多源电磁法中只需要计算一次反射系数。
2 正演模拟给出如图 2所示的三维模型：均匀半空间介质电阻率为1 000 Ω·m，异常体为200 m×200 m×200 m的正方体，埋深100 m。发射源坐标为(-1 000, 0, 0)、长度为200 m，发射电流为1 A。将异常体均匀剖分为8个大小相等、边长为100 m的正方体单元，通过DQL、QA、QL、LQL近似方法及积分方程(IE)法进行模拟，各种方法的剖分网格相同。
|
| ρb为背景电阻率。 图 2 正演模型 Figure 2 Forward model |
|
|
对于高阻异常体，拟线性系列近似方法都能给出精确的结果，所以这里只讨论低阻情况。给定异常体电阻率为100、10、1 Ω·m三种情况，以IE解[12]为准，以1 Hz频率(f)为例，对比DQL、QA、QL、LQL近似方法计算的视电阻率、相位以及误差结果。采用直流电法定义的视电阻率和相位：
(27) 
(28) 式中：ρa(ω)为视电阻率；K为装置系数；Ex(iω)为x方向电场；Δd为接收极距；I为电流；φ(ω)为相位。
首先给定异常体电阻率为100 Ω·m，利用QL和LQL近似方法模拟，计算结果及误差曲线见图 3。对于十倍电导率的变化，QL与LQL近似方法都能给出比较准确的结果，整体误差小于5%(图 3)。
|
| f=1 Hz。 图 3 异常电阻率为100 Ω·m时QL、LQL近似方法计算结果及误差曲线 Figure 3 Calculated results and error curves of QL、LQL approximation method when abnormal resistivity was 100 Ω·m |
|
|
减小异常体电阻率，使其为10 Ω·m，计算结果及误差曲线见图 4。当模型的电阻率变化2个数量级时，QL和LQL近似方法的计算误差在异常位置都较大，QL近似方法精度略高于LQL近似方法(图 4)。因此，这2种方法仅适用于计算电阻率变化在2个数量级以内的地电模型。
|
| f=1 Hz。 图 4 异常电阻率为10 Ω·m时QL、LQL近似方法计算结果及误差曲线 Figure 4 Calculated results and error curves of QL、LQL approximation method when abnormal resistivity was 10 Ω·m |
|
|
下面我们讨论DQL近似方法与QA近似方法。同样异常体电阻率为10 Ω·m，计算结果及误差曲线见图 5。DQL和QA近似方法计算精度较高，对所给的地电模型，其计算误差都在6%以内，而DQL近似方法计算最准确，误差在3%以内(图 5)。
|
| f=1 Hz。 图 5 异常电阻率为10 Ω·m时DQL、QA近似方法计算结果及误差曲线 Figure 5 Calculated results and error curves of DQL、QA approximation method when abnormal resistivity was 10 Ω·m |
|
|
降低异常体电阻率，使其为1 Ω·m，DQL与QA近似方法计算结果及误差曲线见图 6。从图 6中可以看出：QA近似方法计算视电阻率的误差在异常部位已经达到近8%；而DQL近似方法计算结果仍然可靠，其误差与异常体电阻率为10 Ω·m时相差不大。
|
| f=1 Hz。 图 6 异常电阻率为1 Ω·m时DQL、QA近似方法计算结果及误差曲线 Figure 6 Calculated results and error curves of DQL、QA approximation method when abnormal resistivity was 1 Ω·m |
|
|
根据上面的计算结果可以看出，QA和DQL近似方法的计算结果精度较高，且DQL近似方法精度最高，当异常电导率是背景电导率的1 000倍时仍能给出精确的计算结果。
以精度最高的DQL近似方法为例，讨论网格剖分疏密程度对拟线性系列方法计算结果的影响。将图 2模型异常体剖分为125个边长为40 m的正方体，利用IE方法进行模拟；然后，将异常体分别剖分为64个边长为50 m的正方体和8个边长为100 m的正方体通过DQL近似方法进行模拟。通过对比这3个模拟结果来分析网格疏密程度对DQL近似方法计算精度的影响，计算结果及误差曲线见图 7。从图 7中可以看出：网格大小对视电阻率的影响不大，计算误差相差不多；而对相位的影响较大，网格越密，相位的计算精度越高。
|
| 异常电阻率为1 Ω·m，f=1 Hz。 图 7 网格剖分不同时DQL近似方法的计算结果及误差曲线 Figure 7 Calculated results and error curves of DQL approximation method by different grid subdivision |
|
|
对于这些方法的计算速度，QA、QL及LQL近似方法都求解标量形式的反射系数，速度不会有太大差别，因此，我们选择QL近似方法与DQL近似方法进行计算时间对比，对比结果见表 1。从表 1中可以看到DQL与QL近似方法计算时间相差不多，所以DQL近似方法更适合用于三维电磁模拟。
| 异常区域剖分单元数 | DQL近似方法正演时间/s | QL近似方法正演时间/s |
| 500 | 4.30 | 4.20 |
| 800 | 7.12 | 6.69 |
| 1 000 | 9.80 | 8.40 |
| 4 000 | 43.46 | 37.80 |
拟线性系列方法的精度高低取决于反射系数计算是否精确。DQL近似方法考虑到了反射系数的方向性，使反射系数的计算更合理；其他3种方法都属于标量方法，我们仍以QL近似方法为例进行讨论。
对图 2模型通过IE法计算异常单元的λx、λy、λz，与QL近似方法求解的λOL进行对比。将异常体剖分为大小相等的8块，选取其中单元a和单元b的反射系数计算结果作为代表，其中心坐标分别为(-50，-50，150) 和(50，50，150)，这2个单元的λx、λy、λz与λOL随
|
| f=1 Hz。 图 8 QL近似方法与IE方法求解的反射系数对比 Figure 8 Comparison of reflection coefficient between DQL approximation method and IE method |
|
|
接下来分析DQL近似方法计算的对角反射系数，与IE法计算的反射系数对比结果如图 9所示。从图 9以看出，DQL近似方法计算的3个方向的反射系数与IE法计算的基本相同。
|
| f=1 Hz。 图 9 DQL近似方法与IE方法求解的反射系数对比 Figure 9 Comparison of reflection coefficient between DQL approximation method and IE method |
|
|
因为不同方向的背景场存在很大差异以及入射角度变化很大，导致各个位置、各个方向的反射强度不同，DQL近似方法考虑了方向的差异性，这是它比QL近似方法精度高的原因。
3 结论本文对拟线性系列方法的精度进行了对比，确定各种方法的应用范围，得出以下几条结论：
1) 对角拟线性近似方法精度最高，适用范围可达数千倍电导率变化的地电模型，因为它考虑了反射系数的方向性。忽略反射系数方向性是拟线性、拟解析、局部拟线性这些标量近似方法较对角拟线性近似方法精度低的原因。
2) 拟解析近似方法在标量拟线性方法中计算精度最高，适用范围为数百倍电导率变化的地电模型。
3) 拟线性与局部拟线性近似方法计算精度较低，其中拟线性近似方法的精度略高于局部拟线性近似方法，它们对电导率变化在2个数量级以内的电性模型能给出可靠的计算结果。
4) 拟线性系列方法的计算速度都比较快，虽然对角拟线性近似方法增加了反射系数的计算量，但整体模拟时间与其他方法相差不多，所以对角拟线性近似方法在拟线性系列方法中适用性最高。
| [1] | Born M. Optic[M]. New York: Springer, 1933. |
| [2] | Berdichevsky M N, Zhdanov M S. Advanced Theory of Deep Geomagnetic Sounding[M]. Amsterdam: Elsevier Science Publishers, 1984. |
| [3] | Torres-Verdin C, Habashy T M. Rapid 2.5-Dimen-sional Forward Modeling and Inversion via a New Nonlinear Scattering Approximation[J]. Radio Science, 1994, 29(4): 1051-1079. DOI:10.1029/94RS00974 |
| [4] | Torres-Verdin C, Habashy T M. A Two-Step Linear Inversion of Two-Dimensional Electrical Conductivity[J]. IEEE Transactions on Antennas & Propagation, 1995, 43(4): 405-415. |
| [5] | Zhdanov M S, Fang S. Quasi-Linear Approximation in 3-D Electromagnetic Modeling[J]. Geophysics, 1996, 61(3): 646-665. DOI:10.1190/1.1443994 |
| [6] | Zhdanov M S, Fang S. Three-Dimensional Quasi-Linear Electromagnetic Inversion[J]. Radio Sci, 1996, 31(4): 741-754. DOI:10.1029/96RS00719 |
| [7] | Zhdanov M S, Hursán G. 3D Electromagnetic Inversion Based on Quasi-Analytical Approximation[J]. Inverse Problems, 2000, 16(5): 1297-1322. DOI:10.1088/0266-5611/16/5/311 |
| [8] | Zhdanov M S, Tartaras E. Three-Dimensional Inversion of Multitransmitter Electromagnetic Data Based on the Localized Quasi-Linear Approximation[J]. Geophys J Int, 2002, 148(3): 506-519. DOI:10.1046/j.1365-246x.2002.01591.x |
| [9] | Ueda T, Zhdanov M S. Fast Numerical Modeling of Multitransmitter Electromagnetic Data Using Multigrid Quasi-Linear Approximation[J]. IEEE Transactions on Geoscience and Remote Sensing, 2006, 44(6): 1428-1434. DOI:10.1109/TGRS.2006.864386 |
| [10] |
孙建国. 高频渐近散射理论及其在地球物理场数值模拟与反演成像中的应用:研究历史与研究现状概述以及若干新进展[J].
吉林大学学报(地球科学版), 2016, 46(4): 1231-1259.
Sun Jianguo. High-Frequency Asymptotic Scattering Thoeries and Their Application in Numerical Modeling and Imaging of Geophysical Fields: An Overview of the Research History and the State-of-the-Art, and Some New Developments[J]. Journal of Jilin University (Earth Science Edition), 2016, 46(4): 1231-1259. |
| [11] |
王伯年, 李过房, 曹伟武. 爱因斯坦求和约定的推广[J].
上海理工大学学报, 2003, 25(1): 5-7.
Wang Bonian, Li Guofang, Cao Weiwu. On the Extension of Einstein's Summation Convention[J]. Journal of University of Shanghai for Science and Technology, 2003, 25(1): 5-7. |
| [12] |
李建平. 带地形的三维复电阻率电磁场正反演研究[D]. 长春: 吉林大学, 2008.
Li Jianping. Research on Electromagnetic Modeling of 3D Complex Resistivity with Topography [D]. Changchun: Jilin University, 2008. |