结构张量引导的多阶可转向马尔可夫随机场先验波阻抗反演方法

A Structure-Tensor-Guided Multi-Order Steerable Markov Random Field Prior Method for Acoustic Impedance Inversion

摘要: 地震反演是油气勘探与开发中的关键技术之一，而波阻抗反演是利用地震数据、低频地质信息(如低频背景模型、趋势模型)反演出地下介质的波阻抗，从而为地层结构分析、岩性识别与界面追踪提供重要依据。然而，波阻抗反演往往是一个病态问题，易受到噪声、带限效应与初始模型不确定性的影响；同时传统的单道反演无法充分利用采样道之间的构造信息，容易产生横向连续性不足等现象，尤其在地层结构复杂的区域更加明显。因此，为了解决上述问题，我们提出了一种结构张量引导的多阶可转向马尔可夫随机场先验波阻抗反演方法。该方法在贝叶斯框架下将先验建模为各向异性马尔可夫随机场，通过提取采样点的局部主方向，并引入考虑不同距离的多阶邻域约束，构造出沿层与垂直于层的差分约束，以增强对层理连续性与构造边界的刻画能力，同时提高反演稳定性。实验结果表明：在相同设置下，本方法相较于其他基线方法在RMSE、NRMSE、相关系数等指标上取得更优结果，并且对噪声具有更强的鲁棒性。

Abstract: Seismic inversion is one of the key technologies in oil and gas exploration and development. Wave impedance inversion is a method that uses seismic data and low-frequency geological information (such as low-frequency background models and trend models) to invert the wave impedance of underground media, thereby providing important basis for stratigraphic structure analysis, lithology identification, and interface tracking. However, wave impedance inversion is often a pathological problem, which is easily affected by noise, band-limited effects, and uncertainty of the initial model; at the same time, traditional single-channel inversion cannot fully utilize the structural information between sampling channels, and is prone to phenomena such as insufficient lateral continuity, especially in areas with complex stratigraphic structures. Therefore, to solve these problems, we propose a structure tensor-guided multi-order steerable Markov random field prior wave impedance inversion method. This method models the prior in the Bayesian framework as an anisotropic Markov random field, by extracting the local principal direction of the sampling points and introducing multi-order neighborhood constraints considering different distances, constructing differential constraints along the layer and perpendicular to the layer, to enhance the ability to depict layer continuity and structural boundaries, and improve the inversion stability. Experimental results show that, under the same settings, this method achieves better results in terms of RMSE, NRMSE, correlation coefficient, etc., compared to other baseline methods, and has stronger robustness to noise.

1. 引言

油气与矿产资源勘查长期以来是能源与资源邻域的重要研究方向。在勘探工作中，地震反演是对地震资料进行解释的关键环节[1] [2]。根据输入数据类型的不同，反演问题通常可分为叠前反演与叠后反演两类：叠前反演将随偏移距(或入射角)变化的地震道集作为输入，在波动方程或近似正演模型约束下反演地下介质的弹性参数(如$V_p,V_s,\rho$ )，而叠后反演是以叠加后的地震记录为输入。无论是叠前还是叠后反演，都可以看作是一个带限、含噪条件下的病态反问题，具体的，有限的带宽造成地震数据低频信息缺失，并且噪声与子波的不确定性进一步加剧了反演问题的不适定性；此外，在实际收集与处理过程中，地震资料往往存在缺失道、坏点等现象，从而导致反演结果的横向连续性变差，并且偏离真实地下结构。因此，引入恰当的先验信息或正则化约束成为了缓解这一病态问题常用且有效的途径。

而相较于叠前反演，叠后地震资料在弱反射与正入射近似的情况下，可以视作子波与反射系数序列的卷积，而反射系数与波阻抗的变化密切相关，于是利用叠后地震数据进行波阻抗反演，能够提供更接近地质属性的表征，已经广泛应用于油气储层预测、油气藏储量评估、钻井位置与轨迹设计等领域。

