地球物理学报  2017, Vol. 60 Issue (3): 1106-1117   PDF    
基于动态图像变形的PP与PS波层位直接匹配
蒋雪珍1,2, 芦俊3 , 王赟1,4     
1. 中国科学院地球化学研究所, 矿床地球化学国家重点实验室, 贵阳 550081;
2. 中国科学院大学, 北京 101407;
3. 中国地质大学 (北京), 海相储层演化与油气富集机理教育部重点实验室, 北京 100083;
4. 中国地质大学 (北京), 地质过程与矿产资源国家重点实验室, 北京 100083
摘要: 多分量地震资料的矢量偏移、多波地震资料的联合解释与反演均需要估算纵横波的速度比,实现纵波与转换横波在时间或深度域的匹配.基于DTW,本文实现了一种适用于PP与PS波直接匹配的动态图像变形算法.该算法分为三个部分:首先,使用二阶对称动态规划算法逐样点递归计算PP与PS波走时或深度的误差累积和;其次,在以误差累积和为目标函数的回溯阶段设定变形窗,并在纵横波速比约束的变形窗内递归回溯搜索匹配路径;最后,根据最大相关系数判定准则在匹配路径中确定最佳匹配路径,获得使PP与PS波匹配的拉伸或压缩时移量.利用所获得的拉伸压缩时移量计算纵横波速度比就可以实现PP与PS波之间的匹配.模型与实际陆上多分量地震资料测试结果表明:该方法具有较高的匹配精度,且对于信噪比、相似度较低的多分量地震资料,该方法也能产生较好的匹配效果.
关键词: 多分量地震      动态图像变形      纵横波速比      匹配      最大相关系数判定准则     
PP-and PS-waves matching directly based on dynamic image warping
JIANG Xue-Zhen1,2, LU Jun3, WANG Yun1,4     
1. Institute of Geochemistry, Chinese Academy of Sciences, State Key Laboratory of Ore Deposit Geochemistry, Guiyang 550081, China;
2. University of Chinese Academy of Sciences, Beijing 101407, China;
3. China University of Geosciences, Key Laboratory of Marine Reservoir Evolution and Hydrocarbon Accumulation Mechanism, Ministry of Education, Beijing 100083, China;
4. China University of Geosciences, State Key Laboratory of Geological Processes and Mineral Resources, Beijing 100083, China
Abstract: P-to S-wave velocity ratio estimation as well as PP-and PS-waves matching in time or depth domain are the basic procedures for the imaging, inversion, and interpretation of multi-component seismic data. Based on Dynamic Time Warping (DTW), we propose a three-step dynamic image warping algorithm which can be used to match PS-to PP-waves directly in time or depth domain of post-stack sections. With the second order symmetrical dynamic programming method, we first compute the total alignment errors of the double-way travel time or one-way depth for PP-and PS-waves. The warping window on total alignment errors is calculated with the maximum and minimum of the P-to S-wave velocity ratio to constrain the searching zone. We then search the optimal path within the warping window by using the recursive backtracking method. A maximum correlation coefficient criterion is applied to determine the best warping path. We finally calculate time-or depth-shifts between PP-and PS-images and match PS to PP post-stack section with the best shifts or the P-to S-wave velocity ratio calculated from the shifts. Our method is tested by using synthetic and field data from Xin Chang oil field. The results show that our method is effective and accurate in matching the PP-and PS-events, and could be further applied to the seismic data with low signal-to-noise ratio or weak similarity in different components.
Key words: Multi-component seismic data      Dynamic time warping      P-to S-wave velocity ratio      Matching      Maximum correlation coefficient     
1 引言

纵波与转换横波时间或深度域层位匹配是进行多波地震资料联合解释与反演的基础,层位匹配过程中估算的纵横波速度比也是多分量地震偏移成像中的关键参数 (胡晓亚和王赟,2015).以往,纵波和转换横波匹配一般采用人工拾取法以及测井标定法.Gaiser (1996)提出了利用最大相关法求取纵横波速比、平均纵横波速比以及层间纵横波速比.王赟等 (2009)推导了用多波走时参数求取速度比的系列公式.Van Dok和Gaiser (2001)利用最大相似性原理,扫描PP与PS波的速度比谱,拾取平均速度比值,完成PS与PP波在时间域的匹配.Fomel和Backus (2003)采用最小二乘和多次迭代的方法来实现PP与PS波时间域自动匹配.Yuan等 (2008)以最大相似性为迭代收敛准则,采用模拟退火算法实现PP与PS波在时间域上的匹配.Zhang和Wang (2010)采用纵横波速度比最优化估计和子波存储的校准工作流来减小匹配后的子波畸变程度.Chen等 (2014)通过时间匹配工作流来提高多分量地震数据同相轴定位精度以及联合解释反演技术.Pan等 (2016)基于线性平滑模型以及Born公式对HTI介质采用PP与PS AVAZ (azimuthal amplitude variation with offset) 数据联合反演技术来估算裂缝特性.Lu等 (2015, 2016) 利用测井合成记录的标定实现了PP与PS波的匹配,并将PS波的AVO道集压缩至PP波时间,用于后续的多波联合反演与解释.

