稳健稀疏反褶积方法研究

况下 ，通 过地 震波 记录 反求 解 反射 系数 与人 工地 震 子波 的方 法称 为 反褶积 方 法 ． 地 震 波记 录 的过 程 一 般使 用 数 学上 的褶积 描 述［ １ 】 ．如 果 假 设 地 下 介 质 为 古 皮 奥 的 水 平 层 状 介 质 模 型 ， 而 人 工 产 生 的 地 震 子 波 为 雷 克 子 波 ， 则 地 震 波 记 录 可 以看 做 是 由人工 产 生的地 震子 波 与地 下 的波 阻抗变 化 的界 面对 应 的反射 系 数 的褶积 过程 ．
得 到对 于 Ｘ的估 计 ，并构 造地 震波 褶积 矩 阵 Ｘ． 步骤 ３ 重复 步骤 １ 和 步骤 ２直到 算法 收敛 ，终止 ．
（ １ ５ ）
可以看出，该算法 的主要计算量集中在求解优化 问题 （ １ ４ ） 和（ １ ５ ） ．下面我们将说明如
令ｙ （ ｔ ） 为检 波 器记 录 的地 震波 记 录 ，ｘ （ ｔ ） 为人工 产 生的地 震子 波 ，ｒ （ ｔ ） 为反射 系数 ， 则ｙ （ ｔ ） ＝ｘ （ ｔ ） 木 ｒ （ ｔ ） ＋ｎ （ ｔ ） ，其 中ｎ （ ｔ ） 是地震波记录采集时所遇到的噪声，符号 半 表示两
个 向量之 间的褶积 ．
大化 ｒ得到 如 下优化 问题
ｉ ｎ
（ ７ ）
（ ８ ）
所以，根据极大化后验概率准则可以去估计反射系数．因此，对式 （ ７ ） 两边取对数 ，并最
ｉ ｎ Ｊ Ｉ Ｙ—Ｘｒ ｌ ｌ ； ＋刈ｒ
、
其中 ｆ ｊ ． １ Ｊ 表示向量的欧式范数．问 题（ ８ ） 就是文献 ［ ６ 】 定义的贝叶斯反褶积的数学框架． ３ 稳健 稀疏 反褶 积方 法
步骤 如下 ： 步骤 １ 固定地 震波褶 积 矩 阵 Ｘ ，求解
ｍ ｉ ｎ Ｉ Ｉ ｖ —ｘｒ Ｊ Ｊ ｌ ＋￣ Ｊ Ｊ ｒ Ｉ ｌ １ ．
步骤 ２ 得 到反射 系数 的估 计 后 ，构造 反射 系数 褶积 矩 阵 Ｒ ，并求解
（ １ ４ ）
ａ ｒ ｉ ｎ ｆ ｒ Ｙ —ＲＸ ｆ ｌ ｌ ，
第 ３ ２ 卷 第３ 期
２ ０ １ ５ 年０ ６ 月
工
程
数
学
学
报
Ｖ ｏ １ ． ３ ２ Ｎ 。 ． ３
Ｊ ｕ ｎ ｅ ２ ０ １ ５
ＣＨＩ ＮＥＳ Ｅ Ｊ ＯＵＲＮＡＬ ＯＦ ＥＮ ＧＩ ＮＥＥＲＩ ＮＧ Ｍ ＡＴＨＥＭ ＡＴＩ ＣＳ
ｄ ｏ ｉ ： １ ０ ． ３ ９ ６ ９ ￣． ｉ ｓ ｓ ｎ ． １ ０ ０ ５ — ３ ０ ８ ５ ． ２ ０ １ ５ ． ０ ３ ． ０ ０ ２
（ ４ ）
（ ５ ）
即反 射系 数 为高斯 白噪声序 列 ． 同时 ，他还 假 设采 样 的噪音 ｎ为独立 的高斯 白噪 声 ，所 以
） ＝ 南 ｅ ｘ ｐ ｛ 一 雩） ，
即
பைடு நூலகம்
（ Ｙ Ｉ ｒ ） ＝ 面 ｅ Ｘ ｐ ｛ 一 一 （ Ｙ － Ｘ ｒ ） Ｔ （ Ｙ － Ｘ ｒ ） ）
收稿 日期： ２ ０ １ ４ － ０ ４ － ２ ４ ． 作者简介 ：高倩 （ １ ９ ８ ４ 年１ Ｏ 月生 ） ，女 ，博士 ． 研究方 向：地震信 号处理与地理信息系统 基金项 目： 国家 自 然科 学基金 （ １ １ ４ ０ １ ４ ６ ２ ） ．
３ ３ ０
工
程
数
学
学
根据 贝叶斯 公式 ，我 们有 反射 系数 ｒ的后 验概 率 为
（ ６ ）
Ｉ Ｙ） ＝
。 ｃ Ｐ （ Ｙ ｌ
）
。 ｃ ｅ ｘ ｐ ｛ 一 丁 ／ ￣ ｎ Ｔ ｎ ） ｅ ｘ ｐ ｛ 一 — （ Ｙ — — － — Ｘ — ｒ — ） Ｔ （ — Ｙ — 一 － Ｘ ｒ ） ） ．
报
第３ ２ 卷
其中 ｎ （ ｔ １ 是地震波记录采集时所遇到的噪声 ，符 号 木 表示两个 向量之间的褶积 ．进 一步
地，根据式 （ １ ） 有
Ｙ ＝ Ｘ ｒ＋ Ｉ １ ，
（ ２ ）
Ｙ ＝ Ｒ ＋ ｎ，
（ ３ ）
其中 Ｙ ＝ （ １ ， ２ ， … ， Ⅳ ） Ｔ为长度为 Ⅳ 的观测地震波记录，ｒ＝ （ ｒ １ ， ｒ ２ ， … ， ｒ Ｍ） Ｔ为反射 系数序列， Ｘ∈ Ⅳ × Ｍ为地震子波褶积矩 阵 （ 一般为托普利兹矩阵） ， ｘ＝ （ １ ， ２ ， … ， ｘ Ｍ）
ｐ （ Ｙ Ｉ ｒ ） ＝ １ ｅ ｘ ｐ ｛ 一
（ １ ０ ）
另一 方 面 ，大 量 的文 献研 究 表 明 ，反射 系 数序 列 ｒ并不 是来 自经典 的高斯 分 布 ［ ７ ＿ 】 ．
特 别地 ，Ｗ ａ ｌ ｄ ｅ ｎ等 ［ １ ｌ 】 考 察 了来 自不 同地 质条 件下 的多个 检 测 的反射 系 数 的振幅 ，经验 性 地 证 明 了反射 系 数 的分布 一般 是 对称 的， 中心狭 窄 ，尾 部衰 减较 高 斯分 布慢 ．同 时 ，应 该 注 意 到 ，一 般 反射 系数 的非 零 值 表 示 的 是波 阻抗 系 数 变化 的界 面 ．所 以， 反射 系 数 可 以
文章编号： １ ０ ０ ５ — ３ ０ ８ ５ （ ２ ０ １ ５ ） ０ ３ ． ０ ３ ２ ８ — ０ ９
稳健 稀疏反褶 积方 法研 究术
高 倩 ， 刘建朝 ， 常象 宇
（ １ 一 长 安 大 学 地 球 科 学与 资源 学 院 数字 油 田研 究 所 ， 西 安 ７ １ ０ ０ ４ ５ ；
看 成 是稀 疏 信 号 ，即在 波 阻抗 系 数变 化 的界 面上 有脉 冲 值 ，而 其他 信 号位 置 为 ０的稀疏 信
号 ．基于 上述 论述 ，我 们使 用如 下分 布 作为 反射 系数 的先验 分布 ，即
） ＝ １ 唧｛ ＿ ） ７ｌ ｌ ｒ Ｉ ｌ ＝ Ｍ Ｉ ｊ
） ： １ 唧｛ － ） ）ｌ ｌ ｎ Ｉ Ｉ １ ＝ Ｎ
其中 为归一化常数 （ 为 了使得式 （ ９ ） 对应 于分布函数） ．该分布 函数的特点是具有重尾
性 ， 即在 尾 部 的衰减 速 率 明显 慢于 高斯 分 布 ，这样 的分布 是很 多异 常 点 的经 典 的统计 分布 描述 ．故本 文将提 出的稳健 稀 疏反 褶积 方法 也使 用上 述分 布作 为 噪音 的先验 分布 ．因此有
并且 模 型 是有 效 的 ．
关键词 ： 反褶积 ：反射系数 ；稳健估计 ；稀疏性
分类号：ＡＭ Ｓ （ ２ ０ ０ ０ １ ４ ２ — ０ ４
中图分类号： Ｐ ６ ３ １ ． ４
文献标识码 ： Ａ
１ 引言
反褶 积 方法 是 地震 勘探 数 字信 号 处理 领域 中一项 重 要 的技术 ，是地 震勘 探 的 基础 ．它
３３ ２
工
程
数
学
学
报
第３ ２ 卷
３ ． ２ 交 替迭代 算 法
由于在稳 健稀疏反褶积模型 （ １ ３ ） 中反射系数 ｒ 与褶积矩阵 ｘ都是未知的，这给直接
求 解 该 问题 带来 了极 大 的 困难 ．故在 本节 ，我们 将 提 出一种 交 替迭 代算 法 求解 该模 型 ．其
为 地 震 子 波 序 列 ，Ｒ ∈ Ⅳ× Ｍ为 反 射 稀 疏 褶 积 矩 阵 （ 一 般 为 托 普 利 兹 矩 阵） ，ｎ ＝
（ ｎ １ ， 礼 ２ ， … ， ｎ Ⅳ ） Ｔ 表 示观 测噪 声 ．对 于褶积 过程 ，公式 （ ２ ） 和（ ３ ） 是 等价 的 ．在 求解 反
射稀疏时一般使用公式 （ ２ ） ，而在求解地震子波时使用公式 （ ３ ） ．
众 所周 知 ，反 褶积 过程 是 一个 盲 过程 ．所 以在没 有任 何 先验 知 识 的情 况之 下 ，是 不能 够重 构 反射 系 数与 人 工地 震 子波 的 ． 因此 ，Ｃａ ｎ ａ ｄ ａ ｓ ［ ６ 】 建议 使 用 如下 先验 ．首先 对 于 反射 系数 ｒ ，假 设
ｐ （ ｒ ） 一 丽 １ ｅ ｘ ｐ 卜竿 ） ，
的基本 原 理是 ｆ 参见图 １ ） ：使用 仪 器 在地 表 面 的激 发 点 （ 图１ 中 ① 处） 人 工产 生 地震 子波 ．
如 果这 些 人 工 地 震 子 波 在 向地 下 传 播 的过 程 中遇 到 有 波 阻抗 变 化 的界 面 ，则 它 们 就会 产 生 反射 ．其 中一 部 分 能 量作 为 反射 地 震 波 向上 反 射 回 地 表 ，被 地 表 事 先 布 好 的检波 器 接 收ｆ 图 １中 ② 处） ．伴 随着 人工 地震 子 波不 断 的 向下传 播 与反射 ，地表 的检 波器 会 不 间断地 接 收 到地 震 子波 的 反射 信 号 ，并按 照 时 间排 列得 到 一系 列 时 间延迟 的 地震 波记 录 ．一般情
３ ． １ 稳健 稀 疏反 褶积 方法 经 过 多年研 究 ，众 多学 者 发现 在 一般 情 况下 ，地 震 数据 采集 过 程 中往往 会 遇 到异 常 点
的干 扰 ， 即地震 波记 录 的噪 声 不是 常 规 的高斯 白噪 声 ．虽然 ，我们 知道 这 种干 扰 严重 影 响
了利 用 反褶 积方 法 对真 实地 震 信号 重 构效 果 ，但 仍没 有较 好 的方 法 去剔 除 这些 非 高斯 白噪
２ ．西安交通大学数 学与统计 学院信息与系统科学研究所，西安 ７ １ ０ ０ ４ ９ １
摘
要： 地 震勘探数字信 号处理领域 中一个重要 而基本 的方法是地震信 号的反褶积 方法 ．在地 震数据的采集过 程中往往会 遇到异常 点的干扰 ，这种 干扰严重 影响 了利用 反褶积方法 对真 实反射系数 与地震子波 的重构效果 ．本文在 Ｃ ａ ｎ ａ ｄ ａ ｓ 等人提 出的针对 高斯噪声 的 贝叶斯反褶积数 学框架 的基础 之上 ，提 出一种 能够 克服异常 点干扰 的稳 健稀 疏反褶积 方法 ．新方法针对 具有重尾 分布的异常 点噪声与稀疏 的反射系 数建模 ，并使 用交替迭 代 与线性规划 的算法求解 ．最后 ，通过 实验证 明该方法在 克服 异常点 噪声的基础上 ， 能实现对地震子 波与反射 系数的 同步估 计 ，所得 到的估计有效 地消除 了重尾 分布异 常 点噪声的影 响，提高 了地 震信 号反褶积 处理 的精度 ．这也 能证 明所提算 法是收敛 的，
式 ，我们 可 知反射 系 数 ｒ的后验概 率 为
（ １ １ ）
其 中 为归一化 常数 （ 为了使得式子 （ １ １ ） 为分布 函数） ， ＞０ 为常数 ．根据 贝叶斯公
Ｉ ｖ） ＝
。 ｃ
。 ｃ （ ｖ ｌ
ｅ ｘｐ
）
｛ 一 ） ｅ ｘ ｐ ｛ 一
） ＿
第３ 期
高 倩 ， 等 ：稳 健 稀疏 反 褶积 方 法研 究
３ ３ １
声 的异 常 点 ．在统 计 学上 ，一 般对 这 种 噪声 的处 理方 法 是使 用稳 健统 计 量对 于 具有 重尾 分
布 的噪 音进 行建 模 ［ １ ２ ， １ ３ ］ ，例 如可 使用 如下 拉普 拉斯 分布 对 噪声 进行 建模 ，即
（ １ ２ ）
所以，根据极大化后验概率准则可以去估计反射系数．于是，对式 （ １ ２ ） 两边取对数，并最
大化 ｒ 得 到如 下优 化 问题
ｍ
ｒ ． Ｘ
ｉ ｎｌ Ｉ Ｖ —ｘ， Ｉ ＋￣ Ｉ ｌ ｌ ｒ ｌ Ｉ ．
、，
。
一
。 一
（ １ ３ ）
优化模型 （ １ ３ ） 就是本文提出的稳健稀疏反褶积方法．值得注意的是，如果已知褶积矩 阵ｘ，则该 问题等价于统计上的 Ｌ Ａ Ｄ — Ｌ ａ ｓ ｓ ｏ ￣ ［ １ ４ ， １ ５ 】 ．这个模型 由于使用了 ｌ １ ． Ｉ Ｉ １ 惩罚函 数 ，使得这个优化 问题 的解与经典的 Ｌ ａ ｓ ｓ ｏ 方法［ １ ６ ］ 一样具有稀疏解 （ 由参数 入 控制） ．为了 求解模型 （ １ ３ ） ，我们将在下一小节给 出一种交替迭代算法．