0 引言
地震反演能够从地震数据中得到速度和密度等参数,是岩性识别的重要技术[1]。传统的模型驱动地震反演方法通常以Robinson褶积模型理论为基础[2],从参数初始模型出发,以合成地震数据与实际地震数据之间的误差迭代更新参数模型[3-4],主要包括波阻抗反演和AVO反演[5-6]。相比于波阻抗反演利用叠后地震数据,AVO反演利用叠前地震数据,能够提供精度、分辨率更高的速度和密度参数[7],因此被广泛地应用并取得了诸多成果[8]。但此类反演方法通常很难突破地震数据频带带宽的限制[9-10],需要准确的初始低频模型,且其反演结果的分辨率通常不高[11⇓⇓-14]。基于测井约束的模型驱动反演方法将具有高垂向分辨率的测井数据和具有大空间展布范围的地震数据相结合,能够在保证反演结果精度的同时提高反演结果的分辨率[15⇓-17]。但这种方法需要较多的测井数据,且要求井位分布均匀,适合处理构造相对简单的目标储集层,难以用于小样本复杂砂泥岩薄互层的精细识别[18⇓-20]。
为了能在小样本数据的情况下,提高速度、密度反演结果的精度和分辨率,实现薄互层砂体的精细识别,本文提出了一种基于空变目标函数的模型-数据驱动（双驱动）地震AVO反演新方法。该方法利用神经网络拟合地震数据与地震参数之间的强非线性关系,以预测参数与样本参数之间的误差、预测参数低频分量与初始低频模型之间的误差、合成地震数据与实际地震数据之间的误差构建目标函数,以反距离加权思路根据反演目标道所在的位置控制这3个误差在目标函数中的权重,以反向传播算法迭代训练神经网络权系数,实现速度、密度等参数的反演,并在丽水凹陷开展实际数据测试,以期在小样本数据的情况下,得到较高精度、较高分辨率的薄互层砂体的速度和密度参数,实现薄互层砂体的精细识别。
1 反演方法原理
1.1 模型驱动地震AVO反演方法
模型驱动地震AVO反演方法通常以Robinson褶积模型为理论基础[2]。在不考虑随机噪音的情况下,褶积模型可以表示为：
$S=GM* W$
假设G可逆,且不考虑地震子波的影响,反演参数 $M_{inv}$可以表示为
$M_{inv}=G^{-1}S_{real}$
但在AVO反演中,G通常是不可逆的,需要引入其伴随矩阵进行求解,此时反演参数 $M_{inv}$表示为
$M_{inv}=G^{H} G^{-1} G^{H}S_{real}$
AVO反演方法通常利用最小二乘优化理论使合成地震数据与实际地震数据之间的误差达到最小,实现速度、密度等参数的反演,其目标函数为：
$J_{1}=||GM_{inv}-S_{real}||_{F}$
其中 $||·||_{F}$表示F范数（Frobenius范数）,其表达式为：
$\|A\|_{\mathrm{F}}=\left(\sum_{i=1}^{l} \sum_{j=1}^{k} a_{i j}^{2}\right)^{1 / 2}$
由于模型驱动AVO反演方法主要利用地震数据,受限于后者的频带宽度,其反演结果的分辨率普遍不高。
1.2 数据驱动地震预测方法
随着神经网络技术的迅速发展,以神经网络为基础的数据驱动方法被成功应用于速度、密度等地震参数的预测。这类方法通常首先利用测井数据和井旁道地震数据制作训练样本;然后将样本地震数据输入构建好的神经网络中,输出预测参数;接着基于预测参数与样本参数之间的误差迭代更新神经网络权系数（神经网络训练）;最后将实际地震数据输入训练好的神经网络中,输出预测结果。神经网络训练过程中的目标函数一般定义为：
$J_{2}=\left\|\boldsymbol{M}_{\mathrm{pre}}-\boldsymbol{M}_{\text {sample }}\right\|_{\mathrm{F}}$
由于数据驱动地震预测方法主要利用由测井数据生成的训练样本,其预测结果的分辨率较高,但该方法需要大量有代表性的训练样本。受限于实际生产中测井数据的数量,这类方法只在井周围或与测井数据特征类似的区域才能得到预测精度较高的速度和密度参数。
1.3 模型-数据驱动地震AVO反演新方法
在实际生产中进行砂泥岩薄互层识别时,当勘探面积较大的工区内只有几口井时,很难提供大量有代表性的训练样本,属于小样本数据算例。为了能在这种情况下得到较高精度、较高分辨率的薄互层砂体的速度和密度参数,进而实现薄互层砂体的精细解释,本文提出了一种基于空变目标函数的模型-数据驱动（双驱动）地震AVO反演新方法,其流程如图1 所示。在反演之前需要利用测井数据和地震数据等制作初始低频模型和训练样本,并通过在测井数据中添加随机数的方式对训练样本进行增广。新方法为多道预测,每次将相同道数的实际地震数据 $S_{real}$和样本地震数据 $S_{sample}$（基于褶积模型理论,利用测井数据计算得到）输入构建好的神经网络中,输出对应的样本预测参数Mpre-sample和真实预测参数Mpre-real。然后将Mpre-real滤波为与初始低频模型 $M_{initial}$具有同样频带宽度的Mpre-real-filter。同时基于褶积模型,利用Mpre-real正演合成地震数据 $S_{syn}$。最后以Mpre-sample与 $M_{sample}$的误差、Mpre-real-filter与 $M_{initial}$的误差、 $S_{syn}$与 $S_{real}$的误差建立目标函数,并利用反向传播算法迭代更新神经网络权系数直到目标函数趋于收敛。此时神经网络输出的真实预测参数Mpre-real即为反演结果。这个过程中的目标函数可以表示为：
$\begin{aligned}J_{3}=& \mu\left\|\boldsymbol{M}_{\text {pre-sample }}-\boldsymbol{M}_{\text {sample }}\right\|_{\mathrm{F}}+\eta\left\|\boldsymbol{M}_{\text {pre-seal-filer }}-\boldsymbol{M}_{\text {initial }}\right\|_{\mathrm{F}}+\\&(1-\eta-\mu) c\left\|\boldsymbol{S}_{\mathrm{syn}}-\boldsymbol{S}_{\text {real }}\right\|_{\mathrm{F}}\end{aligned}$
（7）式右侧的第1项表示训练样本对目标函数的约束。由于训练样本通常是由测井数据制作的,其在井周围对反演的正面影响较大。因此,越靠近井位,权系数μ值越大。笔者利用反距离加权原理来定义μ,其表达式为：
$\mu=\sum_{p=1}^{n} \frac{1}{d_{p}} / \sum_{q=1}^{m} \frac{1}{d_{q}}$
（7）式右侧第2项表示初始低频模型对目标函数的约束。初始低频模型质量越好,权系数η值越大。需要注意的是 $\mu \ge 0$、 $\eta \ge 0$且 $\mu +\eta \le1$。本文将 $\mu +\eta $赋值为一个小于1的常数（模型数据为0.6;实际数据为0.8）,然后根据μ的值来确定η的值。
（7）式右侧第3项表示地震数据对目标函数的约束。由于地震数据具有波形特征,当实际地震数据与合成地震数据之间存在振幅差异时,基于F范数的目标函数容易陷入局部极小值,影响反演结果的精度。基于零延迟互相关的目标函数考虑了波形之间的相位拟合,可以更好地处理这种振幅差异。此外,与F范数相比,零延迟互相关构建的目标函数能够更好地压制随机噪声。因此,笔者将（7）式改进为：
$\begin{aligned}J_{4}=& \mu\left\|\boldsymbol{M}_{\text {pre-sample }}-\boldsymbol{M}_{\text {sample }}\right\|_{\mathrm{F}}+\eta\left\|\boldsymbol{M}_{\text {pre-real-filter }}-\boldsymbol{M}_{\text {initial }}\right\|_{\mathrm{F}}+\\&(1-\eta-\mu) c\left(1-\frac{\boldsymbol{S}_{\mathrm{syn}} \odot \boldsymbol{S}_{\text {real }}}{\sqrt{\boldsymbol{S}_{\mathrm{syn}} \odot \boldsymbol{S}_{\mathrm{syn}}} \sqrt{\boldsymbol{S}_{\text {real }} \odot \boldsymbol{S}_{\text {real }}}}\right)\end{aligned}$
新方法采用如图2 所示的6层U-net神经网络,包括下采样过程和上采样过程。每个下采样层由3×3卷积层和2×2池化层组成,每个上采样层由2×2上卷积层和3×3卷积层组成。神经网络的输入为3个角度叠加地震数据,输出为纵波速度、横波速度和密度。此外,反演之前需要对所有样本数据进行归一化处理,以保证其处于同一数量级,便于神经网络权系数的更新。归一化公式为：
2 双驱动反演方法试算
为了分析模型-数据驱动（双驱动）地震AVO反演新方法反演结果的精度、分辨率及其对砂泥岩薄互层的识别能力,建立了薄互层地质模型进行反演试算。
薄互层地质模型如图3 所示,共包括201个CDP（共深度点道集）,泥岩背景中发育了4组砂体,第1组为3个厚度为3 m的砂岩层,泥岩夹层厚度为3 m,后面3组分别去掉其中一个砂岩层。其中,砂岩的纵波速度、横波速度和密度分别为3 400 m/s、2 000 m/s和2.65 g/cm3,泥岩的纵波速度、横波速度和密度分别为2 800 m/s、1 500 m/s和2.35 g/cm3。在该模型上建立了4口井,分别代表不同的储集层特征,依次命名为1号井、2号井、3号井和4号井。
实际地震数据的有效频带宽度是有限的,为了真实地模拟实际数据,笔者基于该地质模型的速度、密度,利用Zoeppritz方程计算反射系数[26],并与3,8,32,37 Hz的带通子波褶积合成叠前地震数据。在实际地震反演中,通常会将叠前地震数据按照小、中、大角度进行部分角度叠加作为反演的输入以提高地震数据的信噪比和反演效率。叠加的原则是既覆盖整个角度范围,又包括不同角度区域的地震特征。因此,笔者将1°~15°、16°~30°和31°~45°的叠前地震数据分别叠加为8°、23°和38°叠加地震数据（见图4 ）作为神经网络的输入。实际生产中,测井数据通常较少,一般利用地震数据和两口井的测井数据制作初始低频模型（见图5 ）。此外,为了验证新方法的稳定性,笔者分别利用1号井和3号井的测井数据生成两组训练样本（小样本数据）,并通过在测井数据中添加随机数的方式进行样本增广以用于反演试算。
首先,利用基于Zoeppritz方程的模型驱动方法对该地质模型进行反演,结果如图6 所示。该方法的反演结果具有与地震数据近似的分辨率,无法识别薄互层砂体。然后,利用新方法对该数据进行反演,基于1号井和3号井训练样本的反演结果分别如图7 、图8 所示。上述两种情况下,新方法反演结果的分辨率基本相同,都能够清晰识别出砂泥岩薄互层,表明新方法在小样本数据的情况下能够较稳定地反演得到精度和分辨率较高的速度和密度参数。
为了分析新方法的抗噪能力,在图4 地震数据中加入随机噪声生成信噪比为2的含噪地震数据（见图9 ）,然后利用新方法对含噪地震数据进行反演,结果如图10 所示。由图可见,在地震数据含噪的情况下,新方法的反演结果仍然具有较高的精度和分辨率,能够识别3 m厚度（约1/30波长）的薄互层砂体,表明新方法具有较好的抗噪能力。
3 实例应用
3.1 丽水凹陷储集层特征
研究区位于东海陆架盆地的丽水凹陷,600 km2工区内只包含3口井（A井、B井和C井）,属于典型的小样本数据算例。目的层段古新统深度范围为3 800~4 200 m,呈砂泥岩薄互层特征,单砂体厚度最小为3 m,需要精细识别[27],地震数据主频约为23 Hz,最高有效频率约为41 Hz。利用测井、录井数据开展岩性测井解释,建立砂泥岩解释图版（见图11 ）。由图可见,纵波波阻抗能够识别砂岩（砂岩的纵波波阻抗大于9 300 g·m/（cm3·s））,可以开展基于速度和密度反演结果的薄互层砂体识别。当在某些实际工区内,纵波波阻抗无法识别砂岩时,可以利用速度和密度组合构成其他参数来区分砂岩和泥岩。
3.2 不同反演方法的反演结果对比
本次选取研究区内过A井、B井和C井的一条测线进行试算,共包括801个CDP。反演之前对叠前地震数据进行了去噪处理,以提高其信噪比。8°、23°、38°叠加地震数据和叠后地震数据如图12 所示。以图12 中M、P、N层作为约束,利用A井和C井的测井数据插值得到初始低频模型（见图13 ）。以A井、C井的测井数据和井旁道地震数据制作训练样本并进行增广,以B井作为验证井。
利用基于Zoeppritz方程的模型驱动方法和本文新方法对实际数据进行反演（见图14 、图15 ）。由图可见,新方法反演结果的分辨率高于模型驱动方法反演结果。这是因为模型驱动方法主要利用了地震数据的信息,其反演结果的分辨率受到地震数据分辨率的限制。由于实际地震数据缺乏有效的高频信息,导致该方法的反演结果分辨率不高。而新方法更多利用了来自测井数据的有效高频信息,能够突破地震数据分辨率的限制,进而提高了反演结果的分辨率。
为了清晰对比两种方法反演结果的精度,笔者提取了A井（训练井）和B井（验证井）的数据（见图16 ）,并分别计算了反演结果与测井数据之间的相对误差。由图可见,对于A井,新方法的纵波速度反演结果（相对误差为0.023）、横波速度反演结果（相对误差为0.035）和密度反演结果（相对误差为0.037）与测井数据的吻合度高于模型驱动方法的反演结果（相对误差分别为0.152,0.178和0.218）,表明新方法在训练井处的反演结果具有较高的精度;对于B井,新方法的反演结果（相对误差分别为0.025,0.038和0.040）比模型驱动方法的反演结果（相对误差分别为0.161,0189和0.253）更接近测井数据,表明新方法在验证井处的反演结果精度也高于常规方法,具有一定优势。这是因为模型驱动方法以合成地震数据和实际地震数据之间的误差作为目标函数,当该目标函数满足迭代终止条件时,能够保证合成地震数据与实际地震数据之间的误差收敛于极小值,但不能保证反演参数与真实值之间的误差收敛于极小值;而新方法在目标函数中考虑了预测参数与样本参数之间的误差、预测参数低频分量与初始低频模型之间的误差、合成地震数据与实际地震数据之间的误差,当目标函数满足迭代终止条件时,能够保证这3个误差均较小,反演得到精度较高的纵波速度、横波速度和密度。
此外,笔者还以图12 中的P层位为中心,以20 ms的时窗提取了两种方法反演的纵波速度、横波速度和密度的平均值,并将其展示在图17 中。由图可见,新方法的反演结果具有更高的分辨率,其空间展布情况与已知的地质解释结果相匹配,表明新方法具有较好的实用性。
3.3 双驱动反演方法刻画薄互层砂体
利用模型驱动方法和新方法反演得到的纵波速度和密度计算纵波波阻抗（见图18 ）,由图可见,模型驱动方法得到的纵波波阻抗分辨率不高,只能识别大套的砂体;新方法得到的纵波波阻抗分辨率较高,可以识别薄互层砂体,同时,砂体横向变化特征符合地震特征。为进一步分析两种方法的反演结果对薄互层砂体的刻画能力,提取B井的纵波波阻抗,并展示在图19 中。由图可见,新方法解释岩性与测井解释岩性吻合度高,在4 000 m深度左右可识别最小厚度为6 m（约1/15波长）的砂岩层,表明新方法在小样本数据的情况下能够较精细地识别薄互层砂体。此外,笔者还提取了图12 中M—N层的砂岩厚度图,并将其展示在图20 中。由图可见,A井、B井附近的砂岩厚度较大,C井附近的砂岩厚度较小。这与测井解释结果吻合较好,表明了新方法的可行性和有效性。
4 结论
针对小样本数据,本文提出了一种基于空变目标函数的模型-数据驱动（双驱动）地震AVO反演新方法。该方法利用零延迟互相关和F范数构建目标函数,以反距离加权理论根据反演目标道的位置控制训练样本、初始低频模型和地震数据对目标函数的约束权重,能够在小样本数据的情况下,反演得到较高精度、较高分辨率的速度和密度等参数,可用于薄互层砂体的精细识别。薄互层地质模型试算结果表明,基于小样本数据和合适的初始低频模型,新方法对薄互层砂体的识别精度达到了约1/30波长,且具有较好的稳定性和抗噪能力。丽水凹陷实例应用表明,在小样本数据和合适的初始低频模型情况下,新方法的反演结果具有较高的反演精度和分辨率,能够较为准确地识别薄互层砂体,实现精细化解释。
符号注释：
A——待处理矩阵;aij——待处理矩阵的元素;c——将公式（7）右侧第3项和前两项校正到同一数量级的系数,无因次;dp——实际地震数据的中心道距离第p口井的距离,无因次;dq——实际地震数据的中心道距离第q道的距离,无因次;G——映射矩阵;G-1——G的逆矩阵;GH——G的伴随矩阵;i——待处理矩阵的行序号;j——待处理矩阵的列序号;J1——模型驱动方法的目标函数;J2——数据驱动方法的目标函数;J3——新方法的目标函数;J4——新方法的改进目标函数;k——待处理矩阵的列数;l——待处理矩阵的行数;M——速度、密度构成的参数矩阵;Minitial——初始低频模型;Minv——反演参数矩阵;Mpre——预测参数矩阵;Mpre-sample——样本预测参数矩阵;Mpre-real——真实预测参数矩阵;Mpre-real-filter——真实预测参数低频分量矩阵;Msample——样本参数矩阵;m——工区内的道数; $\max |x|$——样本数据绝对值的最大值;n——井的数量;S——地震数据矩阵;Sreal——实际地震数据矩阵;Ssample——样本地震数据矩阵;Ssyn——合成地震数据矩阵;W——地震子波矩阵;x——样本数据;xscale——样本数据归一化后的值;μ——（7）式右侧第1项的权系数,无因次;η——（7）式右侧第2项的权系数,无因次; $\odot$——矩阵内积; $\|\cdot\|_{F}$——F范数。