全波波形处理中的奇异值分解方法.pdf

24 测井与射孔 1999 年 第 2 期 全波波形处理中的奇异值分解方法 点 至 (江汲石 油学 院 ) I) 摘要在 声 波 测 井垒波 波 形 处 理 的很 多方 法 中 都 会 碰到 求 解 线性 方 程 组 的 问题 ， 井 且往 往都 是复 超 定线性 方 程 组 ， 而 且病 态 程 度较 高 ， 对 于一 般 的 解方 程方 法 大 多失 效 从 而 使 这些 数据 处理 方法 的运 用 受到 了限制 本 文 针 对过 一 问题 ， 介 绍 了用 奇 异 值 分 解 法 求 解病 态 超 定 线 性方 程 组 的 理论 基 础 和 计 算 方 法 +井 针对 复 型情 况 给 出了一 种 解 j兜的变通 方 法 。该方 法褥 出 的解 稳 定性 好 ， 精确 度高 ， 用 于垒 嫂 嫂形 处理 具 有分 辨 率高 、 容 噪 能力 强等 优 点 可 提 高解 释精 度 主题 词垒 波 嫒形 处理 奇 异值 分 解 广义逆 最 小二 乘 前 言 ； 皮 f 声波 全波测井所得到 的声 波全波波形中包含有 纵波、 横波 、 斯通利 波和 伪瑞利波等 各组分 波 ， 如何精确提取 各组 分渡 的信 息是 全波波形处理中 的关键 问题 。 目前 的全波渡形处理方法 主 要有直接相 位法 (D P D )、 时差 时间相关法 (ST C )和 协方差 法等 。在这些方法 中 都要 经常 涉及 到求解 下列线性 方程 组的问题 ： A X B (1 ) 式 中 A为 m n 阶系数矩 阵； x 是 n 1 阶矩 阵 为待求的未 知 向量 ； B 为 m l阶矩 阵 ， 一 般 是与实 际观 测值有关 的向量 。 在 (1 )式 的线性方 程组 中 如果 其系 数矩 阵 A为对称正 定 的非 稀疏矩 阵 ， 则 用 C hol esky 分 解求 解是最快 和最精 确 的算 法 ； 如果 A 为非奇 异矩阵而且 不是方 阵 用 Q R 分 解是 目前 最 安全 的 。但在地球物 理反演中涉及 的系数矩 阵 A 多半是接近奇 异的 因为这时 我们要 用观 测 的信号 去推测场 源或 传播介质的参数 ， 而从信号 向源推移的过 程具有解 析上 的不 稳定性 。 作 为 谱展开 的一种形 式 奇 异值分解 方法在地球物理反演 中有 着重 要意 义。 如在 全波波形 处理的协方差方法 中 ， 应 用快速傅 氏变换将 波形 上所有数据 点变换到频域 ， 然后对 每个频率 处 ， 运 用扩 充的 P rony 法 求取幅度 、 波 数、 衰减和相位 时 就 会涉及到求解 复 型超定 线性 方程 组 。对于超定线性方程组 ， 一般采用 最小 二乘 法求解 ， 即求以下方程 的解 A A XA B (2 ) X一(A A)- LA B (3 ) 其 中 T 表示 转置 ， 一l表示矩 阵求逆 。(A A ) A 称 为矩 阵 A 的 广义逆 。当矩 阵 A 的条件较 好 时 ， 用一般 的方法 由(3)式直接计算 是可 行的。 但在 全波 渡形 处理中 ， 由于参数之 间有 很强的 相关性 ， 在观 测数据中带来的参数信息不足， 以及不可避免的观测误差 ， 在数值上则表现为 A 收 稿 日期 l 1998一l O一07 维普资讯 学兔兔 w w w .x u e t u t u .c o m 第 2 期 高军 ， 等： 全波 波形处理 中的奇异值分解 方法 25 矩 阵某些 奇异值很 小 ， 从而 使反演精度 降低 ， 这就是 矩 阵理 论中所谓 的病 态方程组 。 求 解病 态方程组 ， 用一般的方法 ， 如主元素 消去法 、 I u 分解法 、 H ousehol d 法 以及迭 代法 ， 包括一些 最优化方法 ， 都 不能得到稳定合理 的解。本 文为解决这一 问题 ， 采用了奇异值分解法 求取 广义逆 ， 然后 求取鼎 小二乘 解 的方法 特别地 ， 针 对全波波形处理 中的复型情况 ， 给 出了一 种 变通 的方法 ， 并编 制 了具有 实 用意 义而 又 比较 优化 的程 序 ， 实际运 算表 明 ， 奇异 值分 解方法 具 有克服病 态能力强 、 运 算中不放大误差 的优点 ， 求得 的解稳定性好 、 精确 度高 将 其应用于全 波波形 处理 具有分 辨率高 、 容 噪能 力强等 优点 ， 从 而可 以获得 精确 的 各组分 波 信息 ， 提高解释 精度 。 一 般实 矩阵的奇异值 分解 奇 异 值分 解的 萌芽 思 想至少 可 以追 溯到上 个世 纪 ， 如 Syl vester 于 1899 年 (M esseger of M ath ， 1899 ， V o1 42)发表 的文章 。但作 为数值 计算 的工具被 人们所重 视还是 1965 年 以后 的 事 。G ol ub 和 K ahan (1 965 )首 次指 出了直接 由方程 系数矩 阵计算奇异值 存在的 问题 ， 1970 年 G ol ub 和 R ei nsc h 提 出了下 面要 讨论的奇异值分 解算法 。 理 论基础 方程 (1)中 ， A 矩 阵可以是 方阵 (m = n ) ， 也 可以是长方 阵 (m n )， 不妨 设 m n ， 更具一般 性。理论上讲 ， 解方程组(1)的方法有多种， 但在原始数据不准确 ， 算法本身不精确的计算实践 中， 奇异值分 解方法基 本上是 已经 知道 的唯一 可靠 的方法 。如果 A 矩 阵是满 秩 的， 则解 唯一 ， 并 能用几个不 同的方法 可靠地求 解 ， 其 中某些方法 比奇异值分解法快 ， 但是奇 异值分解法还可 以对 付 A 矩 阵不 满秩 的情 况 ， 除某些大 型 问题 (A 矩 阵的 阶次 特别高 )之外 ， 并不 比其 它方法 多花更多时 间 ， 实际上 ， 从某种意 义上讲 ， 它 比给 出错误解的快速方法 还要 省时 。 对任意矩 阵 A (方 阵或 长方 阵 )进行 奇 异值分 解 ， 是 因 为奇异值 具有 稳 定 的性 质， 可以证 明 ， 当矩阵 A 有一个 扰动 E (波动矩 阵)时 ， 奇异值 的变化 不会超过 E 矩 阵的谱 范数 l l E l l ， 即 E 的最大本征 值 。奇异值分 解的稳定性 ， 对 用广 义逆方法解线性方程 组是 有意 义的。 首先 ， 给 出如下 的实矩 阵奇异值分 解定理 ： 定理 l 癌 A 为 m X n 阶实矩 阵， 则 存在 一个 m X 111 阶的列正交 矩阵 u 和 n X n 阶 的列正 交 矩阵 V ， 使 rn A U l IV 0 (4) L 0 0J 成 立 。 其 中， Y = di ag (6t 52， ， 6T)， r m i n m n )， 且 6l 62 q 0 。 上式称 为实 矩阵 A 的 奇 异值分懈式 ， (i = 1 2 ， ， r)称为 A 的奇异值 。 利用 A 的奇异 值分 解式 ， 可 以计 算 A 的 广义逆 A 。(5) 设 U 一 (U 1， U ： )， 其 中 u l为 u中前 r 列列 正 交 向量 组构成 的 m X r 阶矩 阵 ； V 一 (V ， V 。 )， 其 中 V 1 为 V 中前 r 列列正交 向量组构成 的 n X r 阶矩 阵。则 A 的广义逆 为 A = V 一 U T (5) 利用 A 的广 义逆可 以求解如 (1)式 的线性最 小二乘问题 ： A X B 其 中 A 为 m n 阶矩阵 ， x 为 n 维列 向量 B 为 111 维常数 列向量 。该问题 的最小二乘解 为 i j 维普资讯 学兔兔 w w w .x u e t u t u .c o m 26 测井与射孔 1999 XA B V L一 U TB (6 ) 因此 ， 求 解 问题 归结 为广 义逆矩 阵 A 的求法 ， 具 体说 来 就 是用奇 异值 分解方 法 来求 广 义逆 矩 阵 A 。 计 算方 法 奇 异值分解 的计 算过程分两大步 。第 一步 ， 用豪斯荷 尔簿(H ousehol der)变换 将 A 约化 为 双 对角线矩 阵。即 B O A = O S P一1 e p s p (7 ) 其 中 OU LU 2 一 U k，k =m i n n ， m一 1) 一 v Lv 2 v ，s m i n m ， n 一 2 ) (8 ) O中的每 一个 变换 U (i - 1， 2 ， ， k )将 A 中第 j列主对角 线 以下的元 素变 为零 而 中的每 一 个 向量 V (j一1， 2 ， ， s)将 A 中第 j行 中与主对角 线紧邻的 右次对 角线元素 右边 的元素变 为零 。 对于每一个变换 v j具有如下形式 ： I pV 。 v 7 (9) 其 中 I 为单位 矩 阵； P 为一 比例 因子 ， 以避免计算 过程 中的溢 出现象与误差 的积 累。v j 是一 个 列 向量 V 】 = (v L， 2， v ) (10 ) 则 A V ， = a pA V ，V = A w v i (11) 其 中 W = pAY， = P(v|rd_ ， VIa -， v ) (12) 第二步 ， 用变形 的 Q R 算法进行迭代 ， 计算所有 的奇 异值 。 由第一步 已经得 到了一个双对角线矩 阵 B 。 在这一步 中， 用一 系列 的平 面旋 转变 换将 B 逐 步变成对 角矩 阵。在 每一次的迭代中 ， 用变换 B = U 【 pu u B V L 2V 2r “V 一 】 (13 ) 其 中变 换将 B 中第 j 列 主对角线下 的 一个非零元 素变 为零 ， 同时在第 j 行的 次对角 线元 素 的 右边 出现 一个非零 元素 ； 而变换 V +。 将 第 j l行的次对角 线元素右边 的一 个非零元 素变 为 零 ， 同时在第 j列 的主 对角线元 素的下方 出现一个非零 元素 。 由此 可知 ， 经过 一次迭代 (j= l 2 ， ， p 一1)后 ， B 仍 为双对角线 矩 阵。 但随着迭代 的进行 ， 最后收敛 为对角矩 阵 其 对角线上的元 素 即为奇 异值 。在每 次迭代 时 ， 其初 始变换 V ： 后 ， 将在第 l列 的主对角线下方 出现一个 非零 元素 。在变换 V 中 ， 选 择位移值 的计 算式 如下 ： b =(s 一 】 +Sp) (s 一 】 一 Sp) +e； 一 L 2 维普资讯 学兔兔 w w w .x u e t u t u .c o m 第 2 期 高军 等 ： 全 波波 形处理 中的 奇异 值分解方 法 27 c 一(s e 1) d sg n (b ) b + C =S：一 c (b +d ) 最 后需要对奇异值按 非递 增次序进行 排序 在上 述变换过程 中 ， 若对于 某个 次对角线元 素 ej 满 足 I e I (I S， +。 I +I S， I) 则 可认 为 e； 为零 ， 为预先给定 的精度要 求 。 若对 角线元素 S 满 足 I S I e( Ie I +Ie， I) 则可认 为 S 为零 (即为零奇 异值 )。 复超定线性方程组的最 pJ,Z 乘解 (14 ) (15 ) (16 ) 在声波 测井的全波 波形 处理 中 ， 由于经常要涉 及到 傅立 叶变换 (1)式 中的矩 阵多 为复 型 矩 阵 ， 因此 ， 有 必要在实矩 阵的奇 异值分解 的基础上进 一步研 究复矩 阵的奇异 值分解 ， 以求 得 复超 定线性方程组 的最小二乘解 。 理论 基础 同样 ， 我们先给 出如下定理 ： 定理 2 ： 令 A 是一个 m X n 维复 数矩 阵 ， 则分别 存在 一 个 m X m维和一 个 n n 维 酉阵 u 和 V 使 得 A U 三V (17 ) 其 中上 标 XI表 示矩 阵的共轭转 置 ， 是一个 m n 维 对角 阵 其 主对角 线上 的元 素是 非负 的 并按 下列顺 序排 列 ： IId2 2 O (18 ) 式 中 h m l n m ， n ) 所谓酉 阵 ， 是 指如矩 阵 U 满 足 U = U ， 则 称矩 阵 U 为酉 阵。特 别地 ， 一个实 的酉阵 U 又 叫正交 阵 即有 U 一U 。对角元素 称为矩 阵 A 的奇异值 。 同实矩 阵情形类 似 对复矩 阵情形下 的线性方程 组 (1)式 仍可 以先求对复 矩 阵 A 进行 奇 异值 分解 ， 然后 由(5)式求取 A 的 广义逆 A ， 最 后 由(6)式求 出(1)式 的最 J- 乘解 。 但 这样 计 算 比较烦琐 ， 目前一般 的信 号处 理或矩 阵计算程序库均未包含复矩 阵的奇异值分 解子 程序 ， 因 此 ， 本 文在实矩 阵的奇异值 分解基础上 ， 应用一种变通 的方法 ， 可 以避 免复 矩阵奇异 值分解 的 烦琐 过程 ， 很简单得求取 (1 )的最 小二乘解 。 计算方法 对 于 线 性 方 程组 A X B 其中 A 为 m X 11维 复数矩 阵 ， B 为 m X 1 维复数矩阵 ， x 为 TLX 1 维复数矩 阵。不妨 设 A A t(BR 、 B t X R、 X )分 别表 示矩阵 A (B X )中 各元 素 的实 部和虚部构 成 的矩 阵 ， 显 然 ， A 、 A 为 m X n 维 实矩阵 B 、 B -为 m X 1 维实 矩阵 ， X 、 X -为 111 维 实矩 阵。则上 式可写为 维普资讯 学兔兔 w w w .x u e t u t u .c o m 28 测井与射孔 1999 皋 (A R + i A 1)(X R + i X 】 ) = B R + i B 展 开 为 fA RX R A l X 】 一 B R 【 A RX l A l X R Bl 即为二维复矩 阵方程组 对此方 程组 又可展开如 下 ； 对 于 式 (20 1)， 有 同样 对 于式 (20 2 ) 也有 将两 式合并 ， 即有 A R 一A A R 一A 二 XR】 ： X R X l 1 ： X l n X Rl ： X R X l 】 ： X X R1 ： X R X l 】 ： X In (19 ) (20 1 ) (20 2 ) (2 1 ) (2 2 ) (2 3 ) 即复矩 阵情形下 的超 定线性方 程组 的求解可转化 为如下 实矩 阵超定线性 方程组 的求解 ： C Y D 其 中 ， C 为(2m )(2n )阶实矩 阵 Y为(2n )l阶实矩 阵 ， D 为 (2m )1 阶实矩 阵 c 一 Y 豳 n 一 = 翻 应用 实矩 阵奇异值分 解方 法可 以很 容易地 由式 (24)解 出 Y ， 即可 由(26 )式得 出 X 。 (24 ) (25 ) ( 26 ) ( 27 ) 踟 ； ； ； A A 一 一 一 一 ： 主 ： 一 一 一 一 一 ； ； 一 一 一 一 m n n 耋： 主1： 一 主|： 三 维普资讯 学兔兔 w w w .x u e t u t u .c o m 第 2 期 高军 ， 等 ： 全波波 彤处理 中的奇 异值分解 方法 z9 实 例 协方差方 法是 现代信 号处理应用 于全渡翊l 井解释 的一种 新方法 ， 可 较好地 提取 各组分波 信息 ， 得 到 比较准确 的 各组分 波时差 、 衰 减和 幅度值 。我 们将奇异值分解法应 用于协方差方法 中求解复 型超 定线性 方程组 ， 对塔 里术的塔 中 24 井和塔 中 103 井 的 D SI(偶极 横波成象 测井 ) 声披测井资料进行了处理 ， 得到 了比较理想的纵波、 横波和斯通利波速度和衰 减， 与其它测井 资料所反 映的地层特 征吻合很好 ， 从而说 明了奇异值 分解法所得 解的稳定性和精确性 。 结 论 奇异 值分解法 (SV D )作为一种 信号处理技术 ， 主要 用于求解 线性 方程组 。实例证 明 ， 奇 异 值分 解方 法具有克 服病态能力 强、 运算 中不放 大误差 的优点 ， 求得 的解稳定性好 、 精 确度高 。 将 其应用于全波渡形处理具有分辨率高、 容噪能力强等优点 ， 从而可以获得精确的各组分波信 息 ， 提高解 释精度 。 参考文献 K J E JLe|sen and C H C h en g an d K M T u bm an ，1988 E sd m ati ag phase veLoc i ty and attenuad o n 0f gui d ed s i n ac ou sti c l ogg i n g data ，G eo ph ysic s ， V oL-54 n o 8 p p- 10 5410 59 何