地球物理资料非线性反演方法讲座六共轭梯度法.pdf

第5 卷第4 期 2 0 0 8 年8 月 工程炒球物理荸赧 C H I N E S EJ O U R N A I O FE N G I N E E R I N GG E P H Y S l C S 文章编号 1 6 7 2 7 9 4 0 2 0 0 8 0 3 一0 3 8 1 一0 6 地球物理资料非线性反演方法讲座 T k 共轭梯度法 朱培民 王家映 中国地质大学地球物理与空间信息学院 V 0 1 5 N o 4 A u g 2 0 0 8 摘 要 概要地介绍了在非启发式非线性反演方法 共轭梯度法的原理 算法 优点 以及它的局限性 和 改进型算法 一种全局收敛的随机共轭梯度法 最后通过一个例子来说明共轭梯度法及其改进型算法的优缺 点 关键词 反演 共轭梯度法 随机共轭梯度法 中图分类号 P 6 3 1文献标识码 A 收稿日期 2 0 0 8 0 7 2 0 L e c t u r eo nN o n L i n e a rI n v e r s eM e t h o d si nG e o p h y s i c a lD a t a 6 C o n ju g a t eG r a d i e n tM e t h o d Z h uP e i m i n W a n gJi a y i n g I n s t i t u t eo fG e o p h y s i c sa n dG e o m a t i c s C h i n aU n i v e r s i t yo fG e o s c i e n c e s W u h a n 4 3 0 0 7 4 C h i n a A b s t r a c t T h i sp a p e rb r i e f l yi n t r o d u c e san o n h e u r i s t i ca n dn o n l i n e a ri n v e r s em e t h o d C o n ju g a t eG r a d i e n tM e t h o d i n c l u d i n gi t sp r i n c i p l e a l g o r i t h m a d v a n t a g e sa n dl i m i t a t i o n a l s oi n t r o d u c e sam e t h o d S t o c h a s t i cC o n ju g a t eG r a d i e n tM e t h o dw h i c hi Sa ni m p r o v e m e n to ft r a d i t i o n a lC G ag l o b a l l yc o n v e r g e n c ea n dh y b r i da l g o r i t h mw i t hs t o c h a s t i c h i l l c l i m b i n gt e c h n i q u e A tl a s t f ln u m e r i c a lt e s ti sg i v e nt os h o ws o m ea d v a n t a g e sa n dd i s a d v a n t a g e so fC o n ju g a t eG r a d i e n tM e t h o dc o m p a r e dw i t hS t o c h a s t i cC o n ju g a t eG r a d i e n t M e t h o d K e yw o r d s I n v e r s i o n C o n ju g a t eG r a d i e n t S t o c h a s t i cC o n ju g a t eG r a d i e n t 1 引言 从概率的观点来看 非线性反演方法可以分 成即统计方法和确定性方法两类 1 第一类又称 为启发式反演 如在本讲座的第二 三 四讲中 分 别介绍的蒙特卡洛法 模拟退火法 遗传算法等就 属于这一类 2 4 确定性方法为第二类 又称为 非启发式反演 主要是通过迭代过程来完成非线 性反演的 基金项目 国家自然科学基金 编号 4 0 1 7 4 0 3 3 6 0 4 7 2 0 6 2 资助 作者简介 朱培民 1 9 6 3 一 男 教授 博士生导师 主要研究方向包括地球物理反演 地震勘探和空间信息三维可视化等 E m a i l z h u p m c u g e d u c n 王家映 1 9 3 7 男 教授 博士生导师 主要研究方向为电磁法和地球物理反演理论 E m a i l j y w a n g c u g e d u c a 3 8 2 工程地球物理学报 C h i n e s eJ o u r n a lo fE n g i n e e r i n gG e o p h y s i c s 第5 卷 确定性非线性反演方法采用的反演策略是非 线性问题的线性化 主要利用目标函数的梯度信 息 通过反复迭代 寻找反演的最优解 梯度类方 法包括最速下降法 牛顿法 共轭梯度法 C o n j u g a t eG r a d i e n t C G 变尺度法等E 5 3 也包括纯粹 利用梯度信息的最小二乘 广义逆等变种的反演 方法 梯度类方法的优点是收敛速度快 缺点是 容易陷入局部极值 其解严重依赖于初始猜测 虽然统计类反演的最大优点是不完全依赖于 初始猜测 理论上在反演过程中不会陷入局部极 值 但是 这类方法的缺点也是十分明显的 就是 计算工作量巨大 效率很低 地球物理反演中涉 及的模型参数成百上千 就目前的计算条件来说 统计类方法仍然不能满足大规模地球物理反演的 要求 梯度类方法计算效率高 特别是共轭梯度法 经过不断的改进 并且和统计类反演方法结合形 成了统计加迭代的组合反演方法H 5 消除了依 赖于初始猜测的缺点 成了一种广受欢迎的反演 方案 本文将介绍最常用的梯度类方法 共轭 梯度法及其改进方法 2 共轭梯度法反演的基本原理 对于一个非线性反演问题 要寻找目标函数 的极值以获取对应的极小值解作为反演的结 果口 通常情况下 目标函数都可以近似表示成 一个多维的多项高次代数式 要求取这样目标函 数的极小值对应的解 可以从任意的初始猜测开 始 沿目标函数的负梯度方向 即最速下降方向 向前搜索即可快速到达目标函数的极小点处 沿 目标函数的负梯度方向搜索的方法 称为最速下 降法 S t e e p e s tD e s c e n t 理论上说 最速下降法 沿着负梯度方向搜索应能够很快找到极小解 但 是在实际的计算中发现 当迭代接近于极值点时 会增加成百上千次的搜索也难于到达极值点 这 主要缘于两个问题 在极值点附近 由于梯度值 逐渐减小 迭代步长越来越短 将迅速增加迭代次 数 在极值点附近 由于梯度值逐渐减小 很小 的数值误差 都有可能改变指向极值点的迭代方 向 从而增加迭代次数 这两个因素使得最速下 降法变得不具有实用性 有没有一种方案使得在有限迭代次数内 搜 索到目标函数的极值点呢 回答是肯定的 当目 标函数是咒维二次多项式时 如果迭代的搜索方 向沿共轭方向搜索时 能够在咒次迭代过程中 到达极值点 下面我们来介绍共轭梯度法的原理 及实现 1 引 2 1 最速下降法的搜索方向 设有如下的反演目标函数 妒 x 一百1x T H x b T X c 1 其中 表示模型参数 和b 均是九维向量 f 是 一个常数 H 是一个 z 行的正定对称矩阵 我们 要求取的解就是目标函数对应极值点的解 式 1 是一个标准的二次型 它对应的极值解就是 垆7 x V 甲 x H x 一西一0 2 即 H x b x 一Hb 我们假设在点 处开始沿负梯度方向r o r x 0 一妒7 x 搜索 到达点x 即 x 1 x o 卜a r o 为了使搜索能够快速到达极值点 我们这样 选取a 使甲 x 一甲 x 口r o 达到最小 即取盖P x 一0 或者 五d 9 x 1 d p X o G r o 一9 7 X l T r 0 2 0 上式说明 如果1 7 的取值 使得舻7 x 与r o 正交 即前一次搜索的方向必与下一次的搜索方 向正交 这样使得最速下降的搜索路径成空间锯 齿形 且下降方向是一组交替平行的梯度方向 图 1 这启示我们 虽然最速下降的搜索过程有成 百上千次 但实际上只需要一组咒个 咒维 彼此 正交的梯度方向就可以搜索到极值解 然而 在 设计实际算法时 从一组咒个 T 维 彼此正交的 梯度方向 直接计算出每一方向上的搜索步长是 相当困难的哺 但是 如果把正交搜索方向改为 共轭方向 算法就很容易实现 2 2 共轭方向 设P 和P z 为二个钾维向量 如果存在一个卵 2 的正定对称矩阵H 使得 P H p z 一0 3 就称P 和P 是关于矩阵H 共轭的 如果矩阵H 是一个单位阵 即H I 则 P j H p 2 一P P 2 一P j P z 一0 4 这就是我们熟知的正交向量的定义 所以 我们 可以把共轭看成是正交的推广 或者把正交看成 第4 期朱培民等 地球物理资料非线性反演方法讲座 3 8 3 是共轭的特殊情况 我们把彼此互相共轭的p 1 和p 所代表的咒维空间的两个方向称为共轭方 向 图1 二维空间情形下的最速下降法搜索路径 F i g 1 T h es e a r c hp a t ho fS t e e p e s tD e s c e n t u n d e r2 一Ds p a c e 2 3 共轭向量的构造 沿共轭方向的搜索 关键是要构造出咒个彼 此共轭的方向向量 为了说明 共轭梯度法的原 理 我们先考虑形式如式 1 表示的目标函数 设 有一组以维彼此关于咒 卯的正定对称矩阵H 共 轭的向量P P I I P 一 能够使我们分别沿着 这竹个共轭向量所指的方向各搜索一次 就可以 达到极值点X 设有一组线性无关的向量勖 g L L g 一 可以通过对它们的线性组合构造出一组7 个彼此 共轭的向量 不妨设P o g o 为此构造与P 共 轭的向量P P l g l 口op o 5 其中a 为常数 由于P 与P 关于H 共轭 的向量 故 p 丁 一面H 韪 a m 一一I I g l a o I g 一o 6 所以 Q 一躲 7 因而 旷旷蒜p 8 由于g o p 0 上式也可以写成 旷旷舞叭 在构造了共轭向量p o 与p 的基础上 设我 们已经求出了前k 个彼此共轭的向量p o p L L p k 一1 构造一个向量P k 与p P 1 L L P k 1 都H 共 轭 并先设 P 敫 口 p l 9 欲使p k 与P p L L p t 一 都H 共轭 便有 p H p O i 0 1 L L 忌一1 将式 9 代入上式 得 旷一揣 i 0 1 L L k 1 1 0 从而 旷数一蓦怒p i 依此类推 就可以求得所有一组彼此关于H 共轭 的所有向量P P L L P 一 2 4 共轭梯度反演方法 在上一节中 我们看到构造一组彼此共轭的 向量并不难 只需要找到一组线性无关的向量 g o g I L g 一i 即可 在2 1 节中 我们知道 最 速下降法在搜索的过程中 可以比较容易地计算 出一组n 个互相正交的梯度方向向量 我们不妨 取g o g L L g 一作为这一组行个互相正交的梯 度方向向量 来计算出共轭向量 为此 我们构造 沿共轭方向搜索的反演算法如下 1 给定初始搜索点x o 计算此点的梯度值 甲7 X o 一V 甲 0 一 罄 塞 k 裹 1 2 甚至 我们可以从上式构造一组正交向量 不失一 般性 取 勘一 塞 o L L o o 塞 L L o s n 一 2 o o L 塞 很显然 上述的g o 舀 L L 一 向量 是一组线性 无关的正交向量组 符合构造共轭向量组的条件 首先 令P 一g o 2 置循环变量 七一1 2 L L 行一1 按照公式 9 计算出共轭向量p 到P 一t 3 从初始搜索点x 开始 沿共轭向量p P 3 8 4 工程地球物理学报 C h i n e s eJ o u r n a lo fE n g i n e e r i n gG e o p h y s i c s 第5 卷 L L P 一 所指向的方向依次搜索 每步搜索的计 算公式是 X k X k 一1 做一1 P 一1 忌 1 2 L I 2 1 1 3 每次迭代的修改量为A X X I x 卜l 做一lP 一l 最终 经过n 步搜索 到达极值点 x 一x 一 届阢 1 4 为了求取仅 志一0 1 L L 7 1 我们根据 求极值的定义 在每次沿P t 搜索时 有 荔 x k 户品认 1 陬1 甲 广m 叫 0 根据式 2 有 7 x T p 一l 一 H x l b 7 P 一1 0 H 一l 像一1 P 一1 b 丁P k 一1 0 从而犀一等瓣 啬诂 1 5 如果把式 1 4 代入 2 可以证明X 就是目 标函数币 x 的极小值 同时也是方程H x b 的 解 在n 维目标函数为标准二次型的情况下 共 轭梯度法经过n 次迭代 能获取精确解的性质被 称为二次截止性 2 A 非线性共轭梯度反演 当目标函数是高于二次的连续函数 即目标 函数的梯度存在 时 其对应的解方程是非线性方 程 非线性问题的目标函数可能存在局部极值 并 且破坏了二次截止性 共轭梯度法需要在两个方 面加以改进后 仍然可以用于实际的反演计算 但 共轭梯度法不能确保收敛到全局极值 1 首先是共轭梯度法不能在n 维空间内依靠 咒步搜索到达极值点 需要重启共轭梯度法 继续 迭代 以完成搜索极值点的工作 2 在目标函数复杂 用式 1 0 计算a 时 由于 需要局部线性化 需计算海赛 H e s s i a n 矩阵H 且计算工作量比较大 矩阵H 也有可能是病态 的 一些研究者给出了不同的a 计算方案 但几 乎都抛弃了矩阵H 的计算 其中以F l e t c h e r 和 R e e v e s 1 9 6 4 的方案最为常用 具体形式如下 7 氆一一韭g l g 1 6 式中舀一 和数分别为第志一1 和看k 次搜索时计 算出来的目标函数的梯度 以F l e t c h e r 和R e e v e s 计算方案为基础的共 轭梯度反演算法如下 1 给定n 维目标函数妒 x 初值x x o 最大 迭代次数i 误差e 1 计算 9 0 一一V 甲 X o 置i 一0 k 2 0 9 2g o P 9 0 计算艿一 9 7 9 2 若i e 2 酊g o 执行 3 否则退出 并输出计算结果X 3 计算卢一黜 x x 胁 4 g 一V 中 x 乱d 艿一 艿 一9 1g a a i U n e 一 w p g a a p k 五 1 如果k n 或者g r p o 置p g k O 5 i i 1 返回第 2 步 3 共轭梯度法的局限性及其改 进 3 1 局限性 传统的共轭梯度法在地球物理反演和最优化 等方面 得到了广泛的应用 相对于其它的非线性 反演方法 共轭梯度法效率高 计算稳定 适用于 大规模的地球物理反演问题 从理论上来说 共 轭梯度法几近完美 但自从问世以来的几十年中 共轭梯度法仍然在不断改进之中 在非线性反演 中 共轭梯度法最大的局限是依赖于初值 不能确 保收敛到全局极值 为此一些学者不断探索 提出了两种解决方 案 一种是直接改进共轭梯度算法本身以达到有 条件全局收敛的目的 8 t 9 另一种是利用启发时 非线性反演的具有全局搜索全局极值的能力 加 上共轭梯度法的高效计算的混合算法 10 1 1 下边 我们要介绍的随机共轭梯度法就属于这一类 另外 传统的共轭梯度法也无法支持大规模 的并行反演计算 3 2随机共轭梯度法 随机共轭梯度法是将随机爬山法 1 2 1 S t o c h a s t i cH i l l c l i m b i n gM e t h o d 和共轭梯度法组合 在一起结合两类非线性反演方法的优点 以达到 既快速收敛 又不至于陷入局部极值的反演目的 试验证明这种算法是十分有效的 随机共轭梯度 法的反演步骤如下 1 根据先验信息确定初始模型 0 退出控制 误差e O 及其它控制参数 2 对任何一个模型参数的修改依据下式z P z s i g n d x i r a n d i 一1 2 L L 咒 第4 期朱培民等 地球物理资料非线性反演方法讲座 3 8 5 式中z P 和z 分别表示修改前后x 中的咒个模 型参数中的第i 个参数 d x j 表示第i 个模型参数 的基本偏离值 由人为决定 它的作用是控制模型 参数的修改幅度 s i g n 为符号函数 表示在相等的 概率下 即0 5 的概率 取值为1 和一1 r a n d 为取值范围在 o 1 的随机数 3 计算目标函数 p x 设修改前后的目标函 数分别为 P l d x 和甲 x 其差值为V 甲 牛 I d x 一妒 x 当 V 甲 O 取p O 4 当声一1 时 接受修改后的模型 开始共轭 梯度法搜索 共轭梯度法收敛后 保存此时的反演 解 5 反复2 3 和4 比较每一次接受的反演 解 直到收敛到允许的误差范围 3 3 随机共轭梯度法的并行化 在大规模的地球物理反演中 需要象随机共 轭梯度法这样既计算效率高 又能保证全局收敛 的算法 幸运的是 随机共轭梯度法继承了启发 式非线性反演的能够并行化的优点 合适的并行 化 可以在象微机机群 P c C l u s t e r 这样的计算 机上进行大规模和超大规模的地球物理反演 4反演例子 下边我们以一个二维情况的简单例子 来说 明共轭梯度法与随机共轭梯度法的反演效果 设 有目标函数 舳砧2 一去 e x p 一虹孚 叠 0 6 e x p 一虹孚望 1 7 式中 若取常数m 3 l t 一3 目标函数在点 P 1 3 0 和P 2 一3 0 处有两个极小值 分别是 0 6 0 1 0 5 8 和0 7 6 0 6 3 5 其中点P 1 处是全局极小 值 图2 表1 是传统共轭梯度法的反演试验结果 从 中可以看出反演结果严重依赖于初始猜测 当初 始猜测x X x 分别为 0 5 1 0 达到了全局 极值点 初始猜测为 一0 5 1 o 达到了局部极 值点 而初始猜测为 O 0 2 5 时 在有限的迭代 次数内 既没有搜索到全局极值 也没有搜索到局 部极值 这说明传统共轭梯度法是严重依赖于初 始猜测的 表2 是用随机共轭梯度法 在三次不同的时 间反演试验的结果 从表中可以看出 虽然三次 反演的迭代次数不同 但均得到了相同的反演结 果 都搜索到了全局极值解 这说明随机共轭梯 度法是具有全局收敛能力的反演算法 图2共轭梯度法与随机共轭梯度法反演例子的目标函数 F i g 2O b j e c t i v ef u n c t i o ng i v e ni nn u m e r i c a lt e s t sf o r s h o w i n gC o n ju g a t eG r a d i e n ti n v e r s i o na n d S t o c h a s t i cC o n j u g a t eG r a d i e n ti n v e r s i o n 表1 T a b l e1 共轭梯度法的反演试验结果 R e s u l t sO fi n v e r s i o nt e s t su s i n g C o n j u g a t eG r a d i e n t 注 表中的迭代次数为认为限定的最代迭代次数 表2随机共轭梯度法三次反演试验结果 T a b l e1R e s u l t so ft h r e ei n v e r s i o nt e s t su s i n gS t o c h a s t i c C o n j u g a t eG r a d i e n t 初始猜测反演结果 目标函数值迭代次数 z 1 z 2 z 1 z 2 3 8 6工程地球物理学报 C h i n e s eJ o u r n a lo fE n g i n e e r i n gG e o p h y s i c s 第5 卷 参考文献 1 1 3 王家映 地球物理反演问题概述 J 工程地球物理 学报 2 0 0 7 4 1 1 3 2 王家映 蒙特卡洛法 J 工程地球物理学报 2 0 0 7 4 2 8 1 8 5 3 1 师学明 王家映 模拟退火法I j 工程地球物理学 报 2 0 0 7 4 3 8 1 8 5 1 4 师学明 王家映 遗传算法 J 工程地球物理学报 2 0 0 7 4 3 1 6 5 1 7 4 5 1 席少霖 赵风治 最优化计算方法 J 上海 上海科 学技术出版社 1 9 8 3 9 5 1 2 6 63S h e w c h u kJR A nI n t r o d u c t i o nt ot h eC o n j u g a t e G r a d i e n tM e t h o dW i t h o u tt h eA g o n i z i n gP a i n J C a r n e g i eM e l l o nU n i v e r s i t y P i t t s b u r g h 19 9 4 7 3F l e t c h e rR a n dR e e v e sC F u n c t i o nm i n i m i z a t i o nb y c o n j u g a t eg r a d i e n t s J C o m p u t J 1 9 6 4 71 4 9 1 5 4 1 8 T o u a t i A h m e dD S t o r e yC G l o b a l l yc o n v e r g e n t h y b r i dc o n j u g a t eg r a d i e n tm e t h o d s J J o u r n a lo f O p t i m i z a t i o nT h e o r ya n dA p p l i c a t i o n s 1 9 9 0 6 4 2 3 7 9 3 9 7 9 S h iZa n dG u o J An e wa l g o r i t h mo fn o n l i n e a rc o n j u g a t eg r a d i e n tm e t h o dw i t hs t r o n gc o n v e r g e n c e J C o m p u t a t i o n a l 而且避免了直接求G以及GTG的逆矩阵 也避免了存储G和GTG所需庞大的存储量 另外 由于反演 参数太多 求模型光滑约束的最小构造反演能够有效地消除多余构造信息 得到可靠的反演结果 将这3种方法和技术融合于三维反演中 取得了好的反演效 果 为改进传统最小构造反演收敛慢的问题 还提出了一种新的反演迭代技术 仅需10次左右甚至更少迭代即可收敛 2 学位论文 朱培民 用随机共轭梯度法反演裂隙密度参数的研究 1999 该论文主要是研究如何用纵波AVO数据反演提取裂隙储层参数的问题 为此作者在这一领域开展了一系列的研究工作 论文的研究工作主要集中在以下 三个方面 首先研究了用纵波AVO数据反演裂隙密度参数的可行性 并推导了纵波AVO反演的基本公式 然后提出了一种新的非线性反演算法 随机共轭梯 度法 最后设计并编程实现了随机共轭梯度反演裂隙密度参数的算法 还就多个理论模型进行了试算 3 期刊论文 徐凯军 李桐林 张辉 李建平 XU Kai jun LI Tong lin ZHANG Hui LI Jian ping 基于共轭梯度法的 垂直有限线源三维电阻率反演 煤田地质与勘探2006 34 3 利用垂直有限线源研究油田注水分布和剩余油分布成为一种新型电测方法 本文对该方法实际应用的有效性进行了分析和说明 并利用共轭梯度迭代 技术实现了垂直有限线源三维电阻率反演 结果表明 反演稳定可靠 计算速度快 4 会议论文 刘波 卢涛 姜培学 基于共轭梯度法的三维弯管内壁温度的反演 2009 用外壁测量温度与计算温度构建反问题目标函数 以有限单元法作为导热正问题求解方法 基于共轭梯度算法建立了三维弯管非稳态的导热反问题 模型 利用弯管外壁的测量温度反演出内壁的温度 并与内壁温度实验值对比 对比结果显示 反演值与实验值符合较好 利用该模型可以精确地反演 出弯管内壁面上各点的温度波动值 5 学位论文 李永杰 反射旅行时的自动拾取与共轭梯度法旅行时反演 1996 6 期刊论文 刘喜武 刘洪 李幼铭 LIU Xi wu LIU Hong LI You ming 反射系数与子波同时迭代反演的预条件共轭 梯度法 物探化探计算技术2006 28 3 在反射系数白噪 子波最小相位的假设下 研究基于线性反演的地震反射系数和子波同时估计问题 在Cauchy准则稀疏反演求解中 应用顸条件共轭梯 度法实现反射系数和子波同时迭代反演 在迭代求解正则化方程时 用共轭梯度法求解相应的原问题 初猜子波求解也使用该策略 模型数据试算与比较 表 明了该算法正确而有效 用实际数据检验算法的实用性 经研究表明 预条件共轭梯度法计算的反射系数和子波 要比直接稀疏反演精度高 而且收敛较快 数值稳定 实用性强 7 学位论文 胡厚继 混合遗传算法在大地电磁反演中的应用 2006 遗传算法是近些年来产生和发展的一种模拟生物进化过程的自适应启发式全局优化的搜索算法 它不完全依赖于初始猜测 且具有全局收敛的特点 可以被用来解决各种复杂的实际问题 如工程优化设计 人工智能和决策系统以及地球物理反演等 尽管遗传算法是一种效率很高的全局优化算法 但许多仿真结果表明它有计算时间长 局部搜索能力弱的缺点 而共轭梯度法属于非启发式全局优化搜索方法 收敛速度快 但容易陷入局部极值 且严重依赖初始猜测 根据遗传算法和共轭梯度法的特点 本文提出了一种混合遗传算法 用来进行地球物理反演 该算法既具有遗传算法的全局收 敛性 又有共轭梯度法的快速收敛性 文章首先对大地电磁理论做了详细的介绍 并综述了目前几种应用比较广泛的一维反演方法及其存在的问题 说明还需要对大地电磁一维反演方面 做进一步的完善工作 在大量阅读国内外相关文件的基础上 对遗传算法做了比较详细的文献综述 接下来对遗传算法的工作原理及其基本术语做了比 较详尽的