大地电磁交错采样有限差分二维正反演研究

周武，罗威，蓝星，简兴祥

(1.甘肃省交通规划勘察设计院股份有限公司，甘肃 兰州 730030； 2.成都理工大学 地球物理学院，四川 成都 610000； 3.四川省冶勘设计集团有限公司，四川 成都 610000)

0 引言

大地电磁测深(MT)是探测地球深部电性结构的主要方法，通过在地表同步观测电场、磁场分量，定性或定量地对深部地质构造进行研究，在油气勘探、固体矿产资源勘察、深部地质结构探测、地热和地下水资源调查、地震预测和地质灾害防治等领域应用广泛。在MT二维正演模拟方面，主要计算方法有边界元法、有限单元法和有限差分法等。其中，有限差分法以实现简单、计算量较小和实用性强的优点使用最为广泛[1-3]，但过去的二维有限差分模拟大都将电场和磁场采样点放置于相同位置，而忽略了可能存在的散度问题。交错采样网格的最大特点是能自动保证电磁场分布遵守能量守恒定律[4]，由于三维电磁场散度一般不为零，交错采样网格主要在大地电磁三维正演中被采用[5-6]，而二维情况下的大地电磁交错网格有限差分正演问题研究目前相对较少[7]。

在大地电磁反演问题上，目前基于一阶梯度方向或二阶牛顿方向收敛的各种最优化反演理论已经较为成熟，因此只要拥有良好的正演方法，反演需要考虑的关键问题就是如何计算或处理灵敏度矩阵或海赛矩阵，然后再套用相对成熟的反演方法。目前，大地电磁中主要存在3种主流线性反演方法：OCCAM及其改进型、非线性共轭梯度法(NLCG)和有限内存拟牛顿法(LBFGS)。OCCAM反演[8]及其改进型是基于梯度收敛的最优化法，与最小二乘不同，OCCAM在每次迭代中都要求取最佳的目标函数正则化因子，因此其反演比较稳定，但计算速度很慢，同时由于在反演过程中必须存储全部灵敏度矩阵，内存需求较大。NLCG[9]是介于最速下降法和牛顿法之间的方法，利用梯度的共轭方向提升收敛效率，避免求取二阶海赛矩阵；为避免存储灵敏度矩阵，NLCG利用正演的稀疏矩阵求取梯度向量，因此对内存需求很小，但NLCG在每次迭代中都需进行线搜索求取最佳步长，因此一次反演迭代中有可能进行多次正演。LBFGS[10]是由拟牛顿法发展起来的一种反演方法，拟牛顿法根据更新海赛矩阵的秩不同有非常多的变种，主要有DFP法和BFGS法。LBFGS又是BFGS法的变种，其特点是使用有限次迭代的梯度和修正量信息来更新海赛矩阵，对内存需求小，同时基本和BFGS迭代一样稳定。LBFGS和NLCG相比，二者内存需求差不多，求取梯度的过程一样，唯一的不同是LBFGS在迭代中后期基本上每次迭代的修正方向都为牛顿方向，其修正步长基本为1，这样就可一定程度上节省线搜索时间，其用于线搜索的正演次数约为NLCG的一半。近年来，LBFGS法在大地电磁二三维反演中开展了诸多研究[7,11-13]，已成为目前大地电磁二三维反演的主要方法之一。

本文首先利用交错采样网格实现大地电磁二维正演，然后基于LBFGS法，实现大地电磁二维交错采样网格有限差分LBFGS反演，目标是实现一种精度高、效率高的大地电磁二维正反演算法，并能有效应用于实测大地电磁资料的反演。

1 大地电磁交错网格有限差分二维正演

在大地电磁研究的频率范围内(约105～10-4Hz)，由于频率相对较低，满足似稳电磁场，忽略位移电流后的Maxwell积分形式方程为：

∮H·dl=∬σE·dS，

(1)

∮E·dl=∬iωμ0H·dS。

(2)

式中：E和H分别为电场强度和磁场强度矢量，真空磁导率μ0=4π×10-7H/m，σ为介质电导率，ω为圆频率。

对于二维大地大地电磁，需要垂直测线方向无限延伸。对于第(i,k)个单元的电场分量，根据式(1)可由周边的4个磁场求得，这4个磁场分量，根据式(2)又可由各自周边的2个电场求得。综合起来，即这单元的电场可由其周边的4个电场分量求得，然后按这种关系以此类推，将所有单元的电场与各自周边4个电场的线性方程组合起来，就构建了刚度矩阵。图1为直角坐标系中的第(i,k)个单元差分示意。

1.1 TE模式正演

对于TE模式，以Ey(i,k)为中心，式(1)Maxwell第1个方程展开为：

Hz(i,k)Δz(k)+Hx(i,k+1)Δx(i)-

Hz(i+1,k)Δz(k)-Hx(i,k)Δx(i)=

σ(i,k)Ey(i,k)Δz(k)Δx(i) 。

(3)

对式(3)中的4个磁场，由式(2)又可分别用周围的电场求得：

Ey(i,k)-Ey(i-1,k)=

(4)

Ey(i,k-1)-Ey(i,k)=