在多分量地震数据的偏移处理过程中,为方便PP与PS波的联合解释与反演,经常需要用到层位匹配技术使得PP与PS波具有相同的走时时间刻度.Hale (2013)在DTW算法的基础上提出基于图像的动态变形技术DIW (Dynamic Image Warping, 动态图像变形) 实现PP与PS波的匹配,但在层位匹配过程中,由于PP与PS波子波和频谱特征的差异使得振幅与相位在匹配中产生畸变.Graziano和Hale (2014)提出层位匹配前的PS波子波反褶积方法来消除层位匹配过程中的子波畸变并估算PS子波.Yang等 (2014)提出基于P波速度的S波速度反演图像定位方法,该方法是在PS深度偏移成像过程中基于DIW技术将PS定位到PP图像中从而实现PP与PS波的匹配.

DIW算法是DTW在多维数据上的应用.Hale (2013)使用树序列动态规划 (tree-sequential dynamic programming, TSDP) 方法对2D剖面垂向方向上的时间序列和横向方向上的道序列分别应用DTW算法获得时移量.该方法首先通过一个常数纵横波速度比将PS波的走时压缩到PP波时间上,再采用先垂向方向上的DIW后水平方向上的DIW方法,实现PP时间尺度上的PS与PP波的匹配.Compton和Hale (2013)对此作了进一步改进,对3D PP数据进行强反射粗采样,使用平滑动态图像变形 (Smooth dynamic image warping, SDIW) 算法估算时移量,然后将估算的时移量插值平滑,计算得到纵横波速比.本文引入Hale和Compton提出的DIW、SDIW算法思想,基于二阶DTW算法,采用垂向纵横波速比约束的变形窗与水平斜率约束,并基于最大相关系数判定准则直接将PS与PP波匹配,估算纵横波速比.

2 方法原理 2.1 动态时间变形算法

DTW (Lin et al., 2010) 算法普遍应用于时间序列匹配.该算法以动态规划 (Dynamic programming, DP) 的思想为依据,最初用于解决语音识别中发音长短不一的模板匹配问题.Muñoz和Hale (2013)首次将这一算法引用到测井与地震联合处理技术领域,实现了测井数据与合成地震数据的自动匹配.

当P波从地表入射,经过阻抗差界面反射后,部分能量继续以P波形式被地表检波器接收,部分能量转换为S波被地表检波器接收.假定叠后PP与PS波记录中的某一道地震记录的振幅分别为I(tPP),I(tPS),如图 1a所示,图 1a为接收到两个阻抗差界面反射信号的PP与PS波记录图,图 1b为1a中PS波的频谱图,从图 1b中可看出,PS波的主频约为30 Hz.设对于同一反射界面的PP与PS波,二者双程走时存在一个时间间隔τ(tPP),则有

(1)

其中tPS, tPP为经过某一阻抗差界面反射到地面的双程走时,τ(tPP) 为二者的时间间隔,单位s.图 1a中,地面接收到两个阻抗差界面反射的PP与PS波的时间间隔分别为τ1(tPP) 和τ2(tPP).

图 1 叠后PP与PS记录中提取的某道数据图 (a) 左边为PP波,右边为PS波,τ1(tPP), τ2(tPP) 分别为地表接收到同一阻抗差界面反射的PP与PS波的双程走时间隔.(b) 为PS波的相位谱 (上) 与振幅谱 (下). Fig. 1 One trace of post-stack PP and PS record (a) The left is extracted from post-stack PP record, the right corresponding to the left is extracted from post-stack PS record, τ1(tPP), τ2(tPP) are differences between PP and PS double-way travel time reflected from the same interface; (b) The phase spectrum (the up) and amplitude spectrum (the down) of PS.

