0 引言

地震反演始于20世纪60年代，经过六十多年的发展，现有的反演理论已较为成熟，地震反演结果精度也取得了较大提升，能够较真实地反映地下介质的情况 ^[1]。由于地震反演问题充满不确定性和多解性，通常需要给反演目标泛函施加正则化约束项，根据实现方式可分为稀疏约束反演、地质统计学反演和人工智能反演等 ^[2-4]。渤海湾盆地渤中凹陷古近系广泛发育扇三角洲沉积，具备丰富的油气勘探前景 ^[5-6]，勘探实践表明，渤中凹陷石臼坨凸起构造带扇三角洲沉积发育多期叠置砂体，但砂体单层厚度小、横向展布变化快 ^[7]，同时受到该区钻井较少、目的层地震资料分辨率较低等问题制约，对古近系储层的刻画多以$ 90^{\circ} $相位法 ^[8]、RGB颜色融合法 ^[9]等常规方法为主，储层预测结果多解性强、精度较低 ^[10-11]。目前海上井网密度小，甚至大量无井区域，地震反演的初始模型构建困难，增强了地震反演结果的多解性 ^[12-13]。为解决这一问题，当下主流采用2种方法，一是通过地质统计学提高初始模型建模精度 ^[14]；二是通过研究地震信号时频谱与稀疏性等潜在信息，以降低反演结果对初始模型的依赖程度 ^[15-16]，前者主要适用于开发中后期井网密度大的地区，后者主要适用于油田勘探初期井网密度低的地区 ^[17]。当研究区缺少井信息时，如何充分利用地震信号各频带信息，并将其应用于地震反演中，一直是海上地震勘探的热点 ^[18]。地震信号谱分解技术是通过分析地震数据的时频域特征，挖掘地震数据中的隐藏信息 ^[19]，如被强反射屏蔽的弱反射等。刘文霞 ^[20]在处理深层地震资料使用分频技术，拓宽了地震信号的频带，并保证其信噪比和分辨率。胡咏等 ^[21]针对河流相薄储层时，提出利用分频方法处理地震数据，对识别砂体厚度小于四分之一地震波长的地质体取得了突破。于建国等 ^[22]根据叠后地震信号不同尺度特征，将叠后数据中有效频带分为大、中、小3个尺度，分别进行分频反演处理，取得了优于常规反演结果的分辨率。熊冉等 ^[23]在分频地震反演理论的基础上对测井资料中的先验信息进行研究，利用测井资料约束声波曲线重构过程，并通过分频技术对地震信号中的低频和高频成分进行合理补充，提高了薄层的识别精度。针对地震资料采集过程中因人工震源与接收器等原因的限制，使得地震资料中缺乏有效低频成分的问题，叶云飞等 ^[24]用基于广义S变换时频分析技术补偿了地震数据低频成分，提出了一种基于宽带信号的子波提取方法，获得的子波旁瓣振幅更小。代玲等 ^[25]通过遗传神经网络技术还原最接近真实地层的测井曲线，并以地震波形约束实现分频反演，在泥质砂岩储层预测中取得了较好效果。李丛等 ^[26]在分频反演中引入全频带信息，提出了多子波分频反演技术方法，并在页岩气薄储层识别中取得了较好效果。然而，现有分频约束反演方法大多未能充分考虑不同尺度地震信号的特性，在各尺度信号中采用统一的正则化约束条件，导致反演结果分辨率的提升有限。

以渤中凹陷石臼坨凸起构造带古近系东营组二段为例，基于多尺度地震信号，提出一种地震分频多级稀疏正则化反演方法，利用匹配追踪与Wigner-Ville分布（MP-WVD）时频方法将地震信号分解为三级尺度；基于稀疏理论和贝叶斯理论，对不同尺度地震信号进行分级约束，以期能够提高反演精度，增强薄砂岩储层识别的可靠性。

1 方法原理 1.1 MP与Wigner-Ville分布原理

利用地震分频属性来预测储层精度的关键是谱分解的时频分辨率，匹配追踪方法（MP）-Wigner-Ville（WVD）是一种无Heisenberg不确定原理限制的高时频分辨率谱分解技术，其中MP是关键的信号分解算法，WVD是频谱计算算法，MP将地震信号分解为一系列预定义的小波原子的线性叠加 ^[25]。

$ \boldsymbol{s}=\operatorname{Re}\left[\sum\limits_{n=1}^{N} c_{n} \boldsymbol{m}_{n}+\boldsymbol{r}_{N}\right] $

(1)