(5)

Ey(i+1,k)-Ey(i,k)=

(6)

图1 MT二维交错网格有限差分离散示意Fig.1 MT 2-D staggered grid finite difference separation diagram

Ey(i,k)-Ey(i,k+1)=

。(7)

将这4个关系式代入式(3)中即可消去磁场分量，最终得到Ey(i,k)与其周边4个电场的关系式。用C1～C5的系数表示为

C1Ey(i-1,k)+C2Ey(i,k)+C3Ey(i+1,k)+

C4Ey(i,k-1)+C5Ey(i,k+1)=0 ，

(8)

其中：

C2=-C1-C5-C4-C3-σ(i,k)iωμ0Δz(k)Δx(i) ，

1.2 TM模式正演

对于TM模式，以Hy(i,k)为中心，式(2)Maxwell第2个方程展开为

Ez(i,k)Δz(k)+Ex(i,k+1)Δx(i)-

Ez(i+1,k)Δz(k)-Ex(i,k)Δx(i)=

iωμ0Hy(i,k)Δz(k)Δx(i) 。

(9)

对式(9)中的4个电场，由式(1)又可分别用周围的磁场求得：

Hy(i,k)-Hy(i-1,k)=

(10)

Hy(i,k-1)-Hy(i,k)=

(11)

Hy(i+1,k)-Hy(i,k)=

(12)

Hy(i,k)-Hy(i,k+1)=

(13)

将这4个关系式代入式(9)中即可消去电场分量，最终得到Hy(i,k)与其周边4个磁场的关系式。用C1～C5的系数表示为：

C1Hy(i-1,k)+C2Hy(i,k)+C3Hy(i+1,k)+

C4Hy(i,k-1)+C5Hy(i,k+1)=0 ，

(14)

其中：

C2=-C1-C5-C3-C4-iωμ0Δz(k)Δx(i) ，

1.3 方程组

根据前述TE和TM离散方案，最终可以获得关于Ey和Hy的离散方程：

(15)

式中：Ey和Hy为求解向量，be、bh对应的右端项向量，Ae、Ah分别为式(8)、式(14)中不同单元相关的C系数组成。

要想求解研究区域内所有采样点上的场值，还需给出研究区边界上场值，交错采样网格二维边界条件可简要概括为表1所示形式[2]，其中顶边界为第一类边界条件(其中TE模式顶面位于空气顶层，TM模式顶面位于地表)，侧边界为第二类边界条件，底边界对应第三类边界条件(其中σ为底层介质电导率，Zbottom为底层顶面阻抗)。对构建的方程组，由于直接解法可获得精确解，因此这里采用基于LU分解的直接法求解。

表1 边界条件形式

2 LBFGS反演

根据Tikhonov的正则化反演理论[14]，反演问题可表述为

ψ(m)=ψd(m)+λψm(m)=

(WdΔd)T(WdΔd)+λ(Wmm)T(Wmm) 。

(16)

式中：m为模型;ψd(m)为数据目标函数;ψm(m)为模型约束目标函数;λ为正则化因子;Δd为观测数据与模型响应的差向量;Wd为数据加权矩阵;Wm为模型约束矩阵，目的是使反演得到的模型具备所设定的相应结构特征，一般采用二阶拉普拉斯算子。

在确立反演目标函数后，接下来的工作就是选择合适的方法进行最优化反演。大地电磁反演方法较多，维性、反演效率以及对计算资源需求等因素影响着反演方法的选择。经典的最优化方法基本都是基于梯度下降法和牛顿法的改进或重构，在大地电磁中几乎所有的最优化反演方法都是这两种方法的实用化或变种，其中有限内存拟牛顿法(LBFGS)因其特殊的优势，近年来逐渐应用广泛并受到认可[12-15]。LBFGS算法流程如下:

1)设定k=0，初始模型m0和精度阈值ε；

2)设置初始的正定矩阵H0=I；

5)若‖gk+1‖<ε，则算法结束；

6)利用m个已保存的s、y向量对计算Hk+1：

7)迭代次数k=k+1，跳转到步骤3继续。

3 算例分析

3.1 层状模型正演

由于层状介质模型可通过一维正演求得其解析解，因此二维正演算法的数值解可采用一维解析解进行对比验证。层状模型第一层电阻率100 Ω·m，厚1 km，第二层电阻率1 000 Ω·m，厚9 km，第三层电阻率10 Ω·m。

计算过程中非扩展区域网格数为横向50列，每列宽度10 km，纵向80层，层深度对数递增，网格剖分单元数4 000个，测点均匀分布在地表每列单元中心共100个，计算了104～10-4Hz共61个频率的视电阻率和阻抗相位响应。图2给出了解析解与剖面中心位置测点的视电阻率和阻抗相位对比曲线，可以看出数值结果和解析解曲线基本重合，仅在约1 000 Hz以上高频部分存在细小差异，但这是有限差分数值计算普遍均存在的现象。二维数值解与解析解视电阻率平均相对误差为1.19%，相位平均相对误差为0.72%，表明相对误差较小，二维正演结果精度较高。