由于地震反射能量的差异或不同的处理过程,叠后PP与PS波存在反射振幅动态范围的差异,因此,首先需要分别对PP与PS波做整体归一化处理,设由I(tPP)、I(tPS) 归一化后的振幅分别为In(tPP)、In(tPS),则归一化后PP与PS波在时间剖面上的匹配可表示为

(2)

其中,w为加权系数,无量纲;τ(tPP)≤0,且τ(tPP) 是未知的.

PP与PS波匹配的目的是估算τ(tPP),并求取相对于t=0时刻的平均纵横波速度比,实现叠后PP与PS波的匹配.

图 1a中的PP与PS波归一化处理后,逐样点计算对齐误差,设对齐误差

(3)

计算的结果如图 2a所示,此图表示PP波中的每一采样点与PS波所有采样点的对齐误差.图 2a中PP与PS两个反射信号对应处存在4个交叉点,PP约200 ms处与PS约230 ms处反射信号对应的交叉点,以及PP约420 ms与PS约580 ms的反射对应交叉点即为所选数据的两个界面对应的反射PP与PS波双程走时.

图 2 (a) 归一化后PP与PS波的对齐误差;(b) 由公式 (5) 计算的二阶累积误差和图 Fig. 2 (a) The alignment errors for normalized PP-and PS-waves; (b) The total alignment errors calculated by equation (5)

采用二阶DTW算法估算τ(tPP).DTW算法分为三步:第一步,对公式 (3) 计算结果采用二阶对称DP算法递归计算累积误差和;第二步,利用纵横波速比计算变形窗,并以此作为约束条件;第三步,采用回溯法在变形窗内的累积误差和基于最大相关系数判定准则搜索最佳匹配路径,从而估算出最佳的τ(tPP).

2.1.1 动态规划算法

Herrera和van der Baan (2014)在DP过程中采用公式 (4) 所示的一阶对称递归求和函数,该方法在回溯搜索最优路径时求取的时移量可能会出现剧烈变化的情况.本文采用公式 (5) 所示的二阶对称递归求和函数 (Geppener et al., 2007),降低了时移量变化的剧烈程度.

(4)

(5)

其中,D(tPP, tPS) 为递归函数;δPPδPS分别为PP与PS波时间域的时间采样率,0≤tPPMδPP,0≤tPSNδPSMN分别为PP、PS波时间域的时间粗采样点数.D(tPP, tPS) 为 (tPP=tPP, tPS=tPS) 时的累积误差和,E(tPP, tPS) 为 (tPP=tPP, tPS=tPS) 时的对齐误差;公式 (5) 从 (tPP=0, tPS=0) 逐样点递归计算直到 (tPP=MδPP, tPS=NδPS),计算出所有时间采样点的误差累积和记为D(tPP, tPS),累积误差和越大,说明两道数据相似性越低,匹配程度越低.图 1中的PP与PS波记录,依据公式 (5) 计算出的D(tPP, tPS) 如图 2b所示.在PP与PS两个反射层处分别对应区域值非常小,接近于0,偏离这一区域则数值变大.

2.1.2 变形窗约束条件

根据公式 (5),最佳匹配路径搜索的时间复杂度为O (MN).由于复杂度较高,且可能搜索到无效的匹配路径,定义变形窗为误差累积和内允许匹配路径访问的元素的集合.一般采用约束条件 (即变形窗) 对匹配路径允许访问的范围进行约束.常用两种变形窗为Sakoe-Chiba Band变形窗 (Sakoe and Chiba, 1978) 和Itakura Parallelogram变形窗 (Rabiner and Juang, 1993).Sakoe-Chiba Band变形窗可表示为tPP-rtPStPP+r,其中r为一常数.图 3a为Sakoe-Chiba Band变形窗约束下的误差累积和图,变形窗为沿对角线方向的带形,变形窗外的误差累积和设置为零; Itakura Parallelogram变形窗可表示为tPP-r(tPP)≤tPStPP+r(tPP),其中r(tPP) 为tPP的函数.图 3b为Itakura Parallelogram变形窗约束下的误差累积和图,变形窗为沿对角线方向的平行四边形,变形窗外的误差累积和设置为零.在Sakoe-Chiba Band或Itakura Parallelogram变形窗的约束下,匹配路径将从 (tPP=MδPP, tPS=NδPS) 为起始点递归回溯搜索直至 (tPP=0, tPS=0).而根据纵横波速度比关系 (李庆忠,1992王赟等,2014),Sakoe-Chiba Band或Itakura Parallelogram变形窗在一定程度上提高了计算效率,但增加了大量的病态的访问元素 (Keogh and Ratanamahatana, 2005),同时去除了有效的访问元素.图 3a图 3b中,变形窗内红色区域内的可访问元素皆为病态的无效的访问元素.