式中：$ \boldsymbol{s} $为（实值）地震信号；$ n $为迭代次数；$ N $为最大迭代次数；$ c_{n} $为第$ n $次迭代时的复系数；$ \boldsymbol{m}_{n} $为第$ n $次迭代时的最优小波原子；$ \boldsymbol{r}_{N} $为第$ N $次迭代时的残余误差，该误差仅为噪声；$ \operatorname{Re}[\cdot] $为取实部算子。

$ \boldsymbol{m}_{n} $包含了地震信号的时间和频率特征，可直接从小波原子字典中提取，通常以Morlet小波为原子。

$ \boldsymbol{m}_{n}=\exp \left[-4 \ln 2\left(\frac{f_{n}^{2}}{k_{n}^{2}}\right)\left(t-u_{n}\right)^{2}+\mathrm{i} 2 {\rm{\mathsf{π}}} f_{n}\left(t-u_{n}\right)\right] $

(2)

$ u_{n}=\arg \max\limits_{t}\left\{\operatorname{Env}\left[\boldsymbol{r}_{n-1}(t)\right]\right\} $

(3)

$ \left(k_{n}, f_{n}\right)=\arg \max\limits_{\substack{k \in D_{k} \\ f \in D_{f}}} \frac{\left|\boldsymbol{r}_{n-1}(t) \boldsymbol{m}_{n}\left(t, u_{n}, k, f\right)\right|}{\boldsymbol{m}_{n}\left(t, u_{n}, k, f\right)_{2}} $

(4)

式中：$ t $为时间，s；$ u_{n} $为中心时间，s；$ k $和$ k_{n} $分别为小波尺度、第$ n $个小波的尺度；$ f $与$ f_{n} $分别为信号频率和中心频率，$ \mathrm{Hz} ; \mathrm{i} $为虚数单位；$ \boldsymbol{r}_{n-1} $为第$ n-1 $次迭代时残差，有$ \boldsymbol{r}_{0}=s+\mathrm{iH}[s] $，其中$ \mathrm{H}[\cdot] $为Hilbert变换算子；$ \operatorname{Env}[\cdot] $为取包络算子；$ D_{k} $和$ D_{f} $分别为$ k_{n} $和$ f_{n} $的搜索区间，通常根据经验给定，其中$ D_{k}=[0.3, 6.0] $，

$ D_{f}=\left[f_{u}-R_{f}, f_{u}+R_{f}\right] $

(5)

式中：$ f_{u} $表示$ \boldsymbol{r}_{n-1} $在时间$ u_{n} $处的瞬时频率，$ \mathrm{Hz} ; R_{f} $为$ f_{n} $的搜索半径，Hz，一般设置为20 Hz。

式（1）中，$ c_{n} $表示地震信号的振幅和相位特征的复系数，计算式为

$ c_{n}=\left(\boldsymbol{m}_{n}^{\mathrm{T}} \boldsymbol{m}_{n}+\varepsilon\right)^{-1} \boldsymbol{m}_{n}^{\mathrm{T}} \boldsymbol{r}_{n-1} $

(6)

式中：$ \varepsilon $为阻尼常数。结合WVD，地震信号$ \boldsymbol{s} $的时频谱可表示为

$ \begin{align*} & \boldsymbol{A}(t, f)= \\ & \sum\limits_{n=1}^{N} \frac{4\left\{\left[\operatorname{Re}\left(c_{n}\right)\right]^{2}+\left[\operatorname{Im}\left(c_{n}\right)\right]^{2}\right\} \mathrm{e}^{\left.-\frac{{\rm{\mathsf{π}}}^{2} k_{n}^{2}\left(f-f_{n}\right)^{2}}{2 \ln f_{n}^{2}}-\frac{8 \ln 2 f_{n}^{2}\left(t-u_{n}\right)}{k_{n}^{2}}\right]}}{1+\mathrm{e}^{\left(\frac{-{\rm{\mathsf{π}}}^{2} k_{n}^{2}}{2 \ln 2}\right)} \cos \left\{\arctan \left[\frac{\operatorname{Im}\left(c_{n}\right)}{\operatorname{Re}\left(c_{n}\right)}\right]\right\}} \end{align*} $

(7)

式中：$ \operatorname{Im}\left[c_{n}\right] $为取虚部算子。

1.2 地震分频多级稀疏正则化反演目标方程

基于线性化波阻抗正演理论，地震波单道合成记录$ S_{\mathrm{syn}} $，可通过地震子波与反射系数卷积计算 ^[27]：

$ s_{\mathrm{syn}}(t)=w(t) * r_{\mathrm{s}}(t)+l(t) $

(8)

