基于改进的相位目标函数估算薄层厚度

王 鹏 胡向阳 魏水建

(中国石化石油勘探开发研究院，北京 100083)

1 引言

在地震波传播过程中，地下介质的岩性或物性变化都将引起地震波属性特征的变化，包括速度、振幅、频率、相位等，这些变化已经成为地球物理储层预测、流体识别的重要依据[1]。相位作为地震记录的重要属性之一，携带了丰富的地层反射信息，如地层反射结构特征等。

相位属性在地震解释中有较好的应用，尤其是在断裂解释中。张应波[2]认为地震相位表现了丰富、形象的地质现象，并能解释波形变面积剖面无法解释的多种沉积现象； 毕俊凤[3]、柏冠军等[4]通过观察不同频率的相位滞后分布，较好地实现了宽频带地震资料的断层精细解释和断裂系统识别，获得了较高的分辨率； 王庆华[5]利用分频突出断层响应的最优势频带和相位对地层产状变化敏感的特点，开发了分频相位技术，克服了常规相干技术的不适应性，在识别断层分布特征上起到了独特的作用；姬战怀等[6,7]在分频剖面上利用瞬时相位识别小断层； Kazmi等[8]提出使用连续相位谱的曲率属性识别断层特征的方法，克服了现有方法的局限性并提供了更高的体曲率属性分辨率； Alam等[9]提出了一种利用相位谱对三维数据体中断层进行自动检测的技术； 刘道平[10]、赵子豪等[11]利用地震记录的相位属性进行工区地震层序标定，使砂体标定解释更加可靠； 赵淑红等[12]借鉴图像处理方法，研究利用地震记录的相位谱进行数据恢复。

在薄层识别和厚度估计中，人们一般以地震数据的振幅谱为基础，利用异常特征(如零点周期[13,14]、调谐厚度[15]、振幅谱梯度[16]和主频变化[17-19]等)开展研究。高静怀等[20]、范明霏等[21]定义了广义S变换并将其应用于薄层的地震探测，地层厚度识别能力从λ/4提高至λ/12； Barnes[22]通过分析薄层的振幅谱响应特征研究了频率变化与薄层厚度之间的关系； Puryear等[23]将谱反演应用于检测储层厚度及沉积特征； Nowak等[24]提出了基于AVO响应的地层厚度定量预测方法，提高了谱分解方法对薄层的检测能力； 田鑫等[25]利用地震属性方法进行薄层砂体识别及有效厚度预测。

上述方法虽取得了一定的成果，但均很难规避顶底反射系数组合和地震子波振幅谱的影响。相比于振幅信息，相位对地层厚度的变化同样较为敏感，但由于相位谱较复杂且对地层的响应机制不够明确，相应的研究不多。Georgy等[26]将相位谱应用于薄互层尖灭特征分析，与振幅谱相比，相位谱能够为尖灭点确定提供更准确的信息； 蔡涵鹏等[27]在地震频带范围内应用瞬时相位谱构建了地层厚度估算的目标函数，将顶底反射系数组合项与地层厚度项巧妙地分离，在地层估算过程中不需考虑反射系数大小、极性或地震子波主频的影响。该方法具有较好的理论意义，但存在两个问题： 正切函数有奇异值，影响了目标函数的求解； 受噪声污染的影响较为严重，在一定程度上限制了该方法的应用。

事实上，相位谱具有振幅谱不具备的优势。例如，不论反射振幅强或弱，其相位角基本维持在同一量级；地震记录振幅谱是地震子波振幅谱与反射系数振幅谱的乘积，而前者的相位谱则是后两者相位谱之和，求和与乘积相比相对简单。本文基于双反射系数组合和地震子波常相位的假设，推导了新的地震记录相位与地层厚度关系式，并用该式进行薄层厚度的估算，与已有的相位—厚度关系式相比，新关系式规避了相位正切函数存在的奇异点问题，提高了该方法的抗噪性。模型测试及实际资料处理结果表明，该方法可以有效提高利用相位估算薄层厚度的可靠性。

2 基本原理

假设地层反射是2个反射系数，顶、底反射系数分别为r1和r2，地层厚度为Δt，双程旅行时为2Δt，反射系数所在时刻分别为t-Δt和t+Δt，地震子波频谱为W(f)，复合地震记录频谱为S(f)，不考虑地层吸收、速度频散、入射角度等因素，根据褶积模型，则有

S(f)=W(f)[r1e-i2πf(t-Δt)+r2e-i2πf(t+Δt)]

(1)

假设地震子波相位是零相位，时刻t=0，将式(1)右项改写为实部U(f)与虚部D(f)之和，则实部和虚部分别为