图2 层状模型二维正演中心测点视电阻率和阻抗相位Fig.2 Apparent resistivity and impedance phase of 2-D forward modeling center of layered model

3.2 异常体模型正反演

为进一步测试正演算法的正确性以及反演算法的稳定性和效率，设计了较为复杂的异常体模型。模型的设计参考借鉴Egbert对Modem算法的测试模型[12]，设计由浅至深为2层，横向上3列，共6个交错分布的高低阻异常体，模型共80层×50列=4 000个单元，50个地表测点按100 等距分布，图3为模型的具体形态。模型背景电阻率100 Ω·m，纵向上在150～300 m、800～1 500 m深度存在异常体，位置交错分布，电阻率值分别为1 000 Ω·m、10 Ω·m。计算了104～1 Hz共30个频点的正演响应，图4为正演响应拟断面，可以看到，无论视电阻率还是阻抗相位，都对高阻和低阻异常有较明显的反映。

图3 异常体模型Fig.3 Anomalous body model

图4 异常体模型正演拟断面Fig.4 Forward section of abnormal body model

对异常体模型进行反演计算。反演前在正演数据中加入3%的高斯随机噪音，并将其作为数据协方差，反演正则化因子设为 0.01，初始模型为100 Ω·m的均匀半空间。反演终止条件为拟合差RMS不大于1，这里的拟合差计算公式为：

(17)

分别对合成数据开展TE、TM和TE+TM反演，最终反演的步长曲线和拟合差曲线如图5所示。其中，TE模式反演步长经过8次后稳定为1，TM模式9次后稳定，TE+TM模式29次迭代后稳定为1。3种模式反演的拟合差RMS从初始的约11开始，经过几次迭代后快速衰减到3左右，随后缓慢降低，其中TE、TM和TE+TM分别在第17、23、40次迭代后满足收敛条件。从反演步长曲线和拟合差曲线看出反演较为稳定，符合LBFGS收敛特征。

TE、TM及TE+TM3种模式LBFGS最终反演结果如图6所示。与真实模型对比看出，3种模式的反演结果总体上与真实模型都较为吻合：TE模式对异常体厚度反映较好，但横向边界相对一般；TM模式对异常横向范围反映较好，但纵向分辨相对一般；TE+TM模式综合了TE和TM的优点，对异常体反映较好。上述特征符合MT二维反演规律，表明算法可靠。

另外，为分析LBFGS反演的效率，将该模型同时开展了NLCG二维反演，表2给出了2种算法的时间对比。对于NLCG反演，TE、TM、TE+TM的10次迭代耗时分别为75、57、133 s，而LBFGS的10次迭代耗时分别为59、45、88 s，相对有约27%的效率提升。NLCG拟合差达到1的计算时间分别为175、92、506 s，而LBFGS只用了125、73、371 s，相对有约25%的效率提升。这表明LBFGS反演效率总体上比NLCG高。

图5 异常体模型二维反演步长和拟合差曲线Fig.5 2-D inversion step length and fitting difference curve of abnormal body model

图6 异常体模型二维反演Fig.6 2-D inversion of anomalous body model

表2 LBFGS和NLCG反演效率对比

3.3 实测资料反演

图8为采用本文算法对实测MT资料的TM模式反演结果，分析看出，剖面范围内发育多条断层，断层位置和形态与地表地质调查以及氡气测量结果较为一致。此外，地表的诸多断裂在约1 000 m深度汇集，形成F1和F2两个主要断层，发育深度均大于2 000 m，其中F1倾向南西、 F2倾向北东，2个断裂在深部交汇，断层深部影响带范围较大。F1和F2断层带中间低阻现象较为明显，表明断层带上盘局部岩体破碎带较发育，为地热流体的赋存提供了一定的储集空间，具有一定的导水、导热作用。地热流体沿断层带上升后储存在热储含水层内，在上覆隔水保温盖层的作用下，形成地热资源。通过实测MT资料的反演和解释，可以看出本文算法的反演结果揭露的特征，与地质和相关资料较为一致。

图7 MT测线分布Fig.7 MT line distribution

图8 MT反演结果Fig.8 MT inversion results

4 结论及讨论

本文借鉴三维大地电磁正演模拟中使用的交错采样网格方法，采用交错采样网格开展了大地电磁二维有限差分正演研究，模型计算表明该正演算法具有较高的计算精度。基于交错采样网格大地电磁二维有限差分正演，利用有限内存拟牛顿算法，实现了大地电磁二维快速反演。理论模型计算表明有限内存拟牛顿反演算法具有较好的稳定性，与非线性共轭梯度算法对比发现，LBFGS计算效率至少大于NLCG效率的1.2倍。将反演算法在宕昌县官鹅沟地热勘查项目中开展了大地电磁资料二维反演解释，结果表明本文算法与现有资料和认识较为一致，具有一定的实用价值。

由于算法暂未考虑地形以及各向异性介质的影响，建议下一步开展基于交错采样网格的起伏地形、各向异性正反演研究。此外，尽管LBFGS具有较高的计算效率，但随着数据量和模型单元的增大，并行计算研究仍然值得开展。