基于边界约束有限内存的拟牛顿CSAMT一维反演及应用_唐传章.pdf

第 47 卷 第 5 期 煤田地质与勘探 Vol. 47 No.5 2019 年 10 月 COAL GEOLOGY 2. Key Laboratory of Exploration Technologies for Oil and Gas Resources, Yangtze University, Wuhan 430100, China Abstract CSAMT is characterized by time-consuming calculation and low accuracy. The conventional inversion has simple boundary processing, needs big memory, it is difficult to have accurate search direction, influ- encing severely the effect of inversion. In order to achieve rapid and high-precision forward modeling of CSAMT, to improve the boundary processing , to reduce memory usage, to optimize the searching direction to reduce the times of iterative inversion, LBFGS CSAMT 1D inversion based on boundary constraint and limited memory was proposed. Firstly in CSAMT forward modeling, double-precision and parallel algorithm was adopted, realizing rapid and high-precision calculation of the horizontal electric field excited by a wire source with limited length under conditions of 1D medium model. Secondly in inversion two objective functions based on the relative error and the absolute error were constructed respectively. The constraints of smooth model were intro- duced. Self-adaptive regularization strategy was adopted to renewal the regularization factors. LBFGS-B algorithm ChaoXing 194 煤田地质与勘探 第 47 卷 of limited memory based on boundary constraint was adopted. 1D precise and rapid inversion of CSAMT was re- alized. The amplitude of the horizontal electric field Ex of the geoelectric model of multiple thin layers with low and high resistivity was used as inversion data to conduct direct inversion and model test. The results illustrate that this can reflect the trend of the variation of ation resistivity with depth, and has higher resolution for thin layers with high and low resistivity, the resolving effect is better for a thin layer with low resistivity. The resolution of the inversion profile of the actual CSAMT data is superior to that of the conventional inversion for continuous medium. The geological effect is evident, indicating that this has good application prospect. Keywords CSAMT forward modeling; parallel algorithm; objective function; LBFGS-B inversion 可控源音频大地电磁测深法Controlled Source Audio Magnetotelluric ， 简称 CSAMT是电磁 勘探中一种快速有效的可控源电磁勘探方法，已广 泛用于油气及矿产资源勘查等领域。通常，CSAMT 资料处理采取了远区处理方法，即不考虑源的影响而 计算视电阻率，然后套用 MT 反演与解释方法[1-2]。由 于可控源的三维复杂性，近场影响是客观存在且 不可忽视的。尽管采取了多种校正方法[3]，但效 果有限。最有效的方法是直接利用全区场分量进 行反演[4-7]。尚通晓等[8]将阻尼最小二乘法应用到可 控源频率电磁法数据反演中，克服了近区和中间区 卡尼亚视电阻率畸变的影响；朱威等[9]在此基础上 进行模型约束反演，分析研究了拉格朗日算子对反 演结果的影响；而在反演效率以及灵活性等方面， 周军等[10]，汤井田等[11]也做了相关研究。刘颖等[12] 基于高斯–牛顿算法实现了海洋可控源电磁场一维 反演；堃王鹏等[13]则使用自适应正则化方法反演 CSAMT 全区资料。由于线性逼近算法忽略了目标 函数二阶导数信息，且比较依赖模型初始值，极容 易陷入局部最优，而以遗传算法、模拟退火算法及 拟牛顿算法等为代表的非线性算法具有全局收敛 性，克服了线性迭代算法的缺点，是反演研究的重 点发展方向。李帝铨[14]利用遗传算法反演 CSAMT 全场资料，并在反演时引入最小构造模型作为反演 目标函数的罚函数，提高了反演模型稳定性和收敛 性。拟牛顿法BFGS通过近似 Hessain 矩阵达到减 少计算量的目的，收敛效率更高。当处理大规模数 据时，拟牛顿法需要消耗大量的内存，为解决此问 题，D. C. Liu 等[15]提出了有限内存 BFGSLBFGS 算法，此算法只需储存有限迭代的最新数据便可获 得近似 Hessain 矩阵，极大节省了内存，同时反演 稳定，已成功应用于一维可控源电磁法反演问题中 [16-17]。 基于边界约束有限内存的 LBFGS_B 算法[18-21] 继承了 LBFGS 的优良特性， 可自定义模型参数上下 边界，从而对模型参数进行简单的边界约束处理， 进而减少反演中的非唯一性。 由于 CSAMT 正演的高精度计算十分耗时，本 文对正演中基于高斯勒让德积分的有限长导线源一 维水平电场采取并行计算，然后尝试将 LBFGS_B 算法应用于有限长导线源一维 CSAMT 反演中，并 引入最小构造模型约束，同时采用自适应正则化因 子策略，直接对电场分量进行反演，研究不同地电 模型条件下反演方法的稳定性与可行性。此外，也 比较了基于相对误差和绝对误差的两种不同目标函 数的反演效果。 1 一维正演理论 长度 2L 的发射导线，铺设于地表。若取发射导 线的中心点为坐标原点，有限长导线源在地面上的 电场 x 分量 Ex表达式如下  201 010 2* 11 1 i1 coscos2 /dd d 2π E x Pnm mJmrJmrrmJmrml mnRmnRk R               E 1 式中 PE Edl，PE为偶极矩；r 为收发距；θ 为源方 向与接收点到发射源中点方向的夹角；ω 为圆频 率；μ0 为真空中的磁导率；J0和 J1分别为 0 阶和 1 阶贝塞尔函数；m 为积分变量；k1为波数； 11 1 11 2 2 coth[artanhcothartanh] N N nn Rn hn h nn  ； *111 1 1 11 2 2 2 coth[actanhcothartanh] NN NN nn Rnhnh nn    ； 22 ii nmk， 2 00 i/ ii k 。 式中 hi为第 i 层地层厚度，m；ρi为第 i 层地层电阻 率，Ωm；ki为第 i 层的波数。 有限长接地长导线的电磁场求解可以采用高斯 积分方法计算， 即首先求取长度为 2L 的有限接地长 导线的高斯积分节点，其次以高斯积分对每一个节 点可看作电偶极子进行积分，从而得到对应区间 的电磁场响应，最后将所有区间的电磁场与高斯积 分节点相对应的权重相乘并矢量叠加，得到有限长 接地长导线的电磁场响应。由于有限长导线源一维 正演需要大量的迭代运算，运算时间长，因此为减 ChaoXing 第 5 期 唐传章等 基于边界约束有限内存的拟牛顿 CSAMT 一维反演及应用 195 少反演时间，本文采取了简单的并行策略，利用 OpenMPI 进行并行编程，实现每个频率电磁场响应 的并行计算，其流程如图 1 所示。 图 1 一维有限长导线源 CSAMT 正演流程图 Fig.1 Flow chart of CSAMT forward modeling of one-dimensional limited-length wire source 2 基于 LBFGS_B 一维反演 2.1 目标函数构造 本次反演将直接反演有限长导线源所产生的水 平电场信号。为引入平滑约束条件，其最小构造目 标函数如式2所示  d ,1,, iii lmu iN ≤≤mmm 2 式中 φm为反演目标函数；φdm为数据误差目标 函数；Θm为模型约束函数，α 为正则化因子，mi 为模型向量，N 为模型参数个数，li、ui分别为模型 参数上下边界。 为了研究不同方差函数 φdm构造目标函数的 反演效果，本文分别采取了绝对误差和相对误差目 标函数构建，即  2 d 1 n xixmi i EE    m 3  2 d 1 n xixmi i xi EE E        m 4 式中 n 为频点个数，Exi为归一化的观测数据，Exmi 为模型正演电场水平分量。 关于模型约束函数 Θm 的处理，采取对模型 平滑约束方法[4]，即 22 00 d d dln d dln d m zm z z zz zz    m 5 对式5进行离散有 1 1 TT2 1 11 [ ] N jj jj jjj zz mm hh         ＝mm L Lm 6 式中 L 为N–1N 加权矩阵，其具体表达式为 11 11 0 00 NN          L   其中， 11 / jjjjj zzhh  。反演时按对 数等间隔将地电模型化分N层即与频点个数相等， 这样厚度参数不参与反演计算。 2.2 LBFGS_B 算法原理 假设目标函数二次可微，对已知模型参数进行 二阶泰勒展开，有  T T1 2 kkkkkkk f gmmmmmmBmm 7 式中 gk为目标函数的一阶导数，Bk为二阶导数， k为当前迭代次数。每次迭代开始时，当前迭代 mk，目标函数值φk，梯度值gk，近似矩阵Bk都将 给出。 对约束最小化问题，利用梯度投影方法寻找一 组可行域，随后使目标函数在可行域内参数最小。 为此，先考虑一组分段线性函数，然后获取最速下 降方向在可行域上的投影。  , , kk m tP mtgl u 8 其中， , , ,, , iii iiii iii l ml P m l um lmu u mu       ≤≤，t为最速下 降步长。 计算前迭代一元分段函数 k m t的局部极小 点，可称为Cauchy点 c m。将超过可行域边界的模 型参数定义为集合 c A m，定义非 c A m集合为自 由变量子空间。然后求取集合 c A m的解或者近似 解，利用直接或者迭代算法寻找目标函数在自由变 量子空间中的解。当得到两集合的近似解 1km时， 进行第1k 迭代，沿着 1kk k dmm 方向进行线性 搜索， 并得到满足一定条件下的步长以及新解 1k m  。 在LBFGS_B算法中，近似矩阵Bk计算公式如下式 所示 T kkkk BIW M W 9 其中，θ为当前迭代的缩放因子；I 为单位矩 阵； ][ kkk Wsy； 1 T T kk k k kk          ds M l s s ；[, kk n Yy ChaoXing 196 煤田地质与勘探 第47卷 1 ,] k y； 1 [,,] kk nk Sss，其中n 为LBFGS_B 算法中最大保留向量数； 1kkk smm  ；ykgk1–gk； TT 11 diag,, kk nk nkk sysy     D；   T 11 , if 0else k nik ni k i j syij L          。 此算法已由 Zhu C.等[22]，J. L. Morales等[23]编 程实现，并开源发布。其算法基本流程如图2所示。 图 2 LBFGS_B 反演算法流程图 Fig.2 The inversion flow chart of LBFGS_B algorithm 3 模型检验 采用三层、五层模型进行反演效果检验。模型 正演时，发射源长度为2 km，电流大小为1 A，发 射源中心点坐标为0,0，接收点坐标为100,6 000。 频点41个，频率范围为1104 Hz，呈对数等间隔。 反演过程中，正则化因子α采取陈小斌等[23]提出的 自适应策略，根据第一个高频频点对应的视电阻率 ρnf，利用 nf 1 356H f  计算反演总深度，按照频 点个数对地层进行对数等间隔划分。反演均以均匀 半空间为初始模型，初始模型电阻率皆为最后一个 频点对应的视电阻率，反演时每层电阻率的搜索范 围为1104 Ωm。迭代停止条件如下 以相对误差构造的目标函数，终止条件为 2 1 rms0.001 n xixmi i xi E E n E       10 以绝对误差构造的目标函数，终止条件为  2 1 rms0.001 n xixmi i E n E     11 3.1 三层模型反演 为考察此算法对低阻薄层的分辨能力，三层H 型地电模型中间层设置为低阻薄层，参数如表1所 示，图3为反演结果。根据图3a所示，两种目标函 数的拟合误差皆小于1， 反演的水平电场数据与理 论模型正演场值基本重合。 根据图3b反演结果可知， 能够准确地反映电阻率随深度变化， 说明此方法对低 阻薄层具有较高的分辨能力。图3c表示两种目标函 数收敛情况， 由图可知两种目标函数在反演前期收敛 速度快， 但在后期速度慢。 相对误差构造的目标函数 收敛效果低于绝对误差构造的目标函数。 为考察此算法对高阻薄层的分辨能力，设置三 层K型地电模型，中间层为高阻薄层，参数如表1 所示，反演效果如图4所示。根据图4a所示，两种 目标函数的拟合误差皆小于1， 正反演数据拟合较 好。 根据图4b反演结果看出反演模型能够反映出理 论模型电阻率随深度变化的趋势，对比发现基于式 3的反演效果更好，但两者都不能分辨出高阻薄 层。根据图4c可知，绝对误差构造的目标函数仅需 迭代28次便可退出反演。 两种目标函数调用正演次 数稳定12次，也说明算法的稳定性。 表 1 H 和 K 型层状介质模型参数 Table 1 The layered medium model of H and K type 模型类型 电阻率 ρ/Ωm 厚度 h/m 100 1 000 10 100 H 200 100 1 000 1 000 100 K 500 3.2 复杂层状模型反演 为检验LBFGS_B算法在复杂地层模型的适应 能力，本文设计更为复杂的五层模型，其模型参数 如表2。 ChaoXing 第5期 唐传章等 基于边界约束有限内存的拟牛顿CSAMT一维反演及应用 197 图 3 H 型模型反演结果 Fig.3 The inversion results of H-model 图 4 K 型模型反演结果 Fig.4 The inversion results of K-model 表 2 五五层 KHA 型层状介质模型参数 Table 2 The layered medium model of 5-layers of KHA types 电阻率 ρ/Ωm 厚度 h/m 300 100 900 250 50 350 600 500 1 000 图5所示是KHA层状介质模型反演结果。根 据图5a所示， 两种目标函数都能拟合理论模型观测 数据，相对误差构造的目标函数在高频和低频处拟 合效果相对较差。 根据图5b反演结果看出反演模型 能够反映出理论模型电阻率随深度变化的趋势，分 层能力较好，皆能反演出高阻层间的低阻夹层。根 据图5c所示，算法调用正演次数并不稳定，但两种 目标函数的函数值依然在减小，说明LBFGS_B算 法依然有寻找最优解的能力，而且此算法具有较高 的地层分辨能力。 4 实际资料处理结果 试验测线位于新疆和静县古沦沟。该区域位于 哈萨克斯坦–准噶尔板块与塔里木板块的结合部位， 古代壳幔作用强烈，物质交换频繁，有较好的成矿 条件，是我国重要的成矿带之一。据初步掌握的地 质矿产调查资料，区内已知的金、铁、铜、煤等各 种矿点达20多处，显示该区矿产资源开发潜力巨 大。为全面、系统地掌握古伦沟区域地质矿产资料， 初步查明了区域矿产资源潜力，实施矿产资源综合 开发计划，图6所示为测线位置，在该区域共布设 了5条CSAMT测线。 ChaoXing 198 煤田地质与勘探 第47卷 图 5 KHA 型模型反演结果 Fig.5 The inversion results of KHA-model 图 6 测线位置示意图C1-C5 为 CSAMT 测线 Fig.6 Location of the survey lines 应用LBFGS_B一维可控源电磁反演方法，对 研究区C5测线CSAMT实测资料进行了反演。图7 为C5测线1号测点的LBFGS_B反演模型的正演水 平电场Ex与实测水平电场资料的拟合效果图。从 图7看出，通过对电阻率模型参数进行边界约束， 高频段拟合度较高，总体拟合度优良。为了对比反 演 效 果 ， 同 时 进 行 了 常 规 一 维 快 速 反 演 [24]和 LBFGS_B反演，反演电阻率深度剖面如图8所示。 从图8可以看出，LBFGS_B反演剖面图8b层位清 晰，分层及连续性好，分辨率高。剖面电性结构大 致呈HK型，浅层存在一高阻层，在高阻层下为一 明显低阻薄层，连续性好，电阻率值约20 Ωm。深 部埋深在250550 m存在高阻厚层， 明显看出被多 条断层阻断，形成破碎带，断层附近有明显的低阻 异常体存在，电阻率为100200 Ωm，可解释为矿 体异常。 图 7 实测电场数据拟合对比 Fig.7 The fitting comparison chart of measured electric field data ChaoXing 第5期 唐传章等 基于边界约束有限内存的拟牛顿CSAMT一维反演及应用 199 图 8 实际资料一维反演结果对比 Fig.8 Comparison of one-dimensional inversion result of actual data 上述实际资料处理说明，LBFGS_B反演方法具 有较强的分层能力和对深部异常体的探测能力。可 靠性好，应用效果明显。 5 结 论 a. 改进的CSAMT正演算法，通过引入并行运 算，提高了反演速度与正演精度，为精确反演奠定 了基础。 b. 开 发 的 基 于 水 平 电 场 的CSAMT 一 维 LBFGS_B反演方法， 通过模型检验表明， 拟合度高， 反演稳定，迭代速度快，对薄层有较好的分辨能力。 c. 通过对两种不同目标函数的反演分析，发现 基于绝对误差的目标函数在反演时收敛速度快，调 用正演次数少， 相对基于相对误差的目标函数而言， 优势明显。 d. 实际CSAMT资料反演结果表明，LBFGS_B 反演方法在分层能力和对深部异常体的探测能力方 面优于连续介质反演方法。地质效果明显，具有良 好的应用前景。 参考文献 [1] LU Xinyu， MARTYN U， JOHN B. Rapid relaxation inversion of CSAMT data[J]. Geophysical Journal of the Royal Astronomical Society，1999，1382381–392. [2] 张凯飞. 基于电性约束的 NLCG 反演在 CSAMT 资料中的应 用[J]. 物探与化探，2016，403629–634. ZHANG Kaifei. The application of CSAMT data processed with NLCG inversion based on electrical constraints[J]. Geophysical and Geochemical Exploration，2016，403629–634. [3] 陈明生. 关于频率电磁测深几个问题的探讨一从可控源音 频大地电磁测深原理看解释中的问题[J]. 煤田地质与勘探， 2012，40563–66. CHEN Mingsheng. Problem of data interpretation from the view of principle of controlled-source audio-frequency magnetotellu- rics[J]. Coal Geology Exploration，2012，40563–66. [4] 严良俊，胡文宝，杨绍芳，等. 电磁勘探方法及其在南方碳酸 盐岩地区的应用[M]. 北京石油工业出版社，2001. [5] ROUTH P S， OLDENBURG D W. Inversion of controlled source audio-frequency magnetotellurics data for a horizontally layered earth[J]. Geophysics，1999，6461689–1697. [6] 王若， 王妙月. 一维全资料 CSAMT 反演[J]. 石油地球物理勘 探，2007，421107–114. WANG Ruo， WANG Miaoyue. One-dimensional full-data CSAMT inversion[J]. Oil Geophysical Prospecting， 2007， 421 107–114. [7] 蒋齐平. 基于场值数据的可控源音频大地电磁一维全资料反 演研究[J]. 科学技术与工程，2015，153193–196. JIANG Qiping. Research on 1D full CSAMT data based on filed data[J]. Science Technology and Engineering，2015，153 193–196. [8] 尚通晓，李桐林，关艺晓，等. CSAMT 一维全区反演[J]. 吉 林大学学报地球科学版，2007，37增刊 129–31. SHANG Tongxiao， LI Tonglin， GUAN Yixiao， et al. Full-region inversion of 1-D CSAMT[J]. Journal of Jilin UniversityEarth Science Edition，2007，37S129–31. [9] 朱威，李桐林，尚通晓，等. CSAMT 一维全区反演对比研究[J]. 吉林大学学报地球科学版，2008，38 增刊 112–14. ZHU Wei， LI Tonglin， SHANG Tongxiao， et al. Compared study of 1-D full-region inversion of CSAMT[J]. Journal of Jilin Uni- versityEarth Science Edition，2008，38S112–14. [10] 周军，袁伟，连晨龙. 可控源音频大地电磁法CSAMT一维 全区 Levenberg-Marquardt 反演方法[J]. 科学技术与工程， 2014，1414129–134. ZHOU Jun，YUAN Wei，LIAN Chenlong. Levenberg-Marquard inversion of CSAMT one-dimension region data[J]. Sci- ence Technology and Engineering，2014，1414129–134. [11] 汤井田，张林成，肖晓. 基于频点CSAMT 一维最小构造反演[J]. 物探化探计算技术，2011，336577–581. TANG Jingtian，ZHANG lincheng，XIAO Xiao. One-dimensional minimum structure inversion based on frequency CSAMT[J]. Computing Techniques for Geophysical and Geochemical Ex- ploration，2011，336577–581. [12] 刘颖， 李予国， 柳建新， 等. 海洋可控源电磁场的一维反演[J]. ChaoXing 200 煤田地质与勘探 第47卷 中国有色金属学报，2013，2392551–2556. LIU Ying，LI Yuguo，LIU Jianxin，et al. One-dimension inversion of marine controlled-source electromagnetic fields[J]. The Chinese Journal of Nonferrous Metal，2013， 2392551–2556. [13] 堃王鹏，曹辉，高妍，等. CSAMT 自适应正则化一维全资料 反演[J]. 物探化探计算技术，2013，355530–533. WANG Kunpeng，CAO Hui，GAO Yan，et al. Adaptive regu- larized inversion of 1-D full CSAMT data[J]. Computing Tech- niques for Geophysical and Geochemical Exploration，2013， 355530–533. [14] 李帝铨，王光杰，底青云，等. 基于遗传算法的 CSAMT 最 小构造反演[J]. 地球物理学报，2008，5141234–1245. LI Diquan，WANG Guangjie，DI Qingyun，et al. The ap- plication of genetic algorithm to CSAMT inversion for minimum structure[J]. Chinese Journal of Geophysics，2008， 5141234–1245. [15] LIU D C，NOCEDAL J. On the limited memory BFGS for large scale optimization[J]. Mathematical Programming Series B，1989，453503–528. [16] 孙彩堂，李玲，黄维宁，等. 基于自适应遗传算法的 CSAMT 一维反演[J]. 石油地球物理勘探，2017，522392–397. SUN Caitang，LI ling，HUANG Weining，et al. 1D inversion of CSAMT data based on adaptive genetic algorithm[J]. Oil Geo- physical Prospecting，2017，522392–397. [17] 贾定宇. 基于 L–BFGS 方法的水平电偶极一维反演[D]. 长 春吉林大学，2012. [18] 翁爱华，刘佳音，贾定宇，等. 有限长导线源频率测深有限内 存拟牛顿一维反演[J]. 吉林大学学报地球科学版，2017， 472597–605. WENG Aihua，LIU Jiayin，JIA Dingyu，et al. 1-D inversion for controlled source electromagnetic sounding using limited mem- ory quasi-Newton [J]. Journal of Jilin UniversityEarth Science Edition，2017，472597–605. [19] AVDEEV D， AVDEEVA A. 3D magnetotelluric inversion using a limited-memory quasi-Newton optimization[J]. Geophysics， 2009，74345–57. [20] BYRD R H，LU P，NOCEDAL J，et al. A limited memory algorithm for bound constrained optimization[J]. SIAM Journal on Scientific Computing，1995，1651190–1208. [21] 阮帅. 三维大地电磁有限内存拟牛顿反演[D]. 成都成都理 工大学，2015. [22] ZHU C， BYRD R H， LU P， et al. LBFGS-B Fortran subroutines for large-scale bound constrained optimization[J]. Report NAM-11，EECS Department，Northwestern University，1994. [23] MORALES J L，NOCEDAL J. Remark on “Algorithm 778 L-BFGS-BFortran subroutines for large-scale bound con- strained optimization”[J]. ACM Transactions on Mathematical SoftwareTOMS，2011，3817. [24] 陈小斌，赵国泽，汤吉，等. 大地电磁自适应正则化反演算法[J]. 地球物理学报，2005，4841005–1016. CHEN Xiaobin，ZHAO Guoze，TANG Ji，et al. The adaptive regularized inversion algorithm for magnetotelluric data[J]. Chi- nese Journal of Geophysics，2005，4841005–1016. [25] 严良俊，胡文宝，陈清礼，等. 长偏移距瞬变电磁测深的全区 视电阻率求取及快速反演方法[J]. 石油地球物理勘探，1999， 345532–538. YAN Liangjun， HU Wenbao， CHEN Qingli， et al. The estimation and fast inversion of all-time apparent resistivities in long-offset transient electromagnetic souding[J]. Oil Geophysical Prospect- ing，1999，345532–538. 责任编辑 聂爱兰 ChaoXing