U(f)=W(f)(r1+r2)cos(2πfΔt)

(2)

D(f)=W(f)(r1-r2)sin(2πfΔt)

(3)

式(2)与式(3)相比，可得相位谱φ(f)

(4)

求相位谱的正切值，得

(5)

Ω(f,τ)=tan[φ(f)]·cot(2πfτ)

(6)

τ为扫描变量，当τ=Δt时，中间函数Ω(f,τ)=R为一恒值(标准差为0)，当τ≠Δt时，中间函数Ω(f,τ)则是频率f的函数，不再是恒值(标准差不为0)。

以中间函数的标准差构建目标函数Λ(τ)

(7)

(8)

(9)

(10)

3 仿真数据测试

3.1 奇异值影响分析

假设有一地层处于两个半无限介质中间，双程旅行时为5ms，顶、底反射系数极性均为正，比值为2，地震子波是50Hz主频的Ricker子波(薄层厚度约为λ/8)，时间采样间隔为0.5ms，频率采样间隔为0.5Hz。图1a是反射系数(红色虚线)及其合成地震记录(蓝色实线)的归一化显示。图1b是合成记录的振幅谱，由于受地震子波本身的带宽限制，150Hz以上的振幅较弱，可用信息在0～150Hz之间，且振幅谱上基本看不到反射系数振幅谱的痕迹。图1c是合成记录的相位谱，在0～300Hz之间均有值显示，尤其能够清晰地刻画150～300Hz之间的信息，不受弱振幅的影响，换言之，相位谱可利用信息的频带范围比振幅谱宽。图1d是合成记录频谱的实部(红色虚线)和虚部(蓝色实线)。

图1 合成记录及其频域变换 (a)反射系数及其合成地震记录； (b)合成记录的振幅谱； (c)合成记录的相位谱； (d)合成记录频谱

图2a是根据式(6)估算的中间函数，不同曲线取的τ值不同。可以看到，当τ值与准确值不相等(红线或紫线)时，曲线是非平稳的，其值随频率有一定的变化；当τ值与准确值相等(蓝线)时，曲线呈一恒定值，满足式(7)最佳极值逼近的假设。但蓝线中还存在奇异值及其附近的异常点，这些值突然增大，有可能影响最佳极值的逼近。图2b是根据式(9)估算的中间函数，不同曲线取的τ值不同(同图2a)。可以看到，当τ值与准确值不相等(红线或紫线)时，曲线值较大且波动较强；当τ值与准确值相等(蓝线)时，曲线在零值附近波动，幅度明显变小，需要说明的是，按照理论计算，曲线值应为零，但由于数值离散计算损失了一定的精度，曲线存在一定的误差是合理的，也是数值计算过程中必须考虑的。

不同频带区间呈现的曲线特征不同，在计算目标函数时就存在频带范围的选取，设置固定频带下限为0，频带上限从1Hz逐渐增加到300Hz，考察不同频带区间下目标函数的变化趋势。图2c是由图2a根据式(7)提取的目标函数值(对数显示)，发现τ值取准确值(5ms)时蓝线的值并不总是最小，换言之，5ms并不是在所有频带区间内都是最佳逼近层厚。如果频带选0～50Hz，最佳层厚是6ms(紫线)；如果频带选0～150Hz，最佳层厚是5ms(蓝线)；如果频带选0～225Hz，最佳层厚是6ms(红线)。这种情况在理论上无法预见，且在离散计算过程中是无法避免的，而图1c中的奇异点更是扩大了这种误差。图2d是由图2b根据式(10)提取的目标函数值(对数显示)，在35Hz之前，三条曲线并不能很好地区分开，但之后蓝线(对应准确值5ms)与另两条曲线明显分开且不再交叉。对比结果表明，与式(7)相比，式(10)的最佳逼近值较准确，且较为稳定，受频带区间的影响较小。

图3是顶底反射系数极性相反、比值为-4条件下的中间函数与目标函数，与图2的情况类似，这里不作赘述。图2和图3表明，目标函数的逼近效果与顶底反射系数的极性组合、比值大小无关，最佳逼近值能反映薄层的厚度特征。换言之，基于相位谱的薄层厚度估计不用考虑顶底反射系数比的影响，这是很多振幅谱方法不具备的特点。但相位方法也有假设条件，即要求地震子波为零相位、零时刻正好在两个反射系数中间、信噪比水平不能过低等。另外，在强噪情况下，目标函数容易在奇异点处取得极小值，如何实现全局最佳一致逼近也是需要进一步探索的问题。