图 3 (a) Sakoe-Chiba Band变形窗以及变形窗内的误差累积和图;(b) Itakura Parallelogram变形窗以及变形窗内的误差累积和图;(c) 本文中采用的变形窗图,计算变形窗所设定的纵横波速度比的范围为[1.414, 2.5].其中黑线为变形窗边界,红线为病态访问元素区域边界 Fig. 3 (a) Sakoe-Chiba Band warping window and the total alignment errors; (b) Itakura Parallelogram warping window and the total alignment errors; (c) Warping window calculate by our method, the total alignment errors extracted from Fig. 2b with P-to S wave velocity ratio limits from 1.414 to 2.5. black line stands for the boundary of warping window; red line stands for the boundary of pathological element zone

针对上述问题,根据纵横波速度比关系设定适用于多分量地震资料匹配的变形窗:

(6)

其中γmin, γmax分别为最小和最大纵横波速度比,tPPtPS分别为PP与PS波的双程走时向量.对图 1a中的PP与PS波设定最小、最大纵横波速度比分别为1.414和2.5,计算的变形窗如图 3c所示.从图 3c中可看出,PP与PS波的两个反射信号对应的时间上的误差累积和均在变形窗内,这一时间区域内的误差累积和即为匹配路径所要访问的有效元素.

当采用变形窗约束递归搜索最佳匹配路径时,由公式 (6) 联合 (1) 式得到

(7)

公式 (7) 表明了估算出的时移量的范围.超出这一范围即为无效的时移量.

2.1.3 基于最大相关系数判定准则的回溯法

DTW算法的第三步是在满足连续性、唯一性以及变形窗约束三个条件下,采用回溯法在累积误差和D(tPP, tPS) 上递归回溯搜索匹配路径,并基于最大相关系数判定准则确定最佳匹配路径.

设回溯路径为W={w0, w1, …, wk},回溯搜索匹配路径是采用DP法求取误差累积和的逆计算,因此,回溯搜索的起始位置为

(8)

tPS=(N-1) ·δ(tPS) 时,变形窗在该时刻的大小为

(9)

且有

(10)

(10) 式中.对变形窗内误差累积和以wk为匹配路径起始点遍历回溯搜索,向后回溯搜索一步wk-1=[tPP-δPP, tPS-δPS]或是两步wk-1=[tPP, tPS-δPS],wk-2=[tPP-δPP, tPS-2δPS]或wk-1=[tPP-δPP, tPS],wk-2=[tPP-2δPP, tPS-δPS],直至w0=[0, 0].

在递归回溯搜索过程中,每递归回溯一步,选择误差累积和最小的元素为该匹配路径访问的下一路径.当存在多条选择路径时,若遍历所有的有可能的路径,则耗时较长,这不仅增加了算法的时间和空间复杂度 (Goodrich and Tamassia, 2006),还提高了对计算机硬件方面的要求.由于非对角线路径使得在匹配过程中增加振幅与相位畸变的程度,因此,为解决上述问题,本文采用对角线优先法,即存在多条选择路径时,优先选择wk-1=[tPP-δPP, tPS-δPS]这一路径,同时以纵横波速度比作为路径的约束,搜索出最佳匹配路径.

在地球物理测井与地震资料联合处理过程中,相关系数法常用于判定测井-地震数据匹配质量的好坏 (Hampson-Russell, 1999);此外,相关系数法也用于地震资料处理过程中估算纵横波速比 (Gaiser, 1996; 姜镭等, 2012).

对所有起始点递归回溯搜索,可搜索到3/4S条匹配路径;依次遍历3/4S条匹配路径,计算相关系数,最终依据最大相关系数判定准则就可以确定PP与PS波的最佳匹配路径.如图 4所示,其中红线表示所确定的最佳匹配路径.利用最佳匹配路径计算PP与PS波匹配的时移量τ(tPP),依据τ(tPP) 将归一化前的PS波压缩与PP波匹配,匹配结果如图 5a所示.从图 5a中可看出,压缩后的PS波的两个反射层位处较精确地与PP波匹配.压缩后的PS波的频谱图如图 5b所示.比较压缩前与压缩后的PS频谱分析图可看出:PS波主频由压缩前的30 Hz变为40 Hz,压缩后的PS波主频与PP波一致;压缩后的振幅与相位产生畸变,这是由PP与PS波子波和频谱特征差异导致.