式中：$ w(t), l(t), r_{\mathrm{s}}(t) $分别为震源子波、随机噪音和反射系数。

$ r_{\mathrm{s}}\left(t_{j}\right)=\frac{z\left(t_{j+1}\right)-z\left(t_{j}\right)}{z\left(t_{j+1}\right)+z\left(t_{j}\right)} \approx \frac{1}{2}\left\{\ln \left[z\left(t_{j+1}\right)\right]-\ln \left[z\left(t_{j}\right)\right]\right\} $

(9)

式中：$ j $表示采样点；$ t $为信号时间，$ \mathrm{s} ; z\left(t_{j+1}\right) $为声波阻抗，$ \mathrm{m} / \mathrm{s} \cdot \mathrm{g} / \mathrm{cm}^{3} $。

将式带入式，并写成矩阵形式，可得

$ \boldsymbol{S}_{\text {syn }}=\boldsymbol{W D Z+L} $

(10)

式中：$ \boldsymbol{W} $为子波矩阵；$ \boldsymbol{Z}=\left\{\begin{array}{c}\ln \left[z\left(t_{1}\right)\right], \ln \left[z\left(t_{2}\right)\right], \cdots, \\ \ln \left[z\left(t_{K}\right)\right]\end{array}\right\}^{\mathrm{T}} $，其中$ K $为时间域波阻抗采样总数；$ \boldsymbol{L} $为随机噪音矩阵；$ \boldsymbol{D} $为一阶差分算子，表示为

$ \boldsymbol{D}=\frac{1}{2}\left[\begin{array}{ccccc} -1 & 1 & & & \\ & -1 & 1 & & \\ & & \ddots & \ddots & \\ & & & -1 & 1 \end{array}\right] $

(11)

基于贝叶斯理论，地震波观测数据$ \boldsymbol{S}_{\text {obs }} $，模型数据的后验算分布$ \mathrm{P}\left(\boldsymbol{Z} \mid \boldsymbol{S}_{\text {obs }}\right) $可表示为 ^[28]

$ \mathrm{P}\left(\boldsymbol{Z} \mid S_{\text {obs }}\right)=\frac{\mathrm{P}\left(S_{\text {obs }} \mid \boldsymbol{Z}\right) \mathrm{P}\left(S_{\text {obs }}\right)}{\mathrm{P}(\boldsymbol{Z})} $

(12)

式中：$ \mathrm{P}(\boldsymbol{Z}) $表示为模型参数的先验概率；$ \mathrm{P}\left(\boldsymbol{S}_{\text {obs }} \mid \boldsymbol{Z}\right) $为似然函数；$ \mathrm{P}\left(\boldsymbol{S}_{\mathrm{obs}}\right) $为边缘概率密度。

通常，假设$ \boldsymbol{L} $服从高斯分布，式中（12）似然函数可表示为 ^[5]

$ \mathrm{P}\left(S_{\text {obs }} \mid \boldsymbol{Z}\right)=\frac{1}{(\sqrt{2 {\rm{\mathsf{π}}}})^{K+N_{w}-2} \boldsymbol{C}_{\mathrm{e}}} \mathrm{e}^{\left[-\frac{1}{2}\left(S_{\text {dath }}-W D Z\right)^{\mathrm{T}} \mathrm{C}_{\mathrm{e}}^{-1}\left(S_{\text {doc }}-W D Z\right)\right]} $

(13)

式中：$ N_{\mathrm{w}} $为地震子波长度；$ \boldsymbol{C}_{\mathrm{e}} $为噪声的协方差，$ \boldsymbol{C}_{\mathrm{e}}= $ $ \sigma_{\mathrm{e}}^{2} \boldsymbol{I} $，其中$ \sigma_{\mathrm{e}} $为噪声的标准差，$ \boldsymbol{I} $为单位矩阵。

通常，假设地层中的波阻抗参数符合高斯分布，模型参数的先验分布表示为

$ \mathrm{P}(\boldsymbol{Z})=\frac{1}{(\sqrt{2 {\rm{\mathsf{π}}}})^{k} \boldsymbol{C}_{\mathrm{Z}}} \mathrm{e}^{\left[-\frac{1}{2} \phi(\boldsymbol{Z})^{\mathrm{T}} \boldsymbol{C}_{\mathrm{Z}}^{-1} \phi(\boldsymbol{Z})\right]} $

(14)

式中：$ \phi(\boldsymbol{Z}) $为$ \boldsymbol{Z} $的函数；$ \boldsymbol{C}_{\mathrm{z}} $为背景约束项的协方差，$ \boldsymbol{C}_{\mathrm{Z}}=\sigma_{\mathrm{Z}}^{2} \boldsymbol{I} $，其中$ \sigma_{\mathrm{Z}} $为背景约束项的标准差。