为了获得更准确的地质信息，已有的反演方法可以大致分为三种：确定性模型驱动反演，通过数据拟合项与正则化实现稳定重建；基于统计框架的反演，通过建立似然与先验从而构造后验分布，量化反演过程的不确定性；基于数据驱动的反演，利用训练数据学习观测数据与地下参数之间的映射关系。例如，Wang等人[3]注意到反射系数序列的稀疏特征，在目标函数中引入$L_1$ 正则化并结合量子退火进行全局优化，以提升含噪条件下的稳定性与抗噪能力。Dai等人[4]，构造多道阻抗模型的$L_{2,0}$ 联合约束，并采用分裂Bregman求解，在保持纵向稀疏的同时增强了横向连续性。此外，Bordignon等人[5]将地质统计反演与贝叶斯线性化地震反演相结合，实现了对含烃区带的准确预测。

尽管这些研究取得了卓越的成果，但是波阻抗反演仍然面临三个巨大挑战。1) 部分方法未在与正演模型及噪声统计假设一致的统一框架下构造目标函数，使得正则项及其超参数缺乏明确的物理或统计含义，往往依赖经验设定，从而导致反演结果的可解释性较弱。2) 许多空间约束多基于各向同性或固定方向施加惩罚，易在断层与倾角较大的层理附近产生跨层过平滑或构造不一致的现象，无法体现出沿层连续、跨层突变的结构特征。3) 现有方法通常只考虑一阶或单尺度邻域的变化，无法兼顾多尺度结构的细节刻画与界面保持，并且在地层构造复杂的情况下易产生条纹、阶梯等伪影，导致横向连续性退化。

因此，本文提出一种结构张量引导的多阶可转向马尔可夫随机场先验波阻抗反演方法(简称ST-MOSMRF)。该方法在贝叶斯框架下利用马尔可夫随机场构造结构先验，使目标函数中各项及其权重具有明确的物理含义，缓解仅凭经验“拼接”正则项带来的可解释性与可复现性不足。其次，ST-MOSMRF利用结构张量在每个采样点估计局部法向(垂直层方向)与切向(沿层方向)，将约束方向旋转投影到这两个方向上，从而实现了基于每个点局部方向的各向异性约束，并在机制上减少固定方向惩罚在断层及高倾角层理附近造成的跨层过平滑与结构不一致。进一步地，ST-MOSMRF在方向投影的基础上联合引入多阶邻域差分，具体的，使用一阶差分捕捉局部的梯度变化，二阶差分描述曲率与拐点的特征，三阶差分用于增强结构细节与边界的刻画。最后采用Geman-McClure (GM)饱和势函数对方向差分进行惩罚。

2. 预备知识与理论基础

本章首先给出了波阻抗反演的相关基础知识，介绍了正演模型与基于贝叶斯框架的反演模型，最后引入结构张量并推导局部方向的解析表达，为下一章构造结构张量引导的各向异性马尔可夫先验提供理论基础。

2.1. 叠后波阻抗反演

叠后单道的地震数据被定义为反射系数与带限地震子波的卷积结果，对应的卷积模型可以表示为：

$s\left(t\right)=w\left(t\right)\ast r\left(t\right)+n\left(t\right)$ (2.1)

其中，$s\left(t\right)$ 为地震记录，$w\left(t\right)$ 为地震子波，$r\left(t\right)$ 为反射系数，$n\left(t\right)$ 为噪声，$\ast$ 为卷积操作。

在多道反演中，我们同时考虑多条地震道，将式(2.1)中的变量扩展为矩阵形式，可得多道正演模型：

$d=WR+N$ (2.2)

其中，$d={\left[d_1^{\text{T}},\cdots ,d_M^{\text{T}}\right]}^{\text{T}}$ ，，分别表示地震数据、反射系数与噪声的矩阵形式，$d_j,r_j,n_j\in {ℝ}^N,j=1,\cdots ,M$ ，分别为第$j$ 道地震记录、反射系数与噪声的矩阵形式， $W=\text{diag}\left[w_1,\cdots ,w_M\right]$ 是由多个Toeplitz矩阵组成的块对角矩阵。