图 4 采用本文方法通过回溯变形窗内的误差累积和确定的最佳匹配路径图.其中,红线为最佳匹配路径 Fig. 4 The optimal path calculated by our method. Red line stands for the optimal path
图 5 (a) 左边为PP波,右边为压缩后的PS波;(b) 压缩后的PS波的相位谱 (上) 与振幅谱 (下) Fig. 5 (a) PP (left) and warped PS waves (right) warped by our method; (b) is the phase spectrum (the up) and amplitude spectrum (the down) of warped PS
2.2 动态图像变形算法

粗采样不仅能提高计算效率,还能提高估算的时移序列的精度,因此,首先对2D叠后PP与PS波粗采样.

DIW是DTW在2D数据上的应用,最简单的方法是将整体归一化且粗采样后的2D多分量地震数据看作一系列的垂向时间序列,单独对每一序列采用DTW算法估算时移量τ(tPP).但由于每道PP与PS波记录的信噪比以及在地表接收到的反射强度不一致使得二者的相似度不一致,这就导致了在以时移量压缩2D PS波时横向出现不连续的现象.为避免该现象的出现,本文在基于最大相关系数判定准则的基础上采用如下公式进行横向斜率约束.

(11)

其中dx为空间采样间隔,au表示时移量变化速率的上界.

横向斜率约束后的2D时移量记为τ(x, tPP).估算τ(x, tPP) 时,仅计算了粗采样的位置,但估算纵横波速度比并且将PS与PP波匹配需要所有采样点的时移量.因此,本文采用2D线性插值法插值,最终得到插值后的2D时移量,记为τd_i(x, tPP). 2D线性插值收敛性好,构造简单,且具有较高的数值精度,但在粗采样处时移量不光滑,且2D线性插值的准确度取决于粗采样的密集程度.因此,对插值后的τd_i(x, tPP) 采用简单平均平滑方法做平滑处理,平滑后2D时移量记为τd_i_s(x, tPP).

2.3 纵横波速度比估算

纵横波速度比是纵波与转换波资料匹配的基础 (陈双全和李向阳,2011),也是多分量地震资料偏移成像的关键 (胡晓亚和王赟,2015).采用DIW技术估算出τd_i_s(x, tPP) 后,利用τd_i_s(x, tPP) 计算相对于tPP=0时刻的平均纵横波速度比γ(x, tPP),联合公式 (7),得到平均纵横波速度比

(12)

3 模型测试

为了验证该方法的有效性,本文分别采用加入高斯白噪声前后的简单二层倾斜模型进行测试,如图 6所示.图 6a6b分别为未加噪声的叠后PP与PS剖面,图 6c6d分别为加入能量比为100:1的高斯白噪声后的叠后PP与PS剖面.经计算,二者信噪比分别为1.94 dB和-4.13 dB.该模型地层倾角小于10°;模型CDP间隔为5 m,共计50道;每道251个采样点,采样率为4 ms,记录长度为1 s.为方便显示,图中每5道显示一道地震记录.PP与PS波的地震子波为分别为40 Hz、25 Hz的雷克子波.第一个界面反射较强,第二个界面反射较弱;PP的反射强于PS波的反射.其中PP与PS波的理论纵横波速比γ(x, tPP)=1.732,理论时移量如图 7a所示.

图 6 (a) 和 (b) 分别为未加高斯白噪声时的叠后PP与PS剖面;(c) 和 (d) 分别为加入高斯白噪声后的叠后PP和PS剖面.其中,PP和PS波加入的信噪比分别为1.94 dB和-4.13 dB Fig. 6 Post-stack (a) PP-and (b) PS records without White Gaussian Noise; and post-stack (c) PP-and (d) PS records with White Gaussian Noise. The signal-to-noise ratio (SNR) of PP is 1.97dB, and the SNR of PS is -4.09 dB
图 7 (a) 理论纵横波速度比为1.732计算出的2D时移量图;(b) 和 (c) 分别为未加高斯白噪声时的2D时移量和纵横波速度比图;(d) 和 (e) 分别为加入高斯白噪声后的2D时移量和纵横波速度比图 Fig. 7 (a) Time-shift of module with the oretical P-to S-wave velocity ratio: 1.732; (b) Time-shift and (c) P-to S wave velocity ratio without White Gaussian Noise; and (d) time-shift and (e) P-to S wave velocity ratio with White Gaussian Noise

