当前位置:首页 >> 电力/水利 >>

二维大地电磁正则化共轭梯度法反演算法


 2007 年第 1 期             刘小军等 : 二维大地电磁正则化共轭梯度法反演算法

?71?

二维大地电磁正则化共轭梯度法反演算法
刘小军   王家林   吴健生 (同济大学   海洋地质国家重点实验室 ,上海  200092 )
摘    要 针对大地电磁二维反演中目标函数收敛速度慢而且解的稳定性

较差等问题 , 提出了大地电磁数据的正则化共轭梯 度法反演算法 ( Regularized Conjugate Gradient A lgorithm , RCGA ) 。此算法在构建目标函数时引入正则化的思想 ,利用 共轭梯度法求解最优化问题 。在每次迭代过程中根据目标的收敛情况更新正则化因子 , 有效地解决了迭代时目标 函数发散的问题 。最后对两个典型模型进行了试算 ,对其收敛速度及解的稳定性进行了分析 ,并与传统反演方法的 计算结果进行了对比 ,结果都表明 RCGA 具有收敛速度快 、 精度高 、 结果稳定等优点 。 关键词    大地电磁   目标函数   正则化因子   迭代   共轭梯度法

1      引 言 2  算法原理
实际情况下 ,大地电磁数据反演都是不适定的 , 通常解决这一问题的办法是引入吉洪诺夫的正则化 理论 ,其基本思想是在反演目标函数中加入一个约 束泛函来减小解的奇异性 。现在常用的约束稳定泛 函有 : 模型参数的范数 、 最大平滑稳定泛函 、 最小模 型泛函等 。其中比较著名的反演算法属 Constable
[1]

2. 1  反演目标函数的建立

反演的目的是寻找一个合理的地球物理模型 , 使其正演响应与观测结果相吻合 。但由于地球物理 反演问题固有的非唯一性 , 直接拟合数据目标函数 通常是不稳定的 ,而且不能得到真实的构造特征 ,给 解释带来困难 。为减少解的多解性 , 引入吉洪诺夫 等提出的正则化目标函数 : α P ( m ) = φ ( m ) + as (m )      ( 1 ) 式 中 φ ( m ) 为 数 据 拟 合 差 方 程 ,φ ( m ) =
A (m ) - d
2
D

开发的 OCCAM 法 , 它是寻找在极小可能构造下符 合数据的光滑模型 。在解目标函数的方法上 , 国内 外很多学者从提高反演速度 、 求最优化问题的方法 等方面提出了很多改进方法 , 如 Sm ith 等的快速松 弛法 ( RR I) , 其基本思想是假设电阻率主要以垂向 变化为主 ,这样大大加快了灵敏矩阵的计算速度 ; 国 内王小牧 、 于鹏等提出大地电磁与地震资料仿真退 火约束联合反演 , 此方法不用求目标函数的偏导数 以及解大型矩阵方程组 , 避免了线性化反演方法结 果依赖于初始模型的选取 ; 严良俊 ,胡文宝等将二次 函数全局非线性优化理论引入到大地电磁测深资料 的一维反演 , 防止了解陷入局部极小的问题 ; 陈小 斌
[2]

; m = m ( r) 为描述地下模型的标量

函数 , 本文算法中定义 m 为电层率 σ; d为实测数据 ; α为正则化因子 , 为一个正小数 ; s ( m ) 为稳定器 , 为 模型空间的目标函数 , 本文引入先验模型的最小二 2 范数作为模型目标函数 : s ( m ) = ,式 m - m ref 中 m ref 为参考模型 。 大地电磁测深数据与模型之间的关系是一种复 杂的非线性关系 , 其数学关系可表示为 : d = A ( m ) , 其中 A 为大地电磁正演算子 ; d为大地电磁资料反演 中的测 量 数 据 向 量〗 其 离 散 形 式 : d = { p ( x1 , , α φ φ w 1 )   ( x1 , w 1 )  p ( x1 , w 2 )   ( x1 , w 2 )  p ( xn , α α —————————————— — — — — — — —
基金项目 :中科院知识创新工程项目“ 环渤海 (湾 ) 地区前新 生代海相油气资源研究 ” 资助 。


提出了自适应正则化反演算法 , 实现了最平缓

模型约束下的大地电磁一维反演 , 在实际资料处理 中取得了较好的效果 。 本文在前人的反演算法基础上 , 根据共轭梯度 法稳定性高 ,计算机内存需求少的特点 ,将其应用到 大地电磁二维反演中 , 并在程序设计中引入选取正 则化因子的新方法 , 较好地解决了传统算法收敛速 度慢 、 稳定性差等缺点 ,并论述了算法的实现及应用 情况 。

收稿日期 : 2006 - 06 - 09 第一作者简介 :刘小军 , 1979 生 , 男 , 博士生 , 现主要 从事电磁勘探方法研究工作 。

?72? φ wm )   ( xn , wm ) } 式中 p 为视电阻率 ;φ为阻抗相 α 位 ; xi 为测点坐标 ; w i角速度 。 数据离散化后 , 反演目标函数 ( ( 1 ) 式 ) 可表示 为:
P ( m ) = (W d A ( m ) - W d d )
’ ’ ’