在正入射条件下，界面的反射系数由波阻抗决定，设上下两侧阻抗分别为$Z_1,Z_2$ ，某一层的反射系数定义为：

$r=\frac{Z_2-Z_1}{Z_2+Z_1}\approx \frac{1}{2}\left(\ln Z_2-\ln Z_1\right)$ (2.3)

然而，式(2.3)中的反射系数与波阻抗之间是非线性关系，这种非线性关系会增加反演过程中的计算开销，并且导致反演结果不稳定。当反射系数较小时，通常采用如下式子作为线性近似：

(2.4)

其中，$r_j$ 为单道地震道的反射系数序列，为单道地层的对数阻抗，$S_j$ 为差分算子定义为：

$S_j=\frac{1}{2}\left[\begin{array}{ccccc}-1& 1& 0& \cdots & 0\\ 0& -1& 1& \ddots & ⋮\\ ⋮& \ddots & \ddots & \ddots & 0\\ 0& \cdots & 0& -1& 1\end{array}\right]$

综上，结合式(2.2)与式(2.4)，可以得到多道波阻抗的正演模型：

$d=WSm+N=Gm+N$

其中，$G=WS$ 为正演算子，$S$ 是由$S_j$ 构成的块差分矩阵，为第$j$ 道地震道的对数波阻抗序列。

2.2. 贝叶斯框架

在贝叶斯框架中[6]，将反演视为对地下模型参数的后验概率推断问题：通过最大化后验分布概率从而获得最有可能的估计结果，而后验概率分布可以由似然函数与$m$ 的先验概率分布的乘积表示，形如：

$p\left(m|d\right)=\frac{p\left(d|m\right)p\left(m\right)}{p\left(d\right)}\propto p\left(d|m\right)p\left(m\right)$ (2.5)

其中，$p\left(d|m\right)$ 为似然函数，$p\left(m\right)$ 为先验概率分布。若噪声服从方差为${\sigma }_i^2$ 的高斯分布，似然函数可以表示为每一道独立似然的乘积，于是有：

$p\left(d|m\right)\propto \exp \left(-\frac{1}{2}{\left(Gm-d\right)}^{\text{T}}{C}_d^{-1}\left(Gm-d\right)\right)$ (2.6)

其中$C_d=\text{diag}\left[{\sigma }_1^{2}I,{\sigma }_2^{2}I,\cdots ,{\sigma }_j^{2}I\right]$ 。

因此，从而MAP估计等价于最小化负对数后验得到的目标函数，结合(2.6)式可以得到最终的反演目标函数：

$\begin{array}{c}\widehat{m}=\arg \min f\left(m\right)\\ =\arg \min \frac{1}{2}{\Vert L\left(Gm-d\right)\Vert }_2^2+\left[-\ln \left(p\left(m\right)\right)\right]\end{array}$ (2.7)

其中，$L=\text{diag}\left[{\sigma }_1^{-1}I,\cdots ,{\sigma }_j^{-1}I\right]$ 。

3. 结构张量引导的多阶可转向各向异性马尔可夫先验

在本章节中，我们引入我们的反演方法ST-MOSMRF。如图1所示，ST-MOSMRF方法的目标函数分别由三个部分构成：数据拟合项、低频先验项和结构约束项。数据拟合项由标准贝叶斯框架下的高斯似然得到；低频先验项通过引入低频背景模型$m_0$ (通常由测井数据或低通滤波趋势模型构建)，缓解了地震数据低频缺失、病态性问题并提高反演稳定性。我们首先由结构张量估计每个采样点的局部主方向与法向方向，并由其特征向量构造旋转矩阵$Q$ 。随后将一阶、二阶与三阶差分算子分别旋转到沿层和垂直层方向，构成多阶各向异性结构约束项。

3.1. 马尔可夫随机场