采用本文的方法,在同一条件下分别对加入高斯白噪声前后的PP与PS波匹配并估算纵横波速度比.归一化PP与PS波,设置PP与PS波时间窗为道号1-50,采样时间0~1000 ms;粗采样网格为1个CDP×4 ms,即空间和时间采样间隔都为1个采样点.由于归一化以及粗采样后的PP与PS波所有采样点振幅差较小,因此设置PS波加权系数为1.为使PS与PP波相匹配,设置速度比范围为1.414~2.5,经计算后的变形窗如图 3c所示.横向斜率约束为100%.加入高斯白噪声前后估算的τd_i_s(x, tPP) 以及平均纵横波速度比γ(x, tPP) 如图 7所示.从图 7b7d中结果可知,加入高斯白噪声前后估算的τd_i_s(x, tPP) 为0~300 ms,与γ(x, tPP) 计算出的时移量 (图 7a) 相符.从图 7c7e中可知计算出的γ(x, tPP) 约为1.73, 与γ(x, tPP) 相近.根据γ(x, tPP) 可计算出PS波的1000 ms对应PP的约为730 ms,因此,图 7中730~1000 ms之间的γ(x, tPP) 为外推时移量计算的结果.加入噪声前后未归一化的PS波分别经过图 7b7d中的τd_i_s(x, tPP) 压缩后的结果如图 8b8d所示.对比图 8a8b中的未加高斯白噪声时的PP时间尺度上的PP与PS波,可看出压缩到PP时间尺度后的PS波较好地与PP波匹配.对比图 8c8d中的加入高斯白噪声后压缩到PP时间尺度上的PS波可知,加入噪声后压缩后的PS波同样能够较好地与PP波匹配,且匹配精度与未加噪声时相当.倾斜模型经该算法计算的误差小于等于1个采样率大小.此模型结果很好地说明了本文方法对于PS匹配PP是有效的.

图 8 未加高斯白噪声时的 (a) 叠后PP和 (b) 本文方法压缩后的PS剖面图,加入高斯白噪声后的 (c) 叠后PP和 (d) 本文方法压缩后的PS剖面图 Fig. 8 (a) Post-stack PP-and (b) warped PS-record warped by our method without White Gaussian Noise; and (c) post-stack PP-and (d) warped PS-record warped by our method with White Gaussian Noise
4 实际工区数据应用

本文采用四川新场多分量地震数据进行测试.以Inline1383测线为例,部分叠后PP与PS时间剖面如图 9所示.通过对图 9中PP与PS波进行对比可知,叠后PP与PS剖面浅层成像较差,无明显的同相轴,地层无法一一对应;PP时间2500~5000 ms之间有4层较强的反射,PP与PS波能大致一一对应;PP时间5000~7000 ms之间反射较弱,同相轴不明显,PP与PS波难以对应.采用欧氏距离法计算Inline1383测线PP与PS波的相似度为5.87%,相似度较低,匹配程度差.

图 9 1383测线PP (左) 与PS (右) 剖面图 Fig. 9 Inline 1383 PP-(left) and PS sections (right) in Xin Chang

该测线CDP范围为1002-1779,CDP间距为25 m.每道数据3500个采样点,采样率为2 ms,记录长度为7 s.分别归一化PP与PS波,设置CDP范围1002-1779,采样时间0~7000 ms.粗采样设置为10CDP×6 ms,根据PP与PS波能量分析设置加权系数w=3;加权系数使得归一化后的PS与PP波能量相近.根据已有的地震测井资料可知新场这一工区的纵横波速比在1.9左右,因此,设置的纵横波速比范围为1.6~2.2;斜率约束为100%.

估算出的2D时移量τd_i_s(x, tPP) 如图 10所示;由τd_i_s(x, tPP) 计算的对于tPP=0的平均纵横波速度比γ(x, tPP) 如图 11所示.由图 11可知,由于在0~500 ms处,叠后PP波成像效果极差,且PS波无反射层位,因此,在约500 ms处,γ(x, tPP) 精确度低;500~5000 ms之间,γ(x, tPP) 约为1.8~1.9;CDP为1002-1300,时间为1000~1500 ms区域为一低速度比区,约为1.75;1500~3000 ms区域以及CDP为1340-1779,时间为1200~3000 ms区域为高速度比区域,约为2.0;3000~5000 ms区域平均速度比约为1.85.整个测线剖面γ(x, tPP) 约为1.9.未归一化的PS波经过图 10中的τd_i_s(x, tPP) 压缩后的结果如图 12右图所示.对比图 12中的PP与压缩后的PS波可看出,反射较强的PS波能够很好地与PP匹配;一些反射较弱的难以完全匹配.从整体上看,压缩后的PP时间尺度上的PS与PP剖面匹配较好.