图2 中间函数及目标函数 (a)根据式(6)计算的中间函数； (b)根据式(9)计算的中间函数； (c)根据式(7)提取的目标函数； (d)根据式(10)提取的目标函数

图3 中间函数及目标函数 (a)根据式(6)计算的中间函数； (b)根据式(9)计算的中间函数； (c)根据式(7)提取的目标函数；(d)根据式(10)提取的目标函数

3.2 稳定性分析

采用楔形模型进行目标函数的稳定性分析，地层厚度随着道号的增加而增加，双程旅行时从1ms逐渐增大到30ms，地震子波选取30Hz主频的Ricker子波，第13道的厚度大致对应λ/4。图4a是合成的楔形地震记录，图4b是提取的各道相位谱，窗口选取半径为200ms的Gauss函数，相位信息在200Hz以上开始出现凌乱的情况。

分别利用常规目标函数(式(7))和改进目标函数(式(10))估计各道最佳逼近薄层厚度，采用的相位谱频带分别为0～50Hz、0～100Hz、0～150Hz、0～200Hz、0～250Hz和0～300Hz，估算结果分别如图5a～图5f所示。图5a表明，在0～50Hz频带内，常规式无法较好地估算小于λ/8的薄层厚度；图5b～图5d表明，在0～100Hz、0～150Hz和0～200Hz三组区间上，常规式和改进式结果相近，且都与准确值接近；图5e～图5f表明，当频带区间超出200Hz时，凌乱失真的相位信息影响了薄层厚度的估算，出现较大的误差，尤其是0～300Hz区间提取的各层厚度，改进式估算的结果有三个点与准确值(蓝线)存在明显的偏离，常规式估算的结果则大部分偏离准确值，这与上节讨论的结果类似，常规式更依赖于频带区间的合理选取。

固定频带0～150Hz不变，对图4a所示的合成记录加入一定量的随机噪声，考察含噪情况下常规式与改进式的稳定性。图6是在不同噪声水平下提取的各道最佳逼近薄层厚度，图6a～图6d对应的噪声水平(最大值与无噪记录最大值的比值)分别为0.1%、0.2%、0.3%和0.5%。综合分析发现，噪声的加入对两个目标函数逼近最佳值都有很大的影响，噪声越大，影响越强；相比常规式，改进式具有相对较好的抗噪性，可能的原因是，随机噪声相位与合成记录相位叠加，有增加奇异点数并影响估算结果稳定性的可能。

图4楔形模型 (a)合成记录; (b)相位谱

图5 利用不同频带区间相位信息估算的薄层厚度结果 (a)0～50Hz； (b)0～100Hz； (c)0～150Hz； (d)0～200Hz； (e)0～250Hz； (f)0～300Hz

图6 不同信噪比水平下估算的薄层厚度结果 (a)噪声水平为0.1%； (b)噪声水平为0.2%； (c)噪声水平为0.3%； (d)噪声水平为0.5%

4 实际资料应用

选取的实际资料来自陆上工区，图7a是过工区内A1、A2、A3三口井的过井剖面，目的层为1560～1600ms(图中黑线所示)，是砂泥岩薄层，沉积相对稳定。首先，对地震记录进行插值(插值函数选三次样条)，将时间采样间隔从原先的2ms变为0.2ms。其次，选取Gauss函数作为窗口，窗口宽度参数取15～25ms(约1个波长)，沿同相轴提取目的层地震记录的相位谱，如图7b所示。用于计算目标函数(式(10))的频带分别选取0～50Hz、0～100Hz、0～150Hz和0～200Hz，四组最佳逼近结果取均值输出，作为薄层厚度的初值(如图7c蓝线所示)。显然，薄层厚度的初值存在较大的波动和奇异值(如红圈中部分)，需要对蓝线进行平滑处理。图7c中的红线是蓝线经过3点中值滤波和10点均值滤波得到的光滑处理结果(即最终输出的目的层厚度值)，随着道号的增加，薄层厚度呈先增后减的趋势。图7c中的三个绿色“o”标记，代表井上统计的目的层信息： A1井(第141道)对应的双程旅行时为2.68ms； A2井(第261道)对应的双程旅行时为6.88ms； A3井(第441道)对应的双程旅行时为1.37ms。对比红线发现，A1井和A2井处的结果与井上信息吻合度较好(分别为2.97ms和7.32ms)，

图7 实际叠加剖面处理结果 (a)叠加剖面; (b)目的层相位谱; (c)薄层厚度估算结果

表明基于相位谱能一定程度上获得目的层的厚度信息；但A3井附近的处理结果误差较大(3.18ms)，可能原因是该目的层厚度过薄，相位谱信息的提取不够准确，导致估算结果偏离准确值。

