2011, Vol. 54
2. 中国科学院地质与地球物理研究所, 北京 100029;
3. 中国石油天然气股份有限公司塔里木油田分公司, 新疆 库尔勒 841000
2. Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing 100029, China;
3. Tarim Oilfield Company, PetroChina Company Limited, Korla Xinjiang 841000, China
波在传播过程中与非均匀介质的相互作用是多年来地球物理研究的理论问题之一.地质构造在地震学意义上的复杂性实际上是一个相对的概念.在地震勘探实践中，人们常常根据地震成像的效果，从地层倾角变化和速度横向变化两个方面来定性判断地下地质构造的复杂性.定量分析地质构造在地震学意义上的复杂性具有重要的理论意义和广泛的应用价值，可以准确区分不同地区和不同类型地质构造的复杂性，为地震勘探数据采集方案设计、地震数据处理流程确定和地震解释可靠性分析提供依据.特别是可以定量评价不同地震成像方法对复杂构造的成像品质，改变以往肉眼判断的方式；在地震偏移成像处理之前，准确选择与勘探靶区地质复杂性相适应的地震成像方法，提高成像效率，实现较为准确地预测复杂构造地震勘探的效果，避免地震勘探的盲目性.
用地震波探测地下地质结构时，介质的复杂性是相对于地震波长而言的.Fu[1]提出了一种利用地震波探测地质构造复杂性的定量分析方法，并通过理论模型验证了方法的有效性.首先，将复杂地质构造表示为慢度非均质谱和倾角非均质谱，以此量化速度横向变化和地层倾角变化分布.其次，通过地震成像算子的频散分析建立其角谱函数，以此描述其成像精度随慢度和传播角度的变化规律.最后，将地下复杂构造的地质非均质谱与地震成像算子的角谱进行点积运算，计算该地震成像算子对给定地区复杂构造的成像效率，其值越大，说明地震波的探测能力越强，复杂构造的地质复杂性就越小.本文将通过实际例子推广应用这种地震波探测地质构造复杂性的定量分析方法，研究我国典型地质构造单元的地震勘探复杂性，并在实际应用中进一步完善该方法的各个技术环节.
新疆库车坳陷生烃条件优越，但由于该区陡倾角地层、复杂断裂带与刺穿盐体交织发育，形成异常复杂的地下地质结构，地震资料成像困难，是我国典型的石油勘探高难度区，其复杂性在国内外少见.通过以往多年勘探，库车坳陷发现了多个有利构造，由于构造复杂和储层埋藏深，钻井周期长，勘探费用较高.该区多年的勘探实践表明，准确的地震成像是勘探成功的关键.库车坳陷高陡复杂构造呈现很强的速度横向变化和陡倾角地层分布，地震成像问题多年来一直是困扰该地区勘探开发的关键问题.本文通过定量评估库车坳陷高陡构造地震波勘探的复杂性，为该区新一轮勘探部署和地震勘探数据采集方案设计提供依据，并为准确选择与靶区地质复杂性相适应的地震成像处理方法和准确预测复杂构造地震勘探效果提供技术支撑.
2 方法原理一般地说，通过统计分析可计算地下复杂构造介质的横向速度变化非均质谱和地层倾角变化非均质谱，从地质的角度实现对复杂构造介质在横向速度变化和地层倾角变化两个方面的量化表达.然而，这种地质上的复杂构造量化表征只具有地质意义，与地震波传播探测无关.另一方面，地震波在复杂构造介质中的传播特征可通过波动方程传播算子来描述，利用地震成像算子的频散方程计算地震成像算子的角谱函数，以此定量分析地震成像精度随慢度和传播角度的变化规律.通过将复杂构造介质的横向速度变化和地层角度变化的非均质谱与地震成像算子的角谱函数相关联，可以定义一种关于地质构造相对于地震成像算子的复杂系数，实现对复杂非均匀介质地震探测复杂性的定量评估.
2.1 地下复杂构造介质速度横向变化非均质谱速度横向变化非均质谱是量化表征地下复杂构造介质地质复杂性的重要参数之一，利用统计的方法直接从地震速度模型数据计算速度横向变化的概率分布.采用波动方程偏移逐层延拓的地下介质离散方式，我们首先将一个速度模型剖分成一组叠合的水平非均质薄板，每块薄板由空间坐标定位的若干速度点构成.由于地质沉积在岩性上的延续性，板中速度的横向变化通常满足一定的分布规律，即以某一速度间隔离散化板中的速度值，则板中各速度点值的统计分布一般满足高斯概率分布；其次，用参考速度(例如板中的最小速度值)将板内各点的速度值v(x，z)换算成折射率n(x，z)来表征板内岩性的变化；最后，按照某一折射率离散间隔逐点累计各个折射率值的总数，计算其概率密度[2]，实现以概率密度作为折射率的函数(即横向速度变化非均质谱)来量化该速度板的横向速度变化分布.
速度横向变化非均质谱的具体计算步骤如下:设板中有m个折射率点值，对这些折射率点的数值分布按照大小顺序进行重新排序，用折射率离散间隔Δn对折射率数据进行区间分段，然后统计在n~n+Δn区间内的离散点总数Δm.一般来说Δm随采样间隔Δn的变化而变化，但它与板内的离散点总数之比Δm/m是一个确定的统计数，即区间内的岩性占总岩性的百分率，而Δm/(mΔn)就是单位岩性区间(指折射率在n值附近的单位区间)的分布几率.因此，板内各离散点折射率值为n的概率密度可表示为
|
(1) |
因为概率密度表示板内单位岩性区间各折射率值的分布几率，可以把它视为描述该非均质板速度横向变化分布的非均质谱函数.不同的薄板其非均质谱函数不同，与板内岩性变化密切相关.
图 1a为国际上通用的用于测试叠前深度偏移的Marmousi速度模型[3]，其复杂性主要表现为倾角不断变化的密集薄互层、浅部陡倾角断层、深部不整合面和刺穿盐丘结构，使整个模型呈现很强的速度横向变化和陡倾角分布.图 1b为从Marmousi速度模型上不同深度拾取的两块非均质薄板，两块板的速度横向变化和地层倾角变化差别较大，第一块板描述浅部陡倾角薄互层和断裂分布，第二块板展示了深部不整合面、刺穿盐丘、陡倾角地层和低速的储层.
|
图 1 Marmousi速度模型(a)及其上不同深度拾取的两块非均质薄板Slab1和Slab2(b) Fig. 1 The Marmousi velocity model (a) and two heterogeneous slabs:Slab1 and Slab 2 (b) picked at different depths from the model |
图 2为对应图 1所示两块薄板的速度横向变化非均质谱，其概率密度值的计算采用Δn=0.01的采样间隔.可见在速度横向变化非均质谱上，板1有三个主要的谱分量分布，第一个位于n=0.5附近，为嵌入薄互层中的高速地层，其谱较窄表示岩性较为单一；第二个位于n=0.7附近，概率密度较高，代表中深层大量分布的密集薄互层岩性，其谱较宽表示较强的小尺度岩性非均质分布；第三个位于n=0.95附近，概率密度较低，代表板内浅部少量的低速非均匀介质.与板1 谱分量分布明显不同，板2的三个主要的谱分量分布具有更加明显的岩性特征，第一个窄谱分量位于n=0.4附近，代表单一的盐丘岩性分布；第二个宽谱分量介于n=0.42~0.7 之间，代表中深部大量分布的密集薄互层岩性，谱曲线呈现多个波峰，表明多个特征不同的速度值相近的地层岩性分布，谱较宽，岩性非均质性较强；第三个位于n=0.85 附近，概率密度较低，代表板内嵌入的低速储层介质.对比板1和板2的谱，从速度横向变化来看，板2比板1复杂.
|
图 2 图 1所示两块非均质板的速度横向变化非均质谱 Fig. 2 Slowness-heterogeneity spectra for Slab 1 and Slab 2 showninFigure 1 |
一般而言，不同的地层薄板具有不同的谱结构，往往具有多个分布明显的主分量.地球中实际介质的非均质性可分为两类:大尺度和小尺度.大尺度非均质性确定介质平均特性的空间变化；小尺度非均质性为偏离介质平均值的扰动部分.我们可以从三个方面来解释速度横向变化非均质谱曲线:(1)根据地震成像对速度横向变化的适应条件，可以把主要谱分量沿声折射率轴的分布位置定义为:强扰动区(n=0~0.45)、中扰动区(n=0.45~0.75)和弱扰动区(n=0.75~1.0)；(2)各主要分量的“胖"与“瘦"，即谱的带宽具有明确的物理含义，谱越宽表明介质的非均质性越强；(3)各主要谱分量的幅度，即概率密度值越大表明薄板中含对应的岩性成份越多.因此，速度横向变化非均质谱上，某个谱分量越宽，代表其非均质性很强；某个谱分量幅度越大，表明板内该岩性成份所占比例大.
2.2 地下复杂构造介质地层倾角变化非均质谱非均质谱函数p(n)只量化表征地下地震波传播速度的横向变化分布.地震波探测地下复杂介质的效率既依赖于地下速度的横向变化，又与地层倾角分布密切相关.利用上述统计方法，通过简单扫描技术同样可以从速度模型数据转换的梯度场E(x，z)或地震剖面上自动拾取地下复杂构造的地质界面倾角分布.首先对速度模型数据或地震剖面进行突出弱幅值边缘特征处理，增强地质构造边界，便于地质结构的边缘检测.然后通过扫描的方法计算各点的倾角分布θ(x，z).
除了沉积地层边界外，复杂地质结构还包括大量的具有重要地质特征的不连续界面，如断层或小断裂、河道砂体边界、透镜体边界、礁体以及其他特殊岩性体的轮廓等，这些不连续性反映在图像中为边缘.由于地下地质情况具有相当的复杂性，边缘的表现特征也多种多样，对数据中包含的地质现象进行充分有效的识别、同时兼顾数据中各种潜在的边缘信息也是很有意义的.边缘检测是模式识别和计算机图像算法的重要组成部分，常用的边缘检测算法主要包括Robert算子、Prewitt算子、Laplacian算子、Kirsch 算子、Sobel算子等梯度边缘检测算子[3]，以及针对模糊图像分析的精度更高的小波多尺度边缘检测方法，进行图像的结构信息提取.这些经典的梯度算子对灰度剧烈变化的强边缘检测非常有效.由于地层界面的地震反射强弱取决于层间波阻抗差异和界面倾角，导致地震剖面上的复杂地质构造包含许多灰度变化较小的弱边缘，常规梯度算子得到的梯度值较小，在噪音的干扰下难以检测.一种用于航磁数据和重力数据处理的归一化标准偏差边缘检测算子[5](NSD)可以同时显示大振幅和小振幅梯度的边缘，能有效突出反映褶皱与断层等信息的弱边缘特征.NSD 梯度算子图像处理简单易行，对于二维图像幅度I(x，z)，其定义为:
|
(2) |
在具体实现过程中，利用移动方形窗完成对整个地震图像的边缘增强处理，合理选择移动方形窗的尺寸可以改善NSD 边缘检测算子的处理效果，达到刻画复杂构造的精细结构[6].
对速度模型数据或地震剖面进行NSD梯度算子处理，得到的梯度场E(x，z)上强弱边缘都得到了有效的突出，通过扫描进行多道信号最佳估计计算图像上各点的倾角分布，具体计算步骤如下:对图像上每一个点进行扫描，设扫描角度初值为0°，角度间隔为1°，扫描移动窗为横向Lx道和纵向Lz个样值组成，按照下式计算平均能量谱:
|
(3) |
平均能量最大值对应的扫描倾角值即为该点最佳的倾角值估计.图 3为计算得到的两块薄板的梯度场和地层倾角分布，可见板1陡倾角分布丰富，直观上比板2复杂.
|
图 3 图 1所示两块非均质板的梯度场(a)和地层倾角分布(b) Fig. 3 Gradient fields (a) and dipping-angle fields (b) for Slab 1 and Slab 2 |
我们按照某一角度离散间隔逐点累计板内各个角度值的点数，计算其概率密度，实现以概率密度作为角度值的函数(即地层倾角变化的非均质谱)来量化该非均质板的地层倾角变化分布.具体计算如下:设板内有M个角度点值，对这些角度点的数值分布按照大小顺序进行重新排序，用角度离散间隔对角度数据进行区间分段，然后统计在θ ~θ+Δθ 区间内的离散点总数ΔM.一般来说ΔM随采样间隔Δθ的变化而变化，但它与板内的离散点总数之比ΔM/M是一个确定的统计数，即区间内的角度分布占总角度分布的百分率，而ΔM/(MΔθ)就是单位角度区间(指角度在θ 值附近的单位区间)的分布几率.因此，薄板内各离散点角度值为θ 的概率密度可表示为
|
(4) |
因为该概率密度值表示板内单位角度区间各角度值的分布几率，可视为描述该非均质薄板地层倾角变化分布的非均质谱函数.
图 4为对应图 1所示两块薄板的地层倾角变化非均质谱，其角度分布概率密度值的计算采用Δθ=1°的采样间隔.可见，两块板的地层倾角变化非均质谱都有两个主要的角度分布区，分别以θ=20°(板1)和θ=17°(板2)为分界线.正如直观的判断结果，板1的35°以上的谱分量明显多于板2的谱，即从地层倾角变化来看，板1 比板2 复杂.一般来说，由于沉积环境有一定的稳定性和持续性，具有很大速度和很小速度的岩性为数较少，其百分率较低，而具有中等速度的岩性为数较多，百分率较大.同样对于地层倾角变化也存在类似的分布规律.因此，可以用高斯分布函数来拟合速度横向变化非均质谱和地层倾角变化非均质谱，从而可以获得描述这些地质非均质谱的主要特征，实现从地质上定量分析和估计地下复杂构造介质的复杂性.
|
图 4 图 1所示两块非均质板的地层倾角变化非均质谱 Fig. 4 DippingangleheterogeneityspectraforSlab1andSlab2 |
地震勘探实践表明复杂构造地震成像的效果与地下速度横向变化和地层倾角变化密切相关.我们将地下介质的地质非均质性分解为速度横向变化非均质谱和地层倾角变化非均质谱，这种地质非均质的量化表征解决了地质构造地震勘探复杂性定量估计问题的一个方面.对地震勘探而言，地质的复杂性是相对的，应该相对于地震探测的能力来定义.目前广泛应用的地震成像算子有运算高效的RayKirchhoff算子，适用于速度横向变化弱的非均匀介质；波动方程有限差分算子适用于强速度横向变化介质，但存在网格频散和算子分裂误差，计算效率也较低；半解析的波动方程Fourier变换算子适用于强速度横向变化介质，只采用快速Fourier 变换(FFT)进行波动方程偏移，运算效率高.
对某一实际地质结构，并非精度最高的成像算子就是最佳的成像算子.不同的成像算子具有不同的宽带特性、成像精度和计算效率，适用不同复杂程度的地质构造.最佳的地震成像是指成像算子的宽带特性与地下介质的地质非均质谱分布(即速度横向变化非均质谱和地层倾角变化非均质谱)相互吻合，即地质非均质谱上的大部份谱分量位于成像算子的带通范围内，以确保复杂构造成像过程中所有介质分量被有效“照明".一般来说，由于地震成像涉及海量的数据运算，为了获得高效的运算速度，对于简单介质可采用运算高效的低阶地震成像算子.为了量化表征地震成像算子的这种“照明"特征，我们需要对成像算子进行频散分析，将频散方程表示为角谱函数.一般地说，波动方程地震成像算子的成像精度可以通过其频散方程显式地表示为速度横向变化和传播角变化的函数.
本文采用性能优越的波动方程Fourier变换算子进行地质构造地震勘探复杂性定量分析.这类成像算子可以根据速度横向变化和地震波传播角度进行逐级构造，得到不同精度的成像算子，具有不同计算效率，适应不同的地质复杂性.目前，Fourier变换地震成像算子的研究发展迅速，主要包括相移(PS)、分裂步(SSF)、退化(DP)和Fourier有限差分(FFD)等方法.这些Fourier变换成像方法具有相同的算法结构和简单的频散关系，但计算效率和精度不同，其波场延拓成像过程可统一表示为
|
(5) |
式中kx2+kz2=k02，k0 为背景波数k0 =ω/v0，v0为此薄板中的背景地震速度，ω 为角频率.媒介波场$\hat{u}$(kx，z)可根据不同的Fourier变换成像算子取不同的表达方式.例如，对于零阶Fourier变换成像的相移法，有$\hat{u}$(kx，z)=u(kx，z)，适用于速度横向不变的层状介质.
对于一阶Fourier变换成像的分裂步法(SSF)，媒介波场可表示为
|
(6) |
式中FTx为从x￫kx的正向Fourier变换.该方法只适用于弱的速度横向变化或较小的地层倾角，其频散关系由下式给出
|
(7) |
式中kx=kx/k0，kz=kz/k0.对于二阶Fourier变换成像的退化算子(DP)方法[7-9]，其媒介波场可表示为
|
(8) |
式中系数C(kx)是波数kx的函数，但与折射率n无关.该方法频散关系为一种退化的算子表达:
|
(9) |
式中a1 和b1 为常数.由方程(8)可知，退化算子方法的计算过程与传统分裂步法成像方法相似，它通过在两个分裂步之间作波数域线性插值来实现波场延拓，每延拓一层用三次FFT，比分裂步法多一次快速Fourier变换，但适用于较强速度横向变化和陡倾角地层.
对于强对比介质，还可以采用有限差分(FD)与Fourier变换混合的、精度更高的方法(FFD)[10， 11]，这些混合的方法理论上容许更大传播角度的波场延拓和比纯有限差分法更大的波场延拓步长，其频散关系为有理逼近的算子表达:
|
(10) |
式中，系数aj(n)和bj(n)是折射率的函数，随横向速度变化而变化.
根据上述频散关系可以确定Fourier变换成像算子的宽带特性，从而量化其“照明"能力.我们利用相对相位差e=δΦ-1(δΦ 为波传播的相位扰动)来表示波动方程地震成像的精度.根据频散方程(7)、(9)和(10)，传播角θ、折射率n和地震成像精度e三者的关系可以解析地表示为一种角谱函数，然后通过传播角和折射率的交汇图显示出来.例如，分裂步法地震成像算子的角谱函数可以表示为
|
(11) |
一个地震成像算子的角谱函数有两种表达方式.一种是θ=f(n)，即把传播角θ作为折射率n的函数进行绘图；另一种是δn=g(θ)，即把折射率的变化δn=1-n作为传播角θ 的函数进行绘图.
图 5 比较了三种Fourier 变换地震成像算子(SSF、DP 和FFD)在相对相位差e=10%精度下计算的角谱f(n)和g(θ)曲线.可见，DP和FFD 两种成像算子无论是对折射率或传播角都具有较大的宽带特性，成像精度明显优于SSF 算子.在计算量方面，SSF算子每延拓一层用两次FFT，DP 算子每延拓一层用三次FFT.由于这两种算子只用FFT 进行波场延拓运算，可以充分利用目前发展迅速的向量化快速Fourier变化技术来大幅度提高三维叠前偏移运算效率.FFD 算子每延拓一层用两次FFT和一次FD 校正，其中FD 校正制约了该方法在三维叠前偏移的运算效率.总之，在给定成像精度条件下，最佳的地震成像要求成像算子的角谱f(n)和g(θ)分别“照明"速度横向变化非均质谱上和地层倾角非均质谱上的所有主要谱分量.
|
图 5 三种地震成像算子(SSF、DP和FFD)在相位相对误差e=10%精度下的角谱f(n)和g(θ)曲线 Fig. 5 Angular spectra f(n)and g(θ) of the SSF，DP，and FFD migrators under the relative phase error e=10%: |
复杂构造的地震成像效果实质上取决于地震成像算子的角谱(f(n)和g(θ))与地下介质地质非均质谱(p(n)和q(θ))之间相消或相长的相互作用.大部分地震成像算子由于其平方根方程的数值逼近而不能全局地兼顾所有的速度横向变化谱分量和地层倾角变化谱分量.当波场延拓穿过一块非均质板时，成像算子与地质非均质两种相互独立的谱分量之间的相干作用是指在速度非均质谱和地层倾角非均质谱中，那些位于成像算子角谱通放带内的地质非均质谱分量将得到有效成像；相反那些位于成像算子角谱通放带外的地质非均质谱分量被削弱，成像效果变差并形成偏移噪音.
具体而言就是对于给定的成像精度，有效的地震成像要求成像算子的角谱f(n)“照明"速度横向变化非均质谱p(n)上大部分的谱分量；同时成像算子的角谱g(θ)通放地层倾角非均质谱q(θ)上大部分的谱分量.上述速度横向变化和地层倾角这两种介质特性在成像过程中是相互偶合在一起的，为了量化表征这一相互偶合作用过程的一致性，我们定义如下的地震成像效率:
|
(12) |
式中ηn和ηθ 分别为成像算子对速度横向变化和地层角度变化的成像效率.因此，与地震成像效率相对应的地下介质变化的地震探测复杂系数可表示为φ =1-η，实现对地下复杂介质地震探测复杂性的定量评估.
表 1列出了根据三种地震成像算子(SSF、DP和FFD)计算得到的图 1所示两块非均质板的地震成像效率及其地震探测复杂系数.可见，板1的速度横向变化成像效率略高于板2，但其地层倾角变化成像效率略低于板2.这与两个板的地质非均质谱上主分量的分布相吻合.从图 2 的速度横向变化非均质谱上可见，相对于板1，板2的速度横向变化非均质谱上几个主要的谱分量分布在更强的速度扰动区，并且带宽也更大，即在速度横向变化方面，板2比板1更复杂.从图 4的地层倾角变化非均质谱上可见，相对于板2，板1的谱带宽更大，其35°以上的谱分量明显多于板2 的谱，即从地层倾角变化来看，板1比板2更复杂.因此，从综合考虑速度横向变化和地层倾角变化得到的地震探测复杂系数来看，两块非均质板的地震探测综合复杂性比较接近，板2 略高于板1.表中量化表征地质复杂程度的这些数字与直观的定性判断结果相一致，整个地质模型的复杂性是其所有非均质板地震探测复杂系数的平均或累加.
|
表 1 由三种地震成像算子计算得到的图 1所示两块非均质板的地震成像效率和地震探测复杂系数 Table 1 Imaging efficiencies and complexity coefficients of Slab 1 and Slab 2 by SSF，DP and FFD migrators |
总之，随着速度对比和地层倾角的加大，非均质板变得越来越复杂，其复杂系数取决于所采用的地震成像算子，若采用SSF 地震成像算子，两块板的复杂系数都很高，但采用DP 和FFD 地震成像算子，两块板的复杂系数降低很多.利用复杂系数可识别速度模型中严重破坏地震成像品质的异常复杂区，然后利用更高精度的成像方法对该异常复杂的局部区域进行成像可提高地震成像的效率.根据选择的地震成像算子来定量分析一个实际速度模型或地震剖面地质结构的复杂性是目前地震成像研究的主要问题之一.研究地震成像算子与介质非均质谱相互作用的尺度化特征是定量评估地震成像品质的关键所在，特别是对成像结构的精细程度与地质构造复杂性进行相关性分析.
3 库车坳陷高陡构造地震波勘探复杂性分析本文应用上述介绍的方法对库车坳陷高陡构造地震勘探的复杂性进行定量分析.首先分析库车坳陷高陡构造的地质复杂性，分别计算其横向速度变化非均质谱和地层倾角变化非均质谱，然后通过地震成像算子角谱与这些横向速度变化非均质谱和地层倾角变化非均质谱的相互作用，分别应用SSF、DP 和FFD 三种地震成像算子研究库车坳陷高陡构造地质复杂性的地震可勘探能力.
3.1 库车坳陷复杂高陡构造分布库车坳陷是一个叠加复合型坳陷，发育于天山褶皱带和塔里木板块的结合部，其北紧靠南天山晚古生代陆缘盆地，南邻塔北隆起，东接库鲁克塔格，西与阿瓦提凹陷相邻，为一东西向条带状山前坳陷.山前发育了北西- 南东走向的洪积扇，受水系切割严重，冲沟发育.地表相对高差较大，最高海拔达3500m, 最低约1280 m, 相对高差达到2300 m.库车坳陷是一个以中、新生界沉积为主的前陆盆地，地层结构从老到新具有两种不同的特征，其下为稳定的震旦系结晶变质岩基底，其上披覆着巨厚的中、新生界陆相碎屑岩，地表大部分为第四系，山体局部出露第四系地层.坳陷结构特征为一强烈变形的山前逆冲带，发育一系列不完整的逆冲推覆构造体系，且成排成带分布形成了凸起和凹陷相间的展布特征.库车地区地震勘探具有两重复杂性，首先是该区复杂的地形起伏和近地表岩性横向变化，地震资料的质量和室内处理效果在很大程度上取决于地形起伏、地面地质风化、地表岩性变化和表层地质结构的复杂程度.复杂近地表地质条件对地震勘探而言，表现为近地表强散射噪音环境、近地表大尺度的速度横向变化、近地表小尺度的速度横向变化和近地表吸收衰减环境.这些复杂因素的地震综合反映是复杂畸变的地震信号走时、强烈的信号吸收衰减和极低的资料信噪比.其次是逆冲推覆陡倾角构造，刺穿盐体大面积分布，断裂带破碎等.图 6 为克拉2解释构造模型(图 6a)和过西秋克拉2 构造带的地震偏移剖面，展示了这种陡倾角地层、复杂断裂带与刺穿盐体交织发育的逆冲推覆高陡构造的复杂性，形成强速度横向变化和陡倾角地层分布的地震地质条件，导致严重的地震成像问题，陡倾角地层和大角度断裂带地震成像出不来、模糊、成像错位等.
|
图 6 克拉2解释构造模型(a)和过西秋克拉2构造带(b)的地震偏移剖面 Fig. 6 KL2 geological structure model (a) and seismic migration section (b) crossing the Xiqiu-KL2 structure zone |
针对库车地区复杂的地下地质情况，本文对穿越库车坳陷四个重点构造(大北、却勒、博孜、西秋4和西秋10)的二维地震测线进行了地震波勘探复杂性定量分析.这四个重点构造的具体位置和二维地震测线(虚线)部署如图 7所示.在这些构造上，由于横向速度变化极大，地层倾角陡，存在大角度复杂断裂带和大量的小断块，致使地震资料成像困难.这些地区的地震资料攻关处理结果普遍存在偏移归位不准，大量同相轴收敛不到位，造成控制不同构造带的断面不清，中深部构造形态出不来.另外，叠前叠后的去噪处理导致平滑作用过重，波形特征不明显，断点不清楚，很多小断层被抹掉了.总体上，地震剖面的品质较差，信噪比低，局部层次不清，逆掩推覆体下盘反射很乱，造成地震追踪解释困难.多年来的地震勘探进展缓慢，有必要重新认识和评估库车坳陷地震波勘探的复杂性.
|
图 7 大北、却勒、博孜、西秋4和西秋10四个重点构造的位置和二维地震测线(虚线)部署 Fig. 7 Location of four key structures and 2D seismic survey lines (dashed line) |
库车坳陷高陡构造多为盐相关复合构造，形成强速度横向变化和陡倾角地层分布.盐层是控制盐上、盐下变形的主要因素之一，库车含盐地区构造建模必须综合考虑盐上、盐层和盐下三个构造层的变形耦合.盐上构造层以断层相关褶皱理论为指导，盐下构造层以造山带冲断楔变形理论为指导，以物理模拟实验为手段，综合应用地表露头、地层倾角、钻井和地震等资料，进行构造建模工作.盐上构造层的地层分布特点、厚度变化和变形形态以及盐层的分布都可以为盐下构造层的变形提供制约条件，合理的盐下层变形的物理模拟和数值模拟结果必需同时符合盐上构造层的现今构造形态.
本区以往的地震勘探资料以及钻井资料表明，本区中、新生代地层发育较全，且断裂发育，构造主体部位破碎严重.在构造南、北两翼的山前及戈壁砾石区，地震资料品质较高，T8以上反射层可以连续追踪，反射能量较强；但在山体区，尤其是构造主体部位，内部反射杂乱、地震成像困难，信噪比和分辨率较低，各地层和断层的接触关系和边界不清晰.通过近几年来对库车地区速度场详细的研究，盐相关复合构造存在的问题主要表现为:(1)膏盐层厚度分布不均，速度横向上变化异常，影响盐下构造形态；(2)构造主体部位浅层普遍复杂，速度规律难以把握；(3)深层构造相互叠置，给速度建场带来一定困难；(4)逆掩断层上下盘存在较大速度差，造成隐伏构造在时间域的上拉现象.图 8 为大北和博孜构造四条地震测线速度模型.沿库车坳陷大北至博孜走向，盐下构造发育，盐层的横向岩性组成复杂，盐层厚度变化较大，导致盐层速度横向变化大，在时间域剖面上易形成速度陷井，盐下易出现假构造.特别是在大北地区，盐拱构造普遍发育，形成异常复杂的盐侧翼和盐下构造成像问题.根据实验地质学研究结果，挤压变形过程中盐层的聚集主要受盐下层古地貌形态的影响，往往出现在地形呈坡折或台阶状的位置.这些坡折或台阶状主要是由于造山带向陆内推覆过程中的挠曲变形或盐下构造层早期断裂活动产生.
|
图 8 大北和博孜构造4条地震测线速度模型 Fig. 8 Velocity models for four seismic survey lines across the Dabei and Bozi structures |
从图上可知，盐下小断裂异常发育，地层起伏也较大.由于盐层内速度横向变化和盐层厚度分布不均匀性，对盐下成像、时深转换和构造落实影响较大.图 9为大北构造速度横向变化非均质谱和地层倾角变化非均质谱，由3 条穿越大北构造的二维地震测线速度模型计算统计得到.图 10为博孜构造速度横向变化非均质谱和地层倾角变化非均质谱，由4条穿越博孜构造的二维地震测线速度模型计算统计得到.可见，大北至博孜一线构造带的速度横向变化较大，速度扰动主要分布在0.4~0.6 之间，属于强横向速度变化范畴，由2~3 个峰值组成，表明岩性的变化较为复杂.地层倾角分布第一个主分量在0°~20°之间，表明大部分地层为小角分布；第二个主分量在80°左右，属于盐拱导致的陡倾角分布，其它的倾角分量有介于20°~80°之间，主要代表盐下部分陡倾角地层和小断裂的分布，是地震成像的畸变分量.表 2列出由三种地震成像算子计算得到的大北构造和博孜构造地震成像效率和地震探测复杂系数，从复杂系数来看，博孜构造比大北构造稍微复杂一些，主要表现在，博孜构造的陡倾角分量比大北构造多，二者的横向速度变化相当.总之，就时间偏移而言，盐层厚度越大，成像畸变越明显，导致盐下构造的倾伏程度也越大；盐层与周围地层之间的速度差别越大，构造高点的偏移越明显；定井位时必需避开盐层最厚的地方，这些地方的高点一般都是假高点，工程条件复杂，失利风险很大，一般将井点选在两侧盐厚度较薄的地方.
|
图 9 大北构造速度横向变化非均质谱(a)和地层倾角变化非均质谱(b) Fig. 9 Slowness (a) and dipping-angle (b) heterogeneity spectra for the Dabei structure |
|
图 10 博孜构造速度横向变化非均质谱(a)和地层倾角变化非均质谱(b) Fig. 10 Slowness (a) and dipping-angle (b) heterogeneity spectra for the Bozi structure |
|
|
表 2 由三种地震成像算子计算得到的大北构造和博孜构造地震成像效率和地震探测复杂系数 Table 2 Imaging efficiencies and complexity coefficients of the Dabei and Bozi structures by SSF，DP and FFD migrators |
图 11为却勒-西秋10-西秋4 构造带6 条地震测线速度模型.沿库车坳陷却勒、西秋10 至西秋4构造带走向，多为盐相关逆冲推覆复合构造.却勒地区由于断层滑脱面发育数量不同，可分为两种构造样式.第一种是以膏盐体为界的盐上和盐下构造层，盐上构造层为由北倾沿古近系膏盐岩顶面滑脱的断层控制的断层传播褶皱，背斜的高点位于盐丘的顶部.时间域剖面上，盐下构造层为隐伏于盐丘边缘的低幅度背斜构造.第二种是发育两个滑脱面，即北部大型逆冲推覆底面滑脱面和盐丘顶部滑脱面.以两个滑脱面为界可分为三个构造层，上构造层为由北向南远距离推覆形成的推覆构造；中构造层是两个滑脱面之间的构造层，由于受南北挤压应力和上构造层的压实作用的影响，岩性非常致密，形成一套高速地层；下构造层为盐下受古生界古隆起控制的潜山披覆构造，时间域构造高点位于盐丘北部盐层厚度消失部位.同样，西秋4 至西秋10 构造带多为两个滑脱面为界的三个构造层模式，形成巨厚膏盐岩塑性变形下两背冲断层控制的褶皱构造.
|
图 11 却勒-西秋10-西秋4构造带6条地震测线速度模型 Fig. 11 Velocity models for six seismic survey lines across the Quele-Xiqiu10-Xiqiu 4 structure zone |
盐层横向分布及其速度变化对构造落实影响巨大.盐上普遍发育盐(断层)相关褶皱、盐撤坳陷、盐推覆、塌陷、倒转背斜，资料品质差、波场复杂、速度场建立难，同一套地层的层速度在不同的地区、不同的构造部位均有较大变化；盐层中形成盐墙、盐脊、鱼尾构造，盐刺穿，膏盐岩段表现为杂乱反射；盐下为基底古隆起低幅度背斜构造，小断裂较为发育，盐下地层层速度在很大程度上受膏盐岩厚度变化的影响.总之，由于受构造运动影响，盐上构造层地层比较破碎，地震成像效果很差；盐下构造层受巨厚膏盐岩对地震波的屏蔽的影响，构造部位的成像畸变严重.
图 12至图 14为分别由3 条、4 条和3 条穿越相应构造的二维地震测线速度模型计算统计得到的却勒、西秋10和西秋4构造速度横向变化非均质谱和地层倾角变化非均质谱.可见，沿却勒至西秋一线构造带的速度横向变化大.速度横向变化非均质谱带宽较大，主要分布在0.4~0.7 之间，属于强横向速度变化范畴，特别是0.4 左右的强速度扰动主分量所占比例较大，引起严重的地震成像问题(成像畸变和偏移噪音).谱由多个峰值组成，0.7~1.0之间的岩性分量也相对较强，表明岩性的变化较为复杂.地层倾角分布第一个主分量在0°~20°之间，表明大部分地层为小角分布；第二个主分量在80°左右，相对于大北-博孜构造带，该陡倾角分量所占比重较大.其他的倾角分量介于20°~80°之间，明显比大北-博孜构造带强，而且沿却勒、西秋10 至西秋4构造走向，大于20°的角度分量明显增强，地震成像难度越来越大.总之，西秋4 构造最为复杂，其0.4左右强速度扰动分量所占比重最大，大角度分量也最多，地震成像的难度最大，形成严重的成像畸变和很强的偏移噪音.
|
图 12 却勒构造速度横向变化非均质谱(a)和地层倾角变化非均质谱(b) Fig. 12 Slowness (a) and dipping-angle (b) heterogeneity spectra for the Quele structure |
|
图 13 西秋10构造速度横向变化非均质谱(a)和地层倾角变化非均质谱(b) Fig. 13 Slowness (a) and dipping-angle (b) heterogeneity spectra for the Xiqiu 10 structure |
|
图 14 西秋4构造速度横向变化非均质谱(a)和地层倾角变化非均质谱(b) Fig. 14 Slowness (a) and dipping-angle (b) heterogeneity spectra for the Xiqiu4 structure |
表 3列出由三种地震成像算子计算得到的却勒、西秋10和西秋4构造地震成像效率和地震探测复杂系数，从复杂系数来看，沿却勒、西秋10至西秋4构造走向，复杂系数逐渐增加，这与前面从非均质谱图上的定性分析相吻合.最复杂的西秋4 构造速度横向变化成像效率和地层倾角变化成像效率均比较低.与表 2比较可见，却勒至西秋构造带明显比大北至博孜构造带复杂，勘探的难度也更大，主要表现在盐层厚度较大和盐上断裂褶皱带异常发育，形成陡倾角和强横向速度变化的盐相关逆冲推覆复杂构造，这种复合结构是地震成像长期攻关的目标.
|
|
表 3 由三种地震成像算子计算得到的却勒、西秋10和西秋4构造地震成像效率和地震探测复杂系数 Table 3 ImagingefficienciesandcomplexitycoefficientsoftheQuele，Xiqiu10andXiqiu4structuresbySSF，DPandFFDmigrators |
一般来说，复杂构造地震成像的品质取决于成像算子与地下复杂构造地质非均质谱的相干作用结果.最佳的成像效果往往是成像算子角谱的宽带特性与复杂构造介质的地质非均质谱(横向速度变化和地层角度变化)分布相一致，即地质非均质谱上大部分谱分量介于成像算子角谱的通带范围内.本文中，我们首先利用统计的方法从地质的角度定量表征地下复杂地质构造的非均质变化，即复杂构造介质的横向速度变化非均质谱p(n)和地层角度变化非均质谱q(θ).这些表征复杂构造复杂性的谱函数只有地质意义，与地震探测方法和技术无关.为了量化表征地震成像算子的探测能力，我们从成像算子的频散方程出发，通过构造角谱函数将其成像精度表示为地下横向速度变化和地震波传播角度变化的函数.一个地震成像算子的角谱函数有两种表达方式，一种是θ=f(n)，把传播角θ作为折射率n的函数进行绘图；另一种是δn=g(θ)，即把折射率变化δn=1-n作为传播角θ 的函数进行绘图.
在给定成像精度条件下，最佳的地震成像要求成像算子的角谱f(n)和g(θ)分别“照明"速度横向变化非均质谱上和地层倾角非均质谱上的所有主要谱分量.将地震成像算子的角谱f(n)和g(θ)分别与复杂构造介质的横向速度变化非均质谱p(n)和地层角度变化非均质谱q(θ)作点积来定义该成像算子对给定地区复杂构造介质的成像效率η，从而实现定量表征地震成像中成像算子与地下复杂构造地质非均质谱的相干作用过程.成像效率η 越大，说明地震波的探测能力越强，复杂构造的地质复杂性就越小.因此，与地震成像效率相对应的地下介质变化的复杂系数可定义为φ =1-η，从而实现对地下复杂介质地震探测复杂性的定量评估.
我们以塔里木盆地库车坳陷地区具有代表性的速度模型为例，验证上述方法的可行性.库车坳陷高陡构造多为盐相关逆冲推覆复合构造，盐层是控制盐上、盐下变形的主要因素之一.盐层厚度分布横向变化较大，盐内速度横向变化异常，并且与周围地层速度差异大，严重影响盐下构造形态的成像；盐上普遍发育多条逆冲断层作用下的断裂褶皱构造破碎带，带内地震射线路径十分复杂，对地震波能量的下传产生严重的屏蔽作用，反射信号可追踪性差，成像极其困难；盐下低幅度隐伏构造地层上断裂发育，时间域成像畸变严重，由于盐上和盐中的影响，盐下地震资料的品质比较差.利用穿越大北构造、博孜构造、却勒构造、西秋10构造和西秋4构造的共17条二维地震测线速度模型，本文计算了这些构造的速度横向变化非均质谱和地层倾角变化非均质谱，由此得到这些构造的地震成像效率和地震探测复杂系数，实现定量分析这些高陡构造地震勘探复杂性.结果表明，库车坳陷速度横向变化非均质谱由多个峰值组成，主要分量集中分布在0.4~0.7的强速度扰动区，表明岩性的变化较为复杂；地层倾角变化非均质谱上第一个主分量在0°~20°之间，表明大部分地层为小角分布；第二个主分量在80°左右，主要是盐拱和断裂引起的陡倾角分量.复杂系数计算结果表明，沿却勒至西秋一线构造带明显比大北至博孜构造带复杂，复杂度沿却勒、西秋10 至西秋4 构造走向逐渐增加，主要表现为0.4 左右的强速度扰动主分量所占比重增大，介于0.7~1.0之间的岩性分量和20°~80°之间的倾角分量也相对较强，地震成像的难度逐渐加大，时间域成像易于形成严重的成像畸变和很强的偏移噪音，建议加强速度建模和叠前深度成像攻关研究.
致谢感谢中国石油天然气股份有限公司塔里木油田分公司提供数据和资料.
| [1] | Fu L Y. Quantitative assessment of the complexity of geological structures in terms of seismic propagators. Sci China Earth Sci , 2010, 53: 54-63. DOI:10.1007/s11430-009-0167-z |
| [2] | Wu R S, Xu Z, Li X P. Heterogeneity spectrum and scale-anisotropy in the upper crust revealed by the German Continental Deep-Drilling (KTB) Holes. Geophys| Res| Lett| , 1994, 21: 911-914. |
| [3] | Bourgeois A, Bourget M, Lailly P, et al. Marmousi model and data: Proc. 1900 EAEG workshop on Practical Aspects of Seismic Data Inversion , 1991. |
| [4] | 徐建华. 图像处理与分析. 北京: 科学出版社, 1992 . Xu J H. Image Processing and Analysis (in Chinese). Beijing: Science Press, 1992 . |
| [5] | Cooper G R, Cowan D R. Edge enhancement of potential-field data using normalized statistics. Geophysics , 2008, 73(3): 1-4. |
| [6] | 陈学华, 贺振华, 黄德济. 基于归一化标准偏差的地震资料边缘检测. 北京: 石油工业出版社, 2008 : 245 -249. Chen X H, He Z H, Huang D J. Edge detection in seismic data using normalized statistics. In:Exploration Geophysics Corpus (in Chinese). Beijing: Petroleum Industry Press, 2008 : 245 -249. |
| [7] | Fu LY, Duan Y. Fourier depth migration methods with application to salt-related complex geological structures. 72th Annual International Meeting, SEG, Expanded Abstracts , 2002: 895-898. |
| [8] | Fu L Y. Comparison of different one-way propagators for wave forward propagation in heterogeneous crustal wave guides. Bull| Seismol| Soc| Am| , 2006, 96: 1091-1113. |
| [9] | 符力耘, 孙伟家, 李东平. 退化的Fourier偏移算子及其在复杂断块成像中的应用. 地球物理学报 , 2007, 50: 483–495. Fu L Y, Sun W J, Li D P. Degenerate Fourier migrators for imaging complex fault zones. Chinese J| Geophys|(in Chinese) (in Chinese) , 2007, 50: 483-495. |
| [10] | Collins M D. A split-step Padé solution for the parabolic equation method. J| Acoust| Soc| Am| , 1993, 93: 1736-1742. |
| [11] | Ristow D, Rühl T. Fourier finite-difference migration. Geophysics , 1994, 59: 1882-1893. DOI:10.1190/1.1443575 |