将地层按照采样点离散化后，得到一个二维网格$\Omega$ 上的场，其中相邻网格点之间具有明显的空间相关性。因此，我们将该网格视为一个无向图$\mathcal{G}$ ，其中各网格点作为图上的节点、相邻网格点之间通过无向边连接。由于这个无向图是满足Markov性的，因此可以得到一个关于$\mathcal{G}$ 的马尔可夫随机场。根据Hammersley-Clifford定理，我们将式(2.5)中模型参数$m$ 的联合先验概率分布修改为如下形式[7]：

Figure 1. Flowchart of the ST-MOSMRF method

图1. ST-MOSMRF方法的流程图

$p\left(m\right)=\frac{1}{Z}\exp \left[-{\displaystyle \sum _{C_k}\varphi }\left(D_k\left(m\right)\right)\right]$ (3.1)

其中，$Z$ 为归一化常数，$C_k$ 为第$k$ 阶邻域对应的团集合，$D_k(\cdot )$ 为与无向图中团$C_k$ 相关的梯度函数，$k$ 为MRF邻域的阶数，$\varphi$ 为势函数。

在本文中我们取$k=1,2,3$ ，其中$k=1$ 时，$C_1$ 代表一阶邻域(曼哈顿距离为1的四邻域)，并以此类推。其次选取Geman-McClure函数作为本文使用的势函数(${\varphi }_{GM}$ )，该函数在差分较小时近似二次惩罚，在差分较大时区域饱和，兼顾了地层内的噪声压制与跨层时界面边缘的保持，从而实现多尺度结构刻画的平衡。

3.2. 多阶旋转差分算子

对于地层中的任意采样点$a=\left(x,z\right)\in \Omega$ ，根据叠后地震数据$d$ ，该点的结构张量可以定义为：

$J=K_{\rho }\ast \left[\nabla d\left(a\right){\left(\nabla d\left(a\right)\right)}^{\text{T}}\right]=\left[\begin{array}{cc}{J}_{xx}& J_{xz}\\ J_{xz}& J_{zz}\end{array}\right]$

其中，$\nabla d\left(a\right)={\left[d_x\left(a\right),d_z\left(a\right)\right]}^{\text{T}}$ 为点$a$ 处的梯度向量，$K_{\rho }$ 为高斯核函数，$\rho$ 是高斯核的标准差，使用高斯平滑提高了对噪声的鲁棒性。对半正定矩阵$J$ 进行特征值分解可以得到：

$J=\left[\begin{array}{cc}{e}_1& e_2\end{array}\right]\left[\begin{array}{cc}{\lambda }_1& 0\\ 0& {\lambda }_2\end{array}\right]{\left[\begin{array}{cc}{e}_1& e_2\end{array}\right]}^{\text{T}}$

其中，${\lambda }_1$ 和${\lambda }_2$ 为特征值，$e_1={\left[e_{11},e_{12}\right]}^{\text{T}}$ 与$e_2={\left[e_{21},e_{22}\right]}^{\text{T}}$ 为对应的单位正交特征向量。令${\lambda }_1\ge {\lambda }_2\ge 0$ ，则特征向量$e_1$ 代表该点主方向(垂直于地层的方向)，而$e_2$ 代表次方向(平行于地层的方向)，令$\theta =\text{atan2}\left(e_{22},e_{21}\right)$ 我们可以得到旋转矩阵$Q$ ：

$Q\left(\theta \right)=\left[\begin{array}{cc}\cos \theta & \sin \theta \\ -\sin \theta & \cos \theta \end{array}\right]$ (3.2)

因此，我们利用旋转矩阵$Q$ 计算出沿层与垂直于层的一阶方向差分：

$\left[\begin{array}{c}{D}_t^{(1)}\\ D_n^{(1)}\end{array}\right]=Q\left(\theta \right)\left[\begin{array}{c}{D}_x\\ D_z\end{array}\right]=\left[\begin{array}{cc}\cos \theta & \sin \theta \\ -\sin \theta & \cos \theta \end{array}\right]\left[\begin{array}{c}{D}_x\\ D_z\end{array}\right]$

将上式展开，可以得到一阶邻域上的方向差分：