$ \mathrm{P}\left(\boldsymbol{S}_{\mathrm{obs}}\right) $为常数，将式（14）和式（13）代入（12），波阻抗的后验概率分布可表示为

$ \mathrm{P}\left(\boldsymbol{Z} \mid \boldsymbol{S}_{\text {obs }}\right) \propto \mathrm{e}^{\left[-\frac{1}{2}\left(\boldsymbol{S}_{\text {olm }}-\boldsymbol{W} \boldsymbol{D} \boldsymbol{Z}\right)^{\mathrm{T}} \boldsymbol{C}_{\mathrm{e}}^{-1}\left(S_{\text {dats }}-\boldsymbol{W} \boldsymbol{D} \boldsymbol{Z}\right)-\frac{1}{2} \lambda \phi(\boldsymbol{Z})^{\mathrm{T}} \boldsymbol{C}_{\mathrm{Z}}^{-1} \phi(\boldsymbol{Z})\right]} $

(15)

式中：$ \lambda $为正则化参数；$ \propto $表示正比于。

对式（15）进一步进行整理，得到地震反演的目标函数：

$ \begin{gather*} J(\boldsymbol{Z})=\frac{1}{\sigma_{\mathrm{e}}^{2}}\left(\boldsymbol{S}_{\text {obs }}-\boldsymbol{W} \boldsymbol{D} \boldsymbol{Z}\right)^{\mathrm{T}}\left(\boldsymbol{S}_{\text {obs }}-\boldsymbol{W} \boldsymbol{D} \boldsymbol{Z}\right)+ \\ \lambda \frac{1}{\sigma_{\rm Z}^{2}} \phi(\boldsymbol{Z})^{\mathrm{T}} \phi(\boldsymbol{Z})\end{gather*} $

(16)

式中：传统反演方法中$ \phi(\boldsymbol{Z}) $为背景约束项，由模型参数与模型背景参数差组成$ \boldsymbol{Z}_{\mathrm{ref}} $，具体为

$ \phi(\boldsymbol{Z})=\left(\boldsymbol{Z}-\boldsymbol{Z}_{\mathrm{ref}}\right) $

(17)

基于式，构建地震分频多级稀疏正则化反演目标函数，具体如下：