图8为根据相位属性估算所得的目的层砂体厚度分布图。分析预测结果可知，该地区砂岩薄层的厚度大致分布范围为3～10m，目的层砂体呈连片条带状分布，根据已有地质认识，工区的物源方向为东南—西北向，切片上展示的砂体厚度变化及延展方向符合已知地质规律。表1为预测结果与井上信息进行误差分析，T1、T2两口验证井的预测厚度基本在井上标准值上下浮动，相对误差在可接受范围内(除A3井外)，证明该方法在薄层预测方面应用效果具有一定的可信度和可行性。

图8 目的层砂体厚度预测分布图表1 砂体厚度预测误差表

井名井上厚度/m预测厚度/m相对误差/(%)A13.103.4410.97A26.707.075.52A32.475.73131.98T16.906.831.01T29.608.838.02

5 结论

(1)本文探讨了相位与地层厚度的关系，并尝试基于相位进行薄层厚度的估算，提出新的目标函数，改善常规目标函数存在的奇异值问题，在不同频带和信噪比水平的条件下，改进式能获得比常规式相对更稳定和可信的估算结果。

(2)在基于相位属性考察薄层厚度时，可以较好地规避顶、底反射系数组合的问题，这是振幅谱方法较难规避的难点；但相位方法也有假设条件，即要求地震子波为零相位、零时刻正好在两个反射系数中间、信噪比水平不能过低等。另外，在随机噪声的影响下，目标函数容易在奇异点处取得极小值，如何实现全局最佳一致逼近、提高抗噪能力也是需要进一步探索的问题。

(3)当地层厚度较薄时，相位谱的估计受到干扰，无法较准确地提取薄层厚度，这是本文方法存在的一个问题，也是进一步研究的重点。

(4)地震记录的相位谱是由地震子波相位谱与反射系数相位谱叠加得到的，与振幅谱信息相比，相位谱受弱反射的影响较小，其有效频带比振幅谱宽，所能利用的频率信息更丰富，但地震相位的响应机制需要更多的研究和研讨。

[1] 陈双全,李向阳.应用局部频率属性检测天然气藏.石油地球物理勘探,2012,47(5):754-757. Chen Shuangquan,Li Xiangyang.Gas reservoir detection based on local frequency attributes.OGP,2012,47(5):754-757.

[2] 张应波.地震波相位信息的分析与应用.石油物探,1988,27(4):47-63. Zhang Yingbo.Analysis and application of seismic phase information.GPP,1988,27(4):47-63.

[3] 毕俊凤.优势频带分频相位分析技术的应用.工程地球物理学报,2008,5(2):192-195. Bi Junfeng.The application of analyzing spectral decom-position phase of preponderant frequency channel.Chinese Journal of Engineering Geophysics,2008,5(2):192-195.

[4] 柏冠军,赵汝敏,杨松岭等.优势频带相位分析技术在AA区块断层解释中的应用.石油物探,2011,50(5):513-516. Bai Guanjun,Zhao Rumin,Yang Songling et al.Application of dominant frequency band phase analysis technique on the fault interpretation in AA block.GPP,2011,50(5):513-516.

[5] 王庆华.分频相位分析技术在断层解释中的应用.工程地球物理学报,2013,10(6):896-899. Wang Qinghua.The application of frequency division phase analysis technology to fault interpretation.Chinese Journal of Engineering Geophysics,2013,10(6):896-899.

[6] 姬战怀,严胜刚.分频瞬时相位识别小断层方法及应用.中国石油学会2015年物探技术研讨会,2015,502-505. Ji Zhanhuai,Yan Shenggang.The application of frequency division instantaneous phase analysis to minor fault interpretation.Chinese Petroleum Society 2015 Geophysical Exploration Technology Conference,2015,502-505.

[7] 姬战怀,严胜刚,张莲叶.用分频地震数据的瞬时相位识别小断层.SPG/SEG北京2016国际地球物理会议,2016,162-165. Ji Zhanhuai,Yan Shenggang,Zhang Lianye.Minor fault detection based on instantaneous phase analysis of frequency division seismic data.SPG/SEG Beijing 2016 International Geophysical Conference,2016,162-165.

[8] Kazmi H,Alam A,Ahmad S et al.Fault characterization using volumetric curvature from continuous phase spectrum.SEG Technical Program Expanded Abstracts,2012,31:1-5.

[9] Alam A,Taylor J D.Dip azimuth and fault from continuous phase spectrum.SEG Technical Program Expanded Abstracts,2006,25:998-1002.