上海地质               总第 101 期   Shanghai Geology

α











其中 W m 为模型权重系数 、 d 为数据权重系数 , 为计 W 算简便 , 可令 W d 等于单位矩阵 I、 表示矩阵的转 3
2. 2  求目标函数极小的正则化共轭梯度法

置 。 ( 2 ) 式极小的最优化问题即可得到拟合结果 。 解 本文选用共轭梯度法求 ( 2 ) 式的极小值 , 迭代 在共轭梯度方向连续线性搜索 :
’ ’ ’ ’

过程如下 :

m n +1 = m n +δ = m n - kn l ( m n ) , m



其中 kn = ( l ( m n )


2 α m ) l ( m n ) ], m n +1 为第 n + 1 次 迭后 的模 型 值 , W







l ( m n ) 为第 n 次迭代后的共轭梯度向量 , l ( m n ) 为




第 n 次迭代后的梯度向量 , F 为灵敏度矩阵 , 3 表示 矩阵的转置 。 量相同 :



第 1 次迭代 ( n = 0 ) 的共轭梯度向量与梯度向


l (m 0 ) = l ( m 0 ) = F ( A ( m 0 ) - d ) +α m ( m n W
2

- m ref ) ) ,

第 n + 1次 ( n ≥ 1 ) 迭代共轭梯度向量为此次迭

代的梯度向量与前一次迭代的共轭梯度向量的线性 组合 :


l (m n +1 ) = l ( m n +1 ) +β +1 l ( m n ) , n

’ ’



其 中 : l (m n )


α( m - m ref ) 3 W m 2 ( m - m ref )

’ ’







= F (A (m n )
’ ’ ’ ’

m ref ) ) ; 系数 β +1 = n

度向量为共轭 。

重复上述计算 ( ( 3 ) ~ ( 5 ) 式 ) 可使目标函数收

敛到理想值 。 由于在每一次迭代中都采用抛物线函 数搜索方式 , 因此目标函数总能收敛到极小 。 在每一次迭代过程中都需正演计算理论观测响

应和观测数据关于模型参数的偏导数 , 这是反演中 最重要也是最耗时间的部分 。 为提高正演模拟精度 , 本文采用基于二次场的有限元法计算地表响应函

3

3 ~




l (m n ) ) / [ (m n )

l ( m n +1 ) l (m n )

~3


3











(W d A ( m ) - W d d ) +
’ ’

( 2)

( 3)













- d)
2