$\begin{array}{l}{D}_t^{(1)}=\cos \theta \cdot D_x+\sin \theta \cdot D_z\\ D_n^{(1)}=-\sin \theta \cdot D_x+\cos \theta \cdot D_z\end{array}$ (3.3)

其中，$D_x$ 和$D_z$ 为笛卡尔坐标系下水平与垂直方向的前向差分算子[7]，$D_n^{(1)}$ 与$D_t^{(1)}$ 分别表示在一阶邻域中垂直于层和平行于层的方向差分算子。我们利用海森矩阵的二次型旋转与二项式立方展开分别得到二阶、三阶方向差分算子，具体形式如下：

$\begin{array}{l}{D}_t^{(2)}={\cos }^2\theta \cdot D_{xx}+2\cos \theta \sin \theta \cdot D_{xz}+{\sin }^2\theta \cdot D_{zz}\\ D_n^{(2)}={\sin }^2\theta \cdot D_{xx}-2\cos \theta \sin \theta \cdot D_{xz}+{\cos }^2\theta \cdot D_{zz}\\ D_t^{(3)}={\cos }^3\theta \cdot D_{xxx}+3{\cos }^2\theta \sin \theta \cdot D_{xxz}+3\cos \theta {\sin }^2\theta \cdot D_{xzz}+{\sin }^3\theta \cdot D_{zzz}\\ D_n^{(3)}=-{\sin }^3\theta \cdot D_{xxx}+3{\sin }^2\theta \cos \theta \cdot D_{xxz}-3\sin \theta {\cos }^2\theta \cdot D_{xzz}+{\cos }^3\theta \cdot D_{zzz}\end{array}$ (3.4)

上式中，$D_{xz},D_{xxz},D_{xzz}$ 为笛卡尔坐标系下混合差分[7]，$D_n^{(k)},D_t^{(k)},k=1,2,3$ 分别代表不同邻域中垂直于地层的方向差分和平行于地层的方向差分算子。

3.3. 目标函数

在二维离散网格$\Omega$ 上，对某一网格点记其对数阻抗为$m\left(i,j\right)$ 。如图2所示，我们给定一个位于$\left(i,j\right)$ 的指定点，分别用曼哈顿距离为1、2、3的四邻域作为一阶邻域、二阶邻域和三阶邻域，网格中的数字代表其与中心点的距离。对于一、二、三阶邻域我们同时考虑上下与左右四个方向的差分，并根据(3.3)式与(3.4)式得到四个方向差分算子。例如，对于二阶马尔可夫随机场邻域，全局能量函数$J_{str,2}$ 可以写为：

Figure 2. (i) First-order neighborhood; (ii) Second-order neighborhood; (iii) Third-order neighborhood

图2. (i) 一阶邻域；(ii) 二阶邻域；(iii) 三阶邻域

$\begin{array}{c}{J}_{str,2}\left(m\right)={\lambda }_2{\displaystyle \sum _{c_2}{\varphi }_{GM}}\left(D_2\left(m\right)\right)\\ ={\lambda }_2{\displaystyle \sum _{c_2}\left[{\varphi }_{GM}\left(D_{n,u}^{(2)}\left(m\right)\right)+{\varphi }_{GM}\left(D_{n,d}^{(2)}\left(m\right)\right)+{\varphi }_{GM}\left(D_{t,l}^{(2)}\left(m\right)\right)+{\varphi }_{GM}\left(D_{t,r}^{(2)}\left(m\right)\right)\right]}\end{array}$

其中，$D_{n,u}^{(2)},D_{n,d}^{(2)}$ 分别代表邻域中上下两侧的方向差分算子，同理，$D_{n,l}^{(2)},D_{n,r}^{(2)}$ 分别代表邻域中左右两侧的方向差分算子，${\lambda }_2$ 为二阶邻域全局能量函数的权重参数。Geman-McClure势函数${\varphi }_{GM}$ 的形式如下：

${\varphi }_{GM}\left(t\right)=\frac{t^2}{1+t^2}$