图 10 1383测线压缩时移量图 Fig. 10 Inline1383 time-shift
图 11 1383测线平均纵横波速比图 Fig. 11 The average P-to S-wave velocity ratio of Inline 1383
图 12 1383测线PP (左) 和压缩后的PS (右) 图 Fig. 12 Inline 1383 PP-(left) and warped PS (right) sections by our method

新场Inline1383测线数据测试结果表明,对于相似度较低的PP与PS时间剖面,本文算法能够使得压缩后PP时间尺度的PS与PP波较好地匹配.同时,可估算出相对精细的γ(x, tPP),将有利于多分量地震数据的偏移速度建模.

5 结论

为使得多分量地震数据叠后PS与PP波剖面匹配,本文基于DTW,实现了一种适用于PP与PS波直接匹配的二阶DIW算法.该算法通过最大、最小纵横波速比计算变形窗来约束最佳匹配路径的搜索范围,并基于最大相关系数判定准则确定最佳匹配路径,最后通过计算时移量来估算相对t=0时刻的平均纵横波速度比.模型数据和实际数据测试得到如下的结论:

(1) 本文给出以最大相关系数为判定准则,采用变形窗纵向约束和斜率横向约束的二阶对称DIW算法.在纵横波速比或是纵横波速度未知的条件下,只需根据工区的叠后PP与PS时间剖面就可以估算出较精确的纵横波速比.

(2) 通过简单倾斜模型以及新场工区数据的验证,本文的算法对于高信噪比、相似度较高的叠后PP与PS时间剖面数据,能够得到较高精度的纵横波速度比,且压缩后的PS与PP匹配程度高;该算法在确保精确度较高的情况下不仅适用于时间域与深度域的多分量地震资料,还适用于信噪比、相似度较低的多分量地震资料的匹配.