[10] 刘道平.储层预测中准确的相位标定技术.石油物探,1997,36(增刊):27-32. Liu Daoping.Application of phase demarcation in re-servoir prediction.GPP,1997,36(S):27-32.

[11] 赵子豪,李凌,马跃华等.辫状河沉积储层预测技术——以大港探区孔店油田为例.石油地球物理勘探,2017,52(1):152-159. Zhao Zihao,Li Ling,Ma Yuehua et al.Braided river sedimentary reservoir prediction:an example of Kongdian Oilfield in Dagang.OGP,2017,52(1):152-159.

[12] 赵淑红,朱光明.相位在地震勘探中的用途及研究.石油物探,2005,27(4):87-89. Zhao Shuhong,Zhu Guangming.Application and research of phase in seismic prospecting.GPP,2005,27(4):87-89.

[13] Partyka G A,Gridley J,Lopez J.Interpretational aspects of spectral decomposition in reservoir characterization.The Leading Edge,1999,18(3):353-360.

[14] Partyka G A.Seismic thickness estimation:Three approac-hes,pros and cons.SEG Technical Program Expanded Abstracts,2001,20:503-506.

[15] Widess M B.How thin is a thin bed?Geophysics,1973,38(6):1176-1180.

[16] Sulistyoati T W A,Novitasari L,Winardhi S.Thickness estimation using gradient of spectral amplitude from spectral decomposition.SEG Technical Program Expanded Abstracts,2011,30:1923-1926.

[17] 孙鲁平,郑晓东,首皓等.薄层峰值频率与厚度关系研究.石油地球物理勘探,2010,45(2):254-259. Sun Luping,Zheng Xiaodong,Shou Hao et al.The studies on relationship between thin layer seismic peak frequency and its thickness.OGP,2010,45(2):254-259.

[18] 王云专,郭雪豹,邢小林等.薄层峰值频率特征分析.地球物理学进展,2013,28(5):2515-2523. Wang Yunzhuan,Guo Xuebao,Xing Xiaolin et al.Analysis of peak frequency characteristics of thin bed.Progress in Geophysics,2013,28(5):2515-2523.

[19] 柏冠军,吴汉宁,赵夕刚等.地震资料预测薄层厚度方法研究与应用.地球物理学进展,2006,21(2):554-558. Bai Guanjun,Wu Hanning,Zhao Xigang et al.Research on prediction of thin bed thickness using seismic data and its application.Progress in Geophysics,2006,21(2):554-558.

[20] 高静怀,陈文超,李幼铭等.广义S变换与薄互层地震响应分析.地球物理学报,2003,46(4):526-532. Gao Jinghuai,Chen Wenchao,Li Youming et al.Ge-neralized S transform and seismic response analysis of thin interbeds.Chinese Journal of Geophysics,2003,46(4):526-532.

[21] 范明霏,吴胜和,曲晶晶等.基于广义S变换模极大值的薄储层刻画新方法.石油地球物理勘探,2017,52(4):805-814. Fan Mingfei,Wu Shenghe,Qu Jingjing et al.A new interpretation method for thin beds with the maximum modulus based on the generalized S transform.OGP,2017,52(4):805-814.

[22] Barnes A E.Improving frequency domain thinbed analysis.SEG Technical Program Expanded Abstracts,2004,23:1929-1933.

[23] Puryear C,Castagna J P.Layer-thickness determination and stratigraphic interpretation using spectral inversion:Theory and application.Geophysics,2008,73(2):R37-R48.

[24] Nowak E J,Swan H W,Lane D.Quantitative thickness estimates from the spectral response of AVO measurements.Geophysics,2008,73(1):C1-C6.

[25] 田鑫,王绪本,张铭等.地震属性方法在油田开发阶段薄砂体识别中的应用——以印尼苏门答腊盆地Gemah油田M油层为例.天然气地球科学,2011,22(3):533-538. Tian Xin,Wang Xuben,Zhang Ming et al.Application of seismic attributes identifying thin sand body:An example of M Formation in Gemah Oilfield of Suma-trina Basin.Natural Gas Geoscience,2011,22(3):533-538.

[26] Georgy M,Viatcheslav P,Djalma M S F et al.Phase spectrum applied to pinch out zones analysis.9th International Congress of the Brazilian Geophysical Society,2005,1438-1443.

[27] 蔡涵鹏,龙浩,贺振华等.基于地震数据瞬时相位谱的地层厚度估算.天然气地球科学,2014,25(4):574-581. Cai Hanpeng,Long Hao,He Zhenhua et al.Thickness estimates from instantaneous phase spectrum of poststack data.Natural Gas Geoscience,2014,25(4):574-581.