类似的，可以得到一阶邻域与三阶邻域的全局能量函数$J_{str,1}\left(m\right)$ 和$J_{str,3}\left(m\right)$ 。我们同时考虑多阶邻域，从而将式(3.1)修改为如下形式：

$p\left(m\right)\propto \exp \left[-\left(J_{str,1}\left(m\right)+J_{str,2}\left(m\right)+J_{str,3}\left(m\right)\right)\right]$ (3.5)

综上，结合式(2.7)与式(3.5)，并引入低频约束项从而得到反演目标函数：

$f\left(m\right)=\frac{1}{2}{\Vert L\left(Gm-d\right)\Vert }_2^2+\alpha {\Vert m-m_0\Vert }_2^2+\beta \left(J_{str,1}\left(m\right)+J_{str,2}\left(m\right)+J_{str,3}\left(m\right)\right)$

其中，$\alpha$ 为低频约束的权重系数，用于控制模型对低频信息的依赖程度，$\beta$ 为结构先验项的整体权重系数，对于结构约束项，${\lambda }_k$ 不同邻域的相对权重，而$\beta$ 控制整个先验的强度。

4. 实验结果

4.1. 实验设置

为验证提出方法的有效性与可复现性，本文采用Marmousi2模型作为实验数据集。Marmousi2模型是地震反演邻域中常用的基准数据集，其原始模型包含2721条地震道，每条道包含701个采样点，采样间隔为5米。

本文采用主频为30 Hz的雷克子波作为震源子波，该子波的表达式如下：

$w\left(t\right)=\left(\begin{array}{c}1-2{\left(\pi f_{0}t\right)}^2\end{array}\right)\exp \left(\begin{array}{c}-{\left(\pi f_{0}t\right)}^2\end{array}\right)$

其中，$f_0$ 为主频，子波长度为0.128 s。

考虑到波阻抗反演问题的病态性与低频缺失现象，在ST-MOSMRF方法中，引入低频背景先验以增强反演稳定性。具体的，本文使用阈值为5 Hz的Butterworth低通滤波构造低频背景模型$m_0$ 。目标函数由数据拟合项、低频先验约束以及结构先验约束项三部分构成，实验中低频先验项权重$\alpha$ 设置为$1\times {10}^{-4}$ ，结构约束权重$\beta$ 设置为$1\times {10}^{-3}$ ，结构先验中一、二、三阶邻域对应的全局能量函数权重分别为${\lambda }_1=1$ ，${\lambda }_2=1\times {10}^{-1}$ ，${\lambda }_3=6\times {10}^{-2}$ 。针对目标函数的最小化问题，本文使用Adam (Adaptive Moment Estimation)优化算法进行求解。其中，反演变量初始化为$m=m_0$ ，学习率设置为$\text{lr}=0.1$ ，动量参数为0.9，迭代次数$N_{\text{iter}}=500$ 。为了定量评估反演结果，本文选取如下四个指标：

• 模型RMSE：衡量反演阻抗与真实阻抗间的绝对误差水平，数值越小表示反演精度越高。

• 模型NRMSE：在模型RMSE的基础上进行归一化处理从而消除了量纲的影响，使不同数据集、不同工区或不同尺度实验之间的误差具有可比性。NRMSE越小表示反演结果与真实阻抗越接近。

• 相关系数Corr：衡量反演结果与真实模型的整体一致性或线性相关程度。数值越接近1代表一致性越好。

• 数据RMSE：计算预测数据与观测数据之间的均方根误差，用于衡量地震数据间的拟合质量。数值越小表示数据匹配程度越高。

本文选取了三种典型的反演方法(Tikhonov [8]、TV [9]和GSGTV [10])进行对比实验，在相同设置下对比了反演结果的四个评价指标。

4.2. 反演结果展示与分析