$ \left\{\begin{array}{c} J_{\mathrm{b}}(\boldsymbol{Z})=\frac{1}{\sigma_{\text {eb }}^{2}}\left(\boldsymbol{S}_{\mathrm{b}}-\boldsymbol{W}_{\mathrm{b}} \boldsymbol{D} \boldsymbol{Z}_{\mathrm{b}}\right)^{\mathrm{T}}\left(\boldsymbol{S}_{\mathrm{b}}-\boldsymbol{W}_{\mathrm{b}} \boldsymbol{D} \boldsymbol{Z}_{\mathrm{b}}\right) \\ \quad+\lambda_{\mathrm{b}} \frac{1}{\sigma_{\text {zb }}^{2}}\left\|\boldsymbol{Z}_{\mathrm{b}}-\boldsymbol{Z}_{\text {ref } \_\mathrm{b}}\right\|_{2}^{2}+\beta_{\mathrm{b}}\left\|\boldsymbol{D} \boldsymbol{Z}_{\mathrm{b}}\right\|_{2}^{2} &(A)\\ J_{\mathrm{m}}(\boldsymbol{Z})=\frac{1}{\sigma_{\text {em }}^{2}}\left(\boldsymbol{S}_{\mathrm{m}}-\boldsymbol{W}_{\mathrm{m}} \boldsymbol{D} \boldsymbol{Z}_{\mathrm{m}}\right)^{\mathrm{T}}\left(\boldsymbol{S}_{\mathrm{m}}-\boldsymbol{W}_{\mathrm{m}} \boldsymbol{D} \boldsymbol{Z}_{\mathrm{m}}\right) \\ +\lambda_{\mathrm{m}} \frac{1}{\sigma_{\text {zm }}^{2}}\left\|\boldsymbol{Z}_{\mathrm{m}}-\boldsymbol{Z}_{\text {ref } \_\mathrm{m}}\right\|_{2}^{2}+\gamma_{\mathrm{m}}\left\|\boldsymbol{Z}_{\mathrm{m}}-\boldsymbol{Z}_{\mathrm{b}}\right\|_{2}^{2}+\beta_{\mathrm{m}}\left\|\boldsymbol{D} \boldsymbol{Z}_{\mathrm{m}}\right\|_{1} &(B)\\ J_{\mathrm{s}}(\boldsymbol{Z})=\frac{1}{\sigma_{\text {es }}^{2}}\left(\boldsymbol{S}_{\mathrm{s}}-\boldsymbol{W}_{\mathrm{s}} \boldsymbol{D} \boldsymbol{Z}_{\mathrm{s}}\right)^{\mathrm{T}}\left(\boldsymbol{S}_{\mathrm{s}}-\boldsymbol{W}_{\mathrm{s}} \boldsymbol{D} \boldsymbol{Z}_{\mathrm{s}}\right) \\ +\lambda_{\mathrm{s}} \frac{1}{\sigma_{\text {zs }}^{2}}\left\|\boldsymbol{Z}_{\mathrm{s}}-\boldsymbol{Z}_{\text {ref_} \mathrm{s}}\right\|_{2}^{2}+\gamma_{\mathrm{s}}\left\|\boldsymbol{Z}_{\mathrm{s}}-\boldsymbol{Z}_{\mathrm{m}}\right\|_{2}^{2}+\beta_{\mathrm{s}}\left\|\boldsymbol{D} \boldsymbol{Z}_{\mathrm{s}}\right\|_{0} &(C)\end{array}\right. $

(18)

式中：$ \boldsymbol{S}_{\mathrm{b}}, \boldsymbol{S}_{\mathrm{m}} $和$ \boldsymbol{S}_{\mathrm{s}} $分别为大、中和小尺度地震信号；$ \boldsymbol{W}_{\mathrm{b}}, \boldsymbol{W}_{\mathrm{m}} $和$ \boldsymbol{W}_{\mathrm{s}} $为分别为大、中和小尺度地震信号中提取统计子波；$ \boldsymbol{Z}_{\mathrm{b}}, \boldsymbol{Z}_{\mathrm{m}} $和$ \boldsymbol{Z}_{\mathrm{s}} $分别为大、中、小尺度的波阻抗模型参数；$ \boldsymbol{Z}_{\text {ref＿b }}, \boldsymbol{Z}_{\text {ref＿m }} $和$ \boldsymbol{Z}_{\text {ref＿s }} $分别为大、中和小尺度波阻抗参考模型；$ \sigma_{\mathrm{eb}}, \; \sigma_{\mathrm{em}} $和$ \sigma_{\mathrm{es}} $分大、中和小尺度误差的标准差；$ \sigma_{\mathrm{zb}}, \sigma_{\mathrm{zm}} $和$ \sigma_{\mathrm{zs}} $分大、中和小尺度背景项误差的标准差；$ \lambda_{\mathrm{b}}, \lambda_{\mathrm{m}} $和$ \lambda_{\mathrm{s}} $分别为大、中和小尺度的波阻抗低频背景约束项权重参数；$ \gamma_{\mathrm{m}} $和$ \gamma_{\mathrm{s}} $分别为中尺度和小尺度背景约束项的权重参数；$ \beta_{\mathrm{b}} $，$ \beta_{\mathrm{m}} $和$ \beta_{\mathrm{s}} $分别为$ \mathrm{L} 2, \mathrm{\; L} 1 $和L 0范数正则化约束项的权重参数。

L0范数凸表面积最小（图 1a），稀疏性最强，抗噪性也最强，属于典型非确定性多项式问题，计算复杂，耗时高。L1范数的凸面积介于L0和L2之间（图 1b），常作为L0范数的最优凸近似，比L0范数更加容易优化求解。L2范数凸面积最大（图 1c），稳定性高，最易求解。从反演精度和计算耗时两方面进行综合考虑，低分辨率大尺度信号选用同样精度低但计算快的L2范数约束；中尺度信号稀疏性一般，集中了最多峰值频率信号，选用L1范数约束；小尺度信号信噪比最低，存在大量有效的高频波阻抗信息，选用L0范数约束。

由于目标函数中存在多组正则化参数，联合L曲线法和试错法确定最优的正则化参数组合，以获取准确性高且分辨率高的反演结果。目标函数[式（18-A）]通过共轭梯度法计算；[式（18-B）]中存在L1正则化项，采用迭代重加权最小二乘法求解 ^[29]；[式（18-C）]中存在L0正则化项，采用分裂布雷格曼方法 ^[30]。

首先，对叠后地震数据进行MP-VMD谱分解，根据时频能量分布将原始地震数据重构成大、中和小尺度信号；其次，利用3个不同尺度测井数据分别约束3个尺度的地震信号，从大尺度项目标函数向小尺度目标函数项进行逐级约束求解；最后，地震反演结果以小尺度项目标函数求解结果为准。具体流程如图 2所示。

2 模型数据测试 2.1 多方法时频能量测试对比

根据地震信号的多尺度特征，构建多级旋回复合模型（图 3a），模型中存在2套大尺度旋回体、2套中尺度旋回体、2套小尺度旋回体。根据褶积理论，利用30 Hz雷克子波合成地震记录（图 3b）。

由于地震波在传递过程中会存在调谐效应，难以直接从地震信号中直接剥离出不同旋回体对应的地震响应，而高精度时频分析方法能有效刻画地震信号中各时刻频率能量特征，进而提高不同旋回体的识别精度，因此，时频能量的精度直接决定了分频结果的好坏。建立MP-WVD时频谱（图 3c），将合成记录分为大尺度旋回体信号（频率为6~ 25 Hz）、中尺度旋回体（频率为25~40 Hz）及小尺度旋回体（频率为40~75 Hz）；MP-WVD时频结果的能量团时间宽度与波阻抗模型时间域宽度一致，这也说明MP-WVD时频谱的可靠性。

与连续小波变换时频谱、S变换时频谱相比，MP-WVD时频谱中能量集中性较强，时间和频率方向的分辨率都更高；与WVD变换时频谱相比，MP-WVD时频分析是通过匹配追踪进行原子分解后，再进行WVD变换，有效克服了交叉项的干扰（图 4）。因此，MP-WVD时频谱更适用于地震信号分频处理。

2.2 正则化参数测试

地震信号的时频分析能够有效揭示地震信号中各个频率段对应的信号，并加强其重构的性质，对地震信号中的高低频率信号进行分频处理，最终获得地震信号分频体。基于MP-VMD时频分解和重构原理，对合成记录时频谱高、中和低频率范围的小波原子进行叠加，得到多尺度地震信号分频结果。其中，大尺度信号分频结果能够明显反映模型中的大套厚层地层响应特征，地震波长较长；厚度中等地层的地震响应主要集中在中尺度地震信号分频体中，地震波长适中；厚度较小的地层主要对应地震信号的高频部分，与小尺度地震信号对应较好，地震波长较短（图 5）。

地震反演结果的精度很大程度上取决于地震子波构建，不同尺度地震信号的子波特征存在明显差异。大尺度地震统计子波旁瓣能量较弱，统计子波的主频为20 Hz；中尺度信号中提取的地震统计子波旁瓣能量有所增加，其主频为30 Hz；小尺度信号中提取的地震统计子波旁瓣能量较强，其主频约为47 Hz（图 6）。基于以上分析可知，不同尺度地震信号对应的统计子波均与合成信号使用的30 Hz雷克子波存在差异，因此在反演过程中不能直接使用合成信号中的理论子波，需重新从地震信号中提取地震子波 ^[16]。

基于Hamid等 ^[31]多正则化约束项的权重系数确定思想，联合L曲线法和试错法确定正则化约束项前的权重系数。通常，背景模型约束项的正则化参数过大，会导致反演结果过于平滑，且产生背景模型化，而过小又会导致计算不稳定。当目标函数中包含其他稀疏项时，反演过程不仅依赖背景模型约束项来提升计算的稳定性，还需为背景项设置较小的正则化参数值，进而研究其他2个主导约束项[如式（18-C）中$ \left\|\boldsymbol{Z}_{\mathrm{s}}-\boldsymbol{Z}_{\mathrm{m}}\right\|_{2}^{2} $和$ \left\|\boldsymbol{D} \boldsymbol{Z}_{\mathrm{s}}\right\|_{0} $]的正则化参数取值。以式（18-C）为例来说明式（18）中的正则化参数确定过程，在给定一个较小$ \lambda_{s} $情况下，令$ \beta_{\mathrm{s}}=0 $，用L曲线法确定$ \gamma_{\mathrm{s}} $最优取值范围，即不考虑L0范数的约束，在$ \gamma_{\mathrm{s}} $最优取值范围内，固定$ \gamma_{\mathrm{s}} $取值，继续利用L曲线法寻找最优的$ \beta_{\mathrm{s}} $取值，最后根据反演结果合成记录与地震记录的均方根误差（RMSE）最小的原则，可确定最优$ \gamma_{\mathrm{s}} $和$ \beta_{\mathrm{s}} $。当$ \lambda_{\mathrm{s}}= $ $ 10^{-3} $时，由$ \left\|\boldsymbol{Z}_{\mathrm{s}}-\boldsymbol{Z}_{\mathrm{m}}\right\|_{2}^{2} $前的正则化参数L曲线变化（图 7a）可知，$ \gamma_{\mathrm{s}} $最优取值为$ 5 \times 10^{-4} \sim 5 \times 10^{-3} $，将$ \gamma_{\mathrm{s}} $分为10组，第1组$ \gamma_{\mathrm{s} 1}=5 \times 10^{-4} $，由$ \left\|\boldsymbol{D} \boldsymbol{Z}_{\mathrm{s}}\right\|_{0} $前的正则化参数L曲线变化（图 7b）可知，最优$ \beta_{\mathrm{s} 1}=0.037 $。按照$ \gamma_{\mathrm{s} 2}, \ldots \gamma_{\mathrm{s} 10} $依次为$ 6 \times 10^{-4}, 7 \times 10^{-4}, 8 \times 10^{-4}, \cdots $，$ 5 \times 10^{-3} $等9个值，得到最后一组正则化参数L曲线变化图，可知$ \gamma_{\mathrm{s} 10}=5 \times 10^{-3} $时，$ \beta_{\mathrm{s} 10}=0.089 $，而RMSE最小的测试组为第6组。因此，式（18-C）中的各正则化参数为$ \lambda_{\mathrm{s}}=1 \times 10^{-3}, \gamma_{\mathrm{s}}=1 \times 10^{-3}, \beta_{\mathrm{s}}=0.036 $（表 1）。

在大尺度分频信号地震反演中，由于采用的低频模型精度较低，且目标函数约束使用了稀疏性较弱的L2范数，导致反演结果仅能反映模型变化的整体趋势，而无法准确描述薄层与中厚层特征（图 8a）。相比之下，中尺度分频信号的地震反演通过引入稀疏性较强的L1范数，使得反演结果呈现出明显的块状特征，不仅提升了厚层的识别精度，还准确刻画了中频信号中的中厚层特征（图 8b）。小尺度分频信号地震反演采用稀疏性最强且具备较强抗噪性的L0范数进行约束，这种方法能够充分挖掘小尺度信号中的高频信息，如薄层特征（图 8c）。小尺度反演结果是建立在大尺度和中尺度信号反演精度的基础上，因此，地震分频多级稀疏正则化反演结果以小尺度反演结果为准。

2.3 同类反演方法测试对比

稀疏脉冲反演由于其算法具有较强的稳定性，已成为地震储层预测中最常用的反演方法之一。将地震分频多级稀疏正则化反演与稀疏脉冲反演进行对比（图 9），2种方法得到的反演结果趋势基本一致，稀疏脉冲反演稳定性虽好，但分辨能力弱，无法识别薄层信息（如图 9中黑色箭头所示），而地震分频多级稀疏正则化反演结果精度和准确性更高。

3 实际应用效果

以渤中凹陷石臼坨凸起构造带A区为例，验证本文所提出的地震分频多级稀疏正则化反演方法的实用性。A区位于石臼坨凸起南部（图 10a），具有良好的油气运聚条件和巨大的油气勘探潜力，古近系底受北西向和近东西向断裂带控制，新近系以来受北西向和近东西向2组应力控制形成了较复杂的断裂系统 ^[32-33]。该构造带中潜山为中生界、古生界与太古界形成的大背斜，古近系是一个继承性发育披覆背斜，构造高部位及斜坡区以东营组沉积为主，沙河街组主要分布在边界断层下降盘，受古地貌和沟槽的影响，在古近系形成了以扇三角洲为主的近源沉积体系，古近纪末期的构造抬升使得构造高部位地层剥蚀严重（图 10b）。

古近系自下而上发育孔店组、沙河街组和东营组，其中东营组可自下而上分为东三段、东二段和东一段，本文主要研究层段东二段可细分为东二上亚段和东二下亚段，其中东二上亚段岩性主要为砂岩以及泥岩，东二下亚段岩性为厚层泥岩与砂岩、粉砂岩互层。东二下亚段与其下伏东三段（主要发育大套厚泥岩层，夹有少量薄层砂岩、粉砂岩）为三角洲和中深湖沉积，泥岩中有机质丰度高、类型好，具有较高生烃潜力。东营组上覆馆陶组沉积了大套泥岩，作为盖层（图 10b）。

研究区东二段对应地震反射层为$ \mathrm{T}_{2}-\mathrm{T}_{2 \mathrm{\; m}} $段，地震资料信噪比较低，存在大量杂乱弱反射（图 11），地震资料主频为21 Hz，有效频带宽度为$ 7 \sim 48 \mathrm{\; Hz} $，整体上主频较低，有效频带窄（图 12）。以区域内A井与B井为例，2口井的直线距离约6.3 km，由测井解释可知，2口井东二段岩性横向变化大，中深部的A井以发育水层为主，浅部的B井发育一套9 m厚的含油水层。A井与B井的时频谱整体差异较大，但沿时间方向累加得到的时频能量随频率变化曲线差异较小，主频带均为$ 11 \sim 24 \mathrm{\; Hz} $（图 13，14）。

绘制图 11中所有道的时频能量随频率变化曲线，得到二维剖面的时频能量分布（图 15），根据其变化可以将地震信号分为3个尺度，大尺度信号对应频带$ 0 \sim 10 \mathrm{\; Hz} $；中尺度信号对应频带$ 11 \sim 22 \mathrm{\; Hz} $，具有较强的能量，地震资料主频也主要集中在这一频率范围内；小尺度信号对应的频带$ 23 \sim 50 \mathrm{\; Hz} $，时频能量较弱，主要是包含地震信号中的一些高频弱反射和噪音的时频能量。

在时频谱能量划分尺度范围内，选取符合频率范围的小波原子，利用MP算法对不同尺度的地震信号进行分频重构，得到地震资料分频结果和统计子波提取结果。结果（图 16）显示：不同尺度地震信号同相轴连续性均较好；与原地震波形剖面（参见图 11）相比，3个尺度的地震波形剖面的信噪比均有所提升；大尺度地震信号的分辨率较低、信噪比较好，小尺度地震信号的分辨率最高，地震资料中强反射轴能量有所压制，弱信号能量得到提升，但弱信号的同相轴连续性较差，存在大量的随机噪音；3个尺度地震信号的统计子波长度均为64 ms，采样间隔均为1 ms。

在反演测试中，以B井作为地震反演约束井，A井作为反演验证井。将采用多尺度稀疏约束反演方法得到的大、中、小尺度波阻抗反演结果与单井测井大、中、小尺度的声波阻抗曲线相结合，结果（图 17）显示：不同尺度的反演结果中强波阻抗和弱波阻抗分布位置大致相同；大尺度波阻抗反演结果较稳定可靠，能够较真实地反映砂泥地层波阻抗背景，但仅能刻画地层中的大套砂泥岩体，不能反映中薄层砂体的波阻抗特征；中尺度波阻抗反演结果分辨率有所提升，能够比较明显地刻画储层主体分布；小尺度波阻抗反演结果能够充分利用与挖掘高频信息，精确刻画薄砂体分布特征。将这3个尺度的波阻抗反演结果与A井和B井的对应尺度测井声波阻抗曲线进行对比，井旁波阻抗强度变化与测井曲线变化趋势基本一致，表明不同尺度波阻抗反演结果准确度均较高。A井和B井周围的波阻抗反演体显示出较差的横向连续性，这也说明了石臼坨凸起东二段砂体横向变化快，储层非均质性较强，与现有的工区地质特征认识一致，也证实了地震分频多级稀疏正则化反演方法的可靠性。

通过地震分频多级稀疏正则化反演得到的A井周围不同尺度的波阻抗结果，与相应尺度的地震波形进行叠合。波阻抗的横向变化特征与对应尺度的地震波形的横向变化特征基本一致（图 18a-18c），进一步验证了该方法在波阻抗反演中始终忠实于地震波形。

对研究区分别进行地震分频多级稀疏正则化反演与稀疏脉冲反演，结果（图 18c，18d）显示：稀疏脉冲反演结果虽能和地震波形变化保持一致，但整体分辨率低，介于地震分频多级稀疏正则化约束的大尺度和中尺度反演结果之间；地震分频多级稀疏正则化反演能够充分挖掘不同尺度信号的频率信息，尤其能够充分利用高频弱信号特征，反演结果具有较高的纵向分辨率，对薄层刻画较准确。

4 结论

（1）基于地震分频的多级稀疏正则化反演方法，通过匹配追踪-Wigner-Ville分布技术将地震信号有效分解为大、中、小3个尺度的频段，并在多尺度下应用不同的正则化约束（L2，L1，L0范数），显著提升了薄储层的刻画精度，能够准确挖掘各频段的潜地层特征信息，确保反演结果的可靠信性。

（2）多级稀疏正则化反演方法有效解决了反演过程中不同尺度地震信号的尺度差异问题，通过逐级迭代策略，增强了反演结果的分辨率和抗噪性，单井反演结果与对应的测井声波阻抗曲线吻合度高。

（3）通过渤中凹陷石臼坨凸起构造带A区实际资料应用验证，与传统的稀疏脉冲反演方法相比，地震分频多级稀疏正则化反演对薄储层的识别能力更强，能够准确揭示研究区储层高非均质性和快速横向变化特征。