1、第 37 卷第 1 期2023 年 2 月现代地质GEOSCIENCEVol.37No.1Feb.，2023DOI:10.19657/j.geoscience.1000 8527.2021.156时移可控源音频大地电磁法三维反演研究胡琪璇，谭捍东，于翠(中国地质大学(北京)地球物理与信息技术学院，北京100083)收稿日期:2021-08-27;改回日期:2022-09-30。基金项目:国家自然科学基金项目(41830429)。作者简介:胡琪璇，女，博士研究生，1997 年出生，地球探测与信息技术专业，主要从事地球物理算法研究。Email:huqixuan email cugb edu cn。
2、通信作者:谭捍东，男，教授，博导，1966 年出生，地球探测与信息技术专业，主要从事电法勘探理论及应用研究。Email:thd cugb edu cn。摘要:时移地球物理监测是观测地下物性参数随时间动态变化的有效方法。地球物理反演结果会受到测量误差、测量环境的噪声污染以及地球物理反演多解性等系统性因素的影响。如果对时移地球物理数据进行不同时刻数据的单独反演，可能会导致不同时刻的反演结果可对比性差，从而影响时移地球物理的监测效果。针对可控源音频大地电磁法在监测领域的实际需求，将不同时刻的观测数据放在一起反演，不同时刻的电阻率模型相互约束，实现了时移可控源音频大地电磁法三维反演。设计三种类型模型进
3、行合成数据的单独反演和时移反演对比试算，验证了时移反演算法的有效性。反演结果表明:通过加入相邻时刻模型约束，时移反演可以更好地聚焦异常体的位置，降低噪声以及测量环境变化对反演结果的影响，防止反演伪影掩盖真实的地下变化。研究成果为时移可控源音频大地电磁法的实际应用奠定了良好基础。关键词:时移可控源音频大地电磁法;有限内存拟牛顿反演;时移约束中图分类号:P631.3文献标志码:A文章编号:1000 8527(2023)01 0090 093D Inversion of Time-lapse Controlled SourceAudio-frequency MagnetotelluricsHU Qi
4、xuan，TAN Handong，YU Cui(School of Geophysics and Information Technology，China University of Geosciences，Beijing100083，China)Abstract:Time-lapse monitoring is an effective way to observe underground dynamic changes Geophysical in-version results could be affected by systemic factors，including measure
5、ment errors，noise pollution in the stud-ied environment，and the inversion ambiguity If time-lapse geophysical data is inverted separately at differenttimes，the results may have poor comparability，which would affect the monitoring of time-lapse geophysics Forthe actual needs of the controlled-source
6、audio frequency magnetotelluric(CSAMT)in the monitoring field，weintegrated the observation data at different moments for inversion，and constructed resistivity models at differentmoments(mutually constrained)for three-dimensional inversion of the time-lapse CSAMT Three types of mod-els were designed
7、to perform without time-lapse constraint inversion and time-lapse inversion comparison trialcalculations of synthetic data，which verified the inversion algorithm and its effectiveness The inversion resultsshow that by adding adjacent time model constraints，the time-lapse controllable source inversio
8、n can better fo-cus on the anomaly locations，reduce the noise impact and local environment changes，and minimize the impactof inversion artifacts Our findings have laid a good foundation for the application of time-lapse CSAMTKey words:controlled-source audio-frequency magnetotellurics;limited-memory
9、 quasi-Newton method;time-lapse constraint0引言地球物理方法能够反映地下构造物性参数特点，多用于监测地面及地下条件的变化。在同一位置的不同时间收集相关的数据，通过识别和计算延时图像之间的差异，可以直观评估地表以及地下条件的变化。然而，由于测量误差、测量环境的噪声污染以及地球物理反演的非唯一性的影响，这些延时图像中会出现假异常。此外，这些假异常可能在图像中被显著放大，从而掩盖地质条件的真实变化。为解决这一问题，国内外学者对时移反演进行了深入的研究1 8。学者 Labrecque 和 Yang9 修改了电阻率成像的三维 Occam 反演算法，对背景和后续数
10、据集之间的差异进行反演。Oldenborger等10 11 提出使用初始数据集的反演模型作为后续时间的先验信息。Kim 等12 13 将地下结构和整个监测数据定义在一个时域中，获得一个四维时空模型。Karaoulis 等14 15 提出四维自适应时移电阻率法反演。国内学者马欢等16 利用归一化方法处理时移电阻率法数据，有效分辨局部电阻率结构随时间的微小变化。此前时移反演的研究主要集中在电阻率法，已发表的研究很少有将电磁感应技术考虑在内17 18，基于可控源音频大地电磁法(CSAMT)具有高信噪比、勘探范围广、三维正反演算法发展成熟等特点19 20，本文将时移技术引入可控源音频大地电磁法领域，实
11、现其勘探性能的提高。为体现参数随时间的变化，本文建立多组动态变化的地电模型，分别对每个时刻模型进行常规电磁法正演数值模拟，汇总不同时刻的响应得到时移反演初始数据。在反演目标函数中加入时移项，利用修改的 CSAMT 时移三维有限内存拟牛顿反演程序对合成数据进行反演，通过分析得到的三维反演结果，来讨论时移方法的有效性。1时移可控源音频大地电磁法三维反演理论1.1可控源音频大地电磁法三维正演时移电磁法的正演基本方程同样满足麦克斯韦方程，但是方程中的参数都随时间变化，实现时移电磁法正演算法只需建立多组表征电阻率动态变化的地电模型，然后进行多次常规电磁法正演模拟。正演采用将总场分为一次场与二次场叠加的策
12、略。取谐变因子为 eit时，二次场满足的麦克斯韦方程如下21 22:(Eb)i 0Eb=i 0(0)Ea(1)式中:0为背景电导率，0为背景磁导率，为圆频率，Ea为一次电场，Eb为二次电场。图 1单网格剖分示意图Fig.1Schematic diagram for the single grid division二次电场采用交错网格有限差分法23 求解。网格离散策略采用如图 1 所示的电场单元模式，单元体的棱边中点为电场分量，面中心为磁场分量。将微分方程(1)分别在 x、y、z 方向进行离散，单元网格每个棱边上的电场可由其周围 12 个电场分量表示，并且依据网格剖分数形成 Nx(Ny1)(Nz
13、 1)+Ny(Nx 1)(Nz 1)+Nz(Nx1)(Ny1)个类似的方程，汇总整理得到如下线性代数方程组:KEb=S(2)其中:K 为大型对称稀疏系数矩阵，S 为与边界条件及一次场有关的向量，Eb为待求的二次电场值。通过采用最小残差法(QM)对方程组(2)进行求解，得到每个节点的磁场和电场分量，利用体积加权插值得到地面观测点阻抗(3)，根据计算公式(4)可求得地面观测点的视电阻率值和相位值。Zxy=Ex2Hx1 Ex1Hx2Hx1Hy1 Hx2Hy1(3)xy=1 0Zxy2xy=tan1Im(Zxy)e(Zxy)(4)其中:xy为视电阻率值，xy为相位值，Zxy为张量阻抗。Ex1、Hx1、
14、Hy1为源场1 的作用下产生的电磁场总场的分量值。Ex2、Hx2为源场 2 的作用下产生的电磁场总场的分量值。19第 1 期胡琪璇等:时移可控源音频大地电磁法三维反演研究1.2时移可控源音频大地电磁法三维反演1.2.1时移反演目标函数时移反演的目标函数定义为:(M)=d(M)+m(M)+t(M)(5)其中:d(M)为数据拟合差项，m(M)为模型圆滑约束项，为模型圆滑约束项的正则化因子，t(M)为模型时间约束项(简称时移项)，为时移项对应的正则化因子。M=m1，m2，mntT是时移反演所有时刻的模型参数，nt表示监测时间的个数，mnt对应 nt时刻的模型参数。d(M)的表达式为:d(M)=Wd
15、D F(M)T Wd D F(M)(6)其中，D=d1，d2，dntT是用于反演的数据集，维数为(nd*nt)1，dnt对应 nt时刻的数据，nd表示单一时刻的数据个数;F(M)为正演响应;Wd为数 据 加 权 矩 阵:Wd=diag 1/1，1/2，1/(nt*nd)m(M)的表达式为:m(M)=(M Mref)TC1m(M Mref)(7)其中，Mref为与 M 对应的时移反演所有时刻的初始模型参数向量，两者维数都为(nm*nt)1，nt表示监测时间的个数，nm表示每个模型网格剖分的个数;C1m为模型线性光滑矩阵。t(M)的表达式为:t(M)=ntt=1 mt mt12=MTLTLM(8)
16、其中，mt对应时刻 t 的模型参数，L 为系数矩阵:L=I1 I20000I2 I30000I3 I40000000Int1 Int(9)时移反演通过在目标函数中加入时移项来约束相邻时刻模型的变化。每一次迭代都同时进行若干个时刻剖面数据的反演，并同时获得多个时刻的反演结果，通过相邻时刻数据之间的相互约束模拟地下真实的变化情况。1.2.2时移反演流程反演采用有限内存拟牛顿反演法(LBFGS)，对目标函数(M)求导，在第 k 次迭代处进行泰勒级数展开，并且忽略二阶以上导数，得到:(Mk+1)=(Mk)+(Mk)(Mk+1 Mk)=0(10)通过变换，可以得到模型迭代公式:Mk+1=Mk1(Mk)(Mk)=Mk k Hk1(Mk)(11)由此得到模型每一次反演迭代需要更新的三个参数:目标函数梯度(Mk)(表示为 gk);海森矩阵的逆 Hk1;迭代步长 k。目标函数梯度可表示为:(M)=2(WdJ)T Wd D F(M)+2 C1m(M Mref)+2 LTLM(12)其中，和 为正则化因子。LBFGS 不直接求解雅可比矩阵 J，而是求解雅可比矩阵与数据拟合差的乘积:(WdJ)T Wd D F