图3展示了无噪声基准条件下真实波阻抗模型、低频初始模型与ST-MOSMRF反演结果的剖面对比。可以看出，低频初始模型主要保留地层的大尺度趋势成分，层间阻抗对比与构造细节被明显平滑。相比之下，ST-MOSMRF在保持低频背景趋势的同时有效恢复了部分中高频变化，使主要界面轮廓、构造起伏形态及横向连续性均得到改善；尤其在弯曲构造区域，反演剖面侧向连贯性更稳定，较少出现横向破碎或条带化现象。上述结果表明，所构建的结构先验能够对地层连续性与构造一致性提供有效约束，从而提升对地质结构特征的刻画能力。图4给出了5条随机抽取地震道的波阻抗反演曲线对比。其中，黑色线代表真实波阻抗，红色线代表反演结果，蓝色虚线代表低频模型。如图所示，第252道与第970道反演效果较好，反演曲线在幅值水平与变化趋势上与参考阻抗一致性较高，界面定位较为准确；第1550道与第2150道在主趋势上与参考值吻合，但局部深度段存在一定的高频振荡；第2537道整体趋势与主要界面位置基本对应，但局部幅值略有偏差且细节起伏相对平滑。

Figure 3. Wave impedance profile

图3. 波阻抗剖面

Figure 4. Single-channel impedance curve

图4. 单道波阻抗曲线

表1为相同环境下对比了四种方法的反演结果，总体来看，ST-MOSMRF在精度上表现最优，取得了最低的模型RMSE (98.48)与最低的模型NRMSE (0.92%)，明显优于其他三种方法。具体而言，相比Tikhonov与TV，ST-MOSMRF的模型RMSE约降低50%，相比GSGTV约降低32.5%。同时，ST-MOSMRF取得最高相关系数(99.92%)，表明其在保持模型结构一致性与细节保真方面更具优势。就数据拟合而言，TV的数据RMSE最小(0.0023)，而ST-MOSMRF的数据RMSE仍保持较低水平，并优于Tikhonov与GSGTV，说明本文方法在显著提升模型反演质量的同时，仍能维持良好的数据一致性。

Table 1. Comparison of inversion results under noise-free conditions

表1. 在无噪声环境下的反演结果对比

	Tikhonov [8]	TV [9]	GSGTV [10]	ST-MOSMRF
模型RMSE	198.30	196.05	145.89	98.48
模型NRMSE	1.83%	1.81%	1.35%	0.92%
相关系数	99.72%	99.68%	99.87%	99.92%
数据RMSE	0.0058	0.0023	0.0069	0.0031

在图5中，我们分别向地震数据加入了5%、10%、15%和25%的高斯噪声来测试方法对于噪声的鲁棒性。总体而言，反演波阻抗剖面整体结构稳定，主要界面形态与横向连续性保持较好，高阻抗异常体边界清晰。噪声增至15%后，局部细节开始出现轻微模糊与幅值起伏，但主要层位的连续性仍然较好。到25%噪声时，剖面中随机纹理/条带状扰动明显增强，界面锐度下降、小尺度结构受影响更大，不过反演结果在大尺度上仍能恢复主要构造形态，宏观构造特征保持较好一致性。

Figure 5. Inverted wave impedance profiles at 5%, 10%, 15%, and 25% Gaussian noise

图5. 5%、10%、15%和25%的高斯噪声下的反演波阻抗剖面图

5. 结论

本文针对叠后波阻抗反演中病态性强、噪声敏感以及横向连续性不足等问题，在贝叶斯框架下提出一种结构张量引导的多阶可转向各向异性马尔可夫随机场先验反演方法(ST-MOSMRF)。该方法利用结构张量估计地层的局部主方向，将多阶差分约束旋转并投影到沿层/垂直层方向，同时引入Geman-McClure饱和势函数以鲁棒抑制噪声与异常梯度，从而在保持低频背景趋势的基础上增强对中高频结构的恢复能力。基于Marmousi2等基准模型的定性与定量实验结果表明，与Tikhonov、TV及GSGTV方法相比，ST-MOSMRF在界面刻画、构造起伏恢复与横向连续性保持方面表现更稳定；在RMSE、NRMSE与相关系数等指标上取得更高反演精度，并在5%~25%高斯噪声扰动条件下仍能有效保持宏观构造特征一致，体现出优良的抗噪稳健性。