+ α m (m n W
2

2

, 它保证了两次梯

3 ~




l, ( F F +









3





( 4)

( 5) -



数。 将大地电磁场分为两部分 , 一部分为源在均匀背 景中产生的电磁场 ( Ep或 Hp ) ; 另一部分为将地下异 常体看作散射源 , 其产生的二次场 ( Es 或 Hs ) 。 由麦 克斯韦方程组导出二次场的偏微分方程 , 采用有限 [3 ] 元法 求解二次场值 , 再与解析法求出的一次场求 和后计算视电阻率与阻抗相位 。 由于一次场采用解 析法求得 , 相对减小了有限元法固有模拟误差 , 比总 场模拟效果更好 , 受网格密度影响相对较低 。 反演中求取灵敏度矩阵 F 的速度直接关系到反 演速度的快慢 。 本文根据电磁法中的互易定理 , 使用 [ 9, 10 ] 有限元法计算灵敏度矩阵 。 根据互易定理 源和 接收点可以互换 , 这样在测点放一点源做一次正演 计算 , 对每个单元进行积分就可以得到该测点平行 走向的场对所有单元电导率的偏导数 。 辅助场对单 元电导率的偏导数的基本原理相同 。 这种方法将正 演次数由模型变量个数减少到测点个数 , 大大提高 了计算速度 。 2. 3  正则化因子的选取方法 目标函数中正则化因子 α为数据目标函数与先 验函数之间的权重系数 , 它的大小决定了反演的拟 合对象 , 过大主要拟合先验模型 , 过小主要拟合实测 数据 。 因此 , 反演正则化问题的关键点是正则化因子 的选取方法 , 传统的方法是通过外部输入经验值 , 再 经过几次试算来确定其大小 , 这样增加了计算量 。 为提高收敛速度 , 有效地防止目标函数发散 , 在 [2 ] 编写代码时引入动态的正则化因子 。 其基本原理 是在每次迭代中 , 取初值 α ,α为此次迭代数据目标 0 0 函数 φn ( m ) 与模型目标函数 S n ( m ) 的初值的比 。 按 ( 3 ) 式更新正则化因子 α , 并计算更新 α 后的目标 k k 函数 φn ( m ) , 取 φn (m ) 最小时的 αk 定为一次迭代结 束。 这种自适应方法简单 , 保证了目标函数稳定收 敛 , 而且不需要输入任何经验值 , 不会带来明显的额 外计算量 。 ( 6) ak = α    k = 0, 1, 2, ……n; q > 0. 0

3  算例分析
为验证 RCGA 反演方法的正确性 , 利用所编程 序作了大量的理论模型试算 , 并与传统反演方法的 计算结果进行了对比分析 。 这些工作都表明了 RCGA 方法具有稳定性高 , 而且能更清晰地反映出电性分 布特征 。 限于篇幅 , 下面仅列举其中两个例子 。

 2007 年第 1 期             刘小军等 : 二维大地电磁正则化共轭梯度法反演算法

?73?

模型 Ⅰ为均匀半空间中两异常体 (图 1 中白色 Ω 虚线框 ) , 左边异常体电阻率为 10 ?m 、 边为 右 Ω? , 半空间电阻率为 50 ? , 异常体顶部埋深 Ω m 250 m 都为 500m 。 模拟测点数为 30个 , 采用 14个记录频点 ( 0. 01、. 03、. 1、. 3、、、、 、 、 、 、 、 0 0 0 1 3 5 10 15 30 50 100 300、 1000 ) , 正演数据中加入了 3 % 的高斯噪声 , 有 限元正演网格为 60 × (横向 × 30 纵向 ) 。 取先验模型

Ω m ref = 0, 迭代初始模型为电阻率 50 ?m 均匀半空 间 。 1 为 TE 模式迭代 20 次后的反演结果 , 从中可 图 以看出 , 两种方法都能较好地反映地下两个异常体 的位置和电阻率 , 但 RCGA 反演结果 (图 1 ( a ) ) 反映 异常体形态及大小都比传统反演方法的反演结果 (图 1 ( b) ) 更准确 , 基本没产生多余构造 。

   2 为模型 Ⅰ 图 计算中 , 采用两种方法反演时目 标函数的规一化值随迭代次数的变化情况 , 可以看 出两种方法迭代的目标函数值都是单调下降 , 但

RCGA 反演收敛速度要快 ,并且迭代 20 次后目标函

数值为 0. 016,传统反演方法为 0. 08,很明显 RCGA 反演效果更理想 。

   模型 Ⅱ 为一组合模型 , 如图 3 ( a ) 所示 , 基底为 Ω 一高阻隆起 ,电阻率为 500 ?m , 浅层为隆起后地 Ω 层变形 、 断裂为三个孤立体 ,电阻率分别为 1 ?m、 Ω 1 ?m、 Ω ?m。观测数据为 40 个测点 14 个频点 10 的视电阻率与阻抗相位 , 加入了 3% 的随机噪声 。 采用 矩形单元再剖分为四个三角形单元的有限元 正演计算网格 , Bostick一维反演 结果作为初始模 Ω Ω 型 ,反演中限定电阻率变化范围为 ( 0. 5 ?m , 600 ?m ) 。图 3 ( b ) 为 RCGA 反演结果 , 从中可以比较 准确地分辨出地质体的轮廓 , 高阻隆起界面也较清 晰地反映出来 ,目标函数迭代到稳定值的次数也不 太长 ; 而常规反演剖面图 3 ( c )只反映出了电性分界 面的大致形态 ,对浅部低阻体的勾画及各层的信息 反映很模糊 。结果说明 RCGA 反演方法比常规反 演方法要好 。
[ 11 ]

( a) 模型 Ⅱ 示意图

( b) RCGA 反演结果

?74?

上海地质               总第 101 期   Shanghai Geology

正则化因子的方法 , 较好地提高了解的稳定性和收 敛速度 。与采用传统反演方法进行模型试算的结果 进行对比分析 , 进一步表明了 RCGA 的高稳定性 、 高精度和较快的收敛速度 。 由于大地电磁正演计算速度的原因 , 作者未对 大地电磁三维资料的正则化共轭梯度反演方法进行
( c) 传统方法反演结果

研究 。
参考文献
[ 1 ]   C. deG root - Hed lin and S. Constable O ccam ’ inversion . s to genera te sm ooth, tw o ~ d im ensiona l m odels f rom m agnetotellu ric da ta. Geophysics, 1990, V ol 55 ( 12 ) : .

图 3   组合模型及两种方法反演结果
Fig 3   Subsurface combined model and comparison of the re2 . sults using two different method

4      结 语

1613 - 1624.
[ 2 ]  陈小斌 , 赵国泽等 . 大地电磁自适应正则化反演算法 .

对于数据量大 、 反演参数多的二维大地电磁数 据反演 ,选用正则化共轭梯度法反演电性参数 ,这一 思路与常规方法相比有稳定性高 、 收敛效果好及对 计算机内存要求不高等特点 。在正则化因子的选取 上抛开传统的取值方式 , 采取以目标函数收敛确定

地球物理学报 , 2005, V ol48 ( 4 ) : 937 - 945.
[ 3 ]  徐世浙 . 地球物理中的有限单元法 , 北京 : 科学出版

社 , 1994.
[ 4 ]  王家映 . 石油电法勘探 , 北京 :石油工业出版社 , 1993.

Inver sion A lgor ithm of 2 - D M a gne to te llur ic D a ta U sin g Regula r ized C on juga te Gra d ien t M e thod
L I X iao jun   U WAN G J ia lin   WU J ian sheng

( Sta te Key L abo ra to ry of M a rine Geo logy, Tongji U n ive rsity, Shangha i, 200092 )

A b stra c t: In th is p ap e r, a stab ilized inve rse m e thod of m agne to te llu ric da ta is p re sen ted, wh ich is nam ed R egu la rized Con juga te Grad ien t A lgo rithm , RCGA. The app roach is ba sed on regu la riza tion theo ry and adap tive regu la rized fac to r Con juga te Grad ien t A lgo rithm is u tilized to sea rch fo r the m in im um . The inve rsion re su lts of seve ra l syn the tic da ta a re ana lysed and comp a re w ith the re su lts u sing trad itiona l inve rsion m e thod. A ll of the se ind ica te the effec tivene ss of the p a ram e tric func tiona l . of RCGA.

Key word s:m agne to te llu ric, ob jec t func tion, regu la riza tion p a ram e te r, ite ra tion, con juga te grad ien t m e thod


相关文章:
电磁法
共轭梯度法反演、 拟线性近似反演、 快速松弛反演、...把大地电磁二维反演问题 转化为一系列一维反演问题。...无须直接解正则化方程,可避免偏导数矩阵 G 的计算...
文献检索作业 (2)
大地电磁三维正演,结合正则化反演方案和 共轭梯度...E Ex 广域电磁法已经实现了一维和二维正 反演。 ...反演部分为 John Booker 的 RRI 反演算法。Xinyou ...
共轭梯度法实验报告
下面给出共轭梯度法的具体计算过程: 给定初始向量 x...1 所张成的二维平面 ? 2 ? {x | x ? xk ?...继而重点给出正则化方法的实用共轭梯度算 法并举例...
地球物理反演理论
这种方法叫做正则化方法。 5.多解性 由于观测数据...二、从最速下降法、共轭梯度法、牛顿法、L-bfgs ...六、以地震勘探为例,任选取一种反演算法简述层析...
大地电磁法及其应用
方式来区分的, 并且这种区分一 般只在二维情况下才...搜索技巧来寻找合适的反演模型,如遗传算法、模拟退火...正则化反演方法介绍 3 大地电磁法及其应用 反演的非...
读书报告
大地电磁场亦受地形影响而发生畸变,常使反演结果出现...法、 正则化延拓求水平地 形场分布、 空间滤波法...作者还在共扼梯度算法的实现过程 中,用到了黄金分割...
大地电磁设计方案-王
计算, 研究所需要探测的主要电性标志层在大地电磁测...3)基于 TE 模式和 TM 模式的二维联合反演成像方法...本项目拟采用非线性共轭梯度反演方法,把大 地电磁...
大地电磁测深法
大地电磁测深法_冶金/矿山/地质_工程科技_专业资料...(spectre)分析, 获得各个场分量的频谱,然后计算它们...反演得到的是沿每个测量剖面的地下的二维电性结构(...
大地电磁资料静态校正方法研究
大地电磁测深法资料处理中的关键问题, 如果校正不当, 会使后续的反演解释得出...二维偏离度 0.6 0.4 -45 0.2 -90 103 10 2 101 100 10-1 10-2 ...
地球物理反演成像方法综述
地 球物理反演是地球物理资料定量解释的理论和算法...地震波与地球波速结构、大地电磁和重磁反问题中迅速...(Homotopy Method),非线性共轭梯度法(Non-Linear ...
更多相关标签:
共轭梯度 迭代正则化 | 共轭梯度法反演 一维 | 梯度正则化 | 梯度下降 正则化 | 共轭梯度法 | 共轭梯度算法 | 共轭梯度法 matlab | 共轭梯度 |