参考文献
Chen S Q, Li X Y, Tang J M. 2014. Converted-wave time domain registration using the inverted pseudo-PS-wave attribute section. J. Geophys. Eng., 11: 015007. DOI:10.1088/1742-2132/11/1/105007
Chen S Q, Li X Y. A study on the event time registration of PP-and PS-wave based on the prestack inversion method.Changsha: Chinese Geophysical Society, 2011.
Compton S, Hale D. 2013. Estimating Vp/Vs ratios using smooth dynamic image warping.//SEG Technical Program Expanded Abstracts 2013. SEG, 1639-1643, doi:10.1190/segam2013-1300.1.
Fomel S, Backus M M. 2003. Multicomponent seismic data registration by least squares.//SEG Technical Program Expanded Abstracts 2003. SEG, 781-784, doi:10.1190/1.1818052.
Gaiser J E. 1996. Multicomponent Vp/Vs correlation analysis. Geophysics, 61(4): 1137-1149. DOI:10.1190/1.1444034
Geppener V V, Simonchik K K, Haidar A S. 2007. Design of speaker verification systems with the use of an algorithm of Dynamic Time Warping (DTW). Pattern Recognition and Image Analysis, 17(4): 470-479. DOI:10.1134/S1054661807040050
Goodrich T M, Tamassia R. 2006. Algorithm Design:Foundations, Analysis, and Internet Examples (in Chinese). Huo H W Trans. Beijing:Posts & Telecom Press, 8-12.
Graziano C, Hale D. 2014. Wavelets and warping PS seismic images.//SEG Technical Program Expanded Abstracts 2014. SEG, 1848-1852, doi:10.1190/segam2014-0362.1.
Hale D. 2013. Dynamic Warping of seismic images. Geophysics, 78(2): S105-S115. DOI:10.1190/geo2012-0327.1
Hampson-Russell. 1999. Theory of the strata program. Technical report. CGGVeritas Hampson-Russell.
Herrera R H, van der Baan M. 2014. A semiautomatic method to tie well logs to seismic data. Geophysics, 79(3): V47-V53. DOI:10.1190/GEO2013-0248.1
Hu X Y, Wang Y. 2015. The new multi-component seismic technology progress-Analysis and commentary of thesis about multi-component seismic technology on 2013 SEG Annual Meeting. Progress in Geophysics (in Chinese), 30(1): 391-400. DOI:10.6038/pg20150157
Jiang L, Ding W N, Li S. 2012. Influential factors of PP-wave and PS-wave matching. Geophysical Prospecting for Petroleum (in Chinese), 51(4): 367-370. DOI:10.3969/j.issn.1000-1441.2012.04.008
Keogh E, Ratanamahatana C A. 2005. Exact indexing of dynamic time warping. Knowledge and Information Systems, 7(3): 358-386. DOI:10.1007/s10115-004-0154-9
Li Q Z. 1992. Velocity regularities of P and S-waves in formations. OGP (in Chinese), 27(1): 1-12.
Lin J, Cervone G, Franzese P, et al. 2010. Assessment of error in air quality models using dynamic time warping.//Proceedings of the 1st ACM SIGSPATIAL International Workshop on Data Mining for Geoinformatics. San Jose, California:ACM, 38-44, doi:10.1145/1869890.1869895.
Lu J, Yang Z, Wang Y, et al. 2015. Joint PP and PS AVA seismic inversion using exact Zoeppritz equations. Geophysics, 80(5): R239-R250. DOI:10.1190/GEO2014-0490.1
Lu J, Meng X H, Wang Y, et al. 2016. Prediction of coal seam details and mining safety using multicomponent seismic data:A case history from China. Geophysics, 81(5): B149-B165. DOI:10.1190/GEO2016-0009.1
Muñoz A, Hale D. 2013. Automatically tying well logs to seismic data. CWP Report 725, 253-260.
Pan B, Sen M K, Gu H M. 2016. Joint inversion of PP and PS AVAZ data to estimate the fluid indicator in HTI medium:a case study in Western Sichuan Basin, China. J. Geophys. Eng., 13: 690-703. DOI:10.1088/1742-2132/13/5/690
Rabiner L, Juang B H. Fundamentals of speech recognition.Englewood Cliffs: Prentice Hall, 1993: 146-147.
Sakoe H, Chiba S. 1978. Dynamic programming algorithm optimization for spoken word recognition. IEEE Transactions on Acoustics, Speech, and Signal Processing, 26(1): 43-49. DOI:10.1109/TASSP.1978.1163055
Van Dok R, Gaiser J E. 2001. Stratigraphic description of the Morrow Formation using mode-converted shear waves:Interpretation tools and techniques for three land surveys. The Leading Edge, 20(9): 1042-1047. DOI:10.1190/1.1487310
Wang Y, Yang D Y, Lu J. 2009. Study of three-component converted wave velocity ratio in TI media. Journal of China Coal Society (in Chinese), 34(6): 752-755.
Wang Y, Xu X K, Yang D Y. 2014. Ultrasonic elastic characteristics of five kinds of metamorphic deformed coals under room temperature and pressure conditions. Science China:Earth Sciences, 57(9): 2208-2216. DOI:10.1007/s11430-014-4922-4
Yang D, Shang X F, Malcolm A, et al. 2014. Image registration guided wavefield tomography for shear wave velocity model building.//SEG Technical Program Expanded Abstracts. SEG, 1200-1205, doi:10.1190/segam2014-0847.1.
Yuan J J, Nathan G, Calvert A, et al. 2008. Automated C-wave registration by simulated annealing.//SEG Technical Program Expanded Abstracts. SEG, 1043-1047, doi:10.1190/1.3059105.
Zhang F, Wang Y H. 2010. Wavelet-preserved PP-and PS-wave registration. J. Geophys. Eng., 7(4): 395-403. DOI:10.1088/1742-2132/7/4/006
陈双全, 李向阳. 2011.基于叠前反演方法的纵波与转换波时间匹配研究.//中国地球物理学会第二十七届年会.长沙:中国地球物理学会.
古德里奇, 塔玛西亚. 2006.算法分析与设计.霍红卫译.北京:人民邮电出版社, 8-12.
胡晓亚, 王赟. 2015. 多分量地震技术新进展--SEG2013年会多分量地震技术论文分析与评述. 地球物理学进展, 30(1): 391–400. DOI:10.6038/pg20150157
姜镭, 丁蔚楠, 李珊. 2012. 纵、横波匹配影响因素分析. 石油物探, 51(4): 367–370. DOI:10.3969/j.issn.1000-1441.2012.04.008
李庆忠. 1992. 岩石的纵、横波速度规律. 石油地球物理勘探, 27(1): 1–12.
王赟, 杨德义, 芦俊. 2009. TI介质中三分量转换波的速度比研究. 煤炭学报, 34(6): 752–755.
王赟, 许小凯, 杨德义. 2014. 常温压条件下五种变质程度构造煤的超声弹性特征. 中国科学:地球科学, 44(11): 2431–2439. DOI:10.1007/s11430-014-4922-4