基于复杂构造区的采集观测系统设计与优化方法

何宝庆，吕盼盼，何永清，宁宏晓，马 兰，卢秀丽，蔡锡伟.

(东方地球物理公司采集技术中心，河北涿州 072751)

近年来，“两宽一高”技术在地震勘探中发展迅猛，宽频、宽方位、高密度等技术相结合的方式逐渐代替传统采集方法，使得地震资料品质有了明显的提高，采集成本也成倍增加。采集观测系统的成本和对地下介质的探测能力受很多参数影响，例如炮距、道距、最大偏移距、纵横比、覆盖次数、面元尺寸和覆盖密度等[1-2]。长期以来，采集参数设计大部分沿用20年前基于水平层状介质的假设，在面元和道距设计方面，仍然从单个论证点出发，根据地层倾角、绕射收敛等信息进行采集方案设计；而在最大偏移距的设计方面，仍然根据速度分析精度、动校拉伸等方面的要求，这些技术都是基于水平叠加时代地震资料处理技术相对落后的产物。叠前偏移技术的广泛应用，意味着在地震数据叠前偏移时不需要做动校正，这样就不必考虑动校拉伸等因素的限制。随着勘探技术的发展和计算机水平的提高，我们应该发展更能满足当前地震勘探要求的采集设计方法。

复杂构造地区，特别是逆掩推覆区，以及含盐丘构造地区的观测系统设计与优化[3]，是近年来很多专家学者重点研究的问题。一方面，高速构造体和盐丘体阻碍了地震波的穿透，形成照明阴影区；另一方面，由于地下构造的复杂造成了地震波场的复杂，常规的速度分析手段很难适应。这势必导致使用同一观测系统采集时，工区某些部位成像较好，但其他某些部位成像则较差。对此我们需要针对不同的构造部位采用不同的手段进行观测系统设计。基于复杂模型的波动方程照明分析和叠前偏移分析是目前常用的观测系统评价方法。

波动方程照明分析方面，谢小碧[4]等从照明矩阵出发，研究了入射波和反射波在角度域与反射面的相互作用，并讨论了照明分析在观测系统设计中一些可能的应用。董良国[5]等提出利用波动照明确定针对目的层的地面最优炮点分布，并利用射线追踪和波动方程模拟，分析目的层的CRP点覆盖次数等，得到针对目的层的最优检波器排列方式和排列长度。赵虎[6]等基于地震波照明结果采用局部加密炮点规则布置接收排列的方式解决目的层反射能量不均的问题。许银坡[7]等提出了提高目的层阴影区成像质量的观测系统优化设计方法，达到使用最少的有效激发点提高阴影区照明强度，改善地震剖面质量的目的。这些方法在一定程度上提高了目的层阴影区的照明能量，但是仍然无法判断出哪种观测系统能够满足我们的勘探要求。

叠前偏移成像方面，周旭[8]等利用复杂构造区三维物理模拟数据和实际采集的宽方位三维地震数据，进行了观测系统参数的退化处理试验，分析了道密度、观测宽度、线距变化对叠前偏移的影响。姚江、秦广胜[9-10]等通过分析观测系统参数对叠前偏移成像效果的影响，利用双聚焦、叠前偏移响应、采集脚印分析等技术合理选择观测系统参数。但是当偏移速度模型准确时，叠前深度偏移成像对最大偏移距、道距、炮密度等参数的敏感性较差，所以利用叠前偏移成像来分析观测系统不是一个有效的手段。

本文在全波形反演方法研究的基础上，利用全波形反演的高分辨率特性，提出了将全波形反演应用于复杂高陡构造区采集观测系统设计的观点。试验表明，只有达到全波形反演能够收敛的最小的“最大偏移距”才能保证全波形反演的准确性。

1 方法原理

自Tarantola 1984年提出反演方法以来[11-13]，很多学者投入全波形反演方法的研究中，并不断地将其完善。目前，对于模型数据而言，已经可以得到满意的结果；对于采集得到的地震数据，通过逐步进行全波形反演，也可以得到相对精确的速度。

在地质模型存在的情况下，地震波正演数据可以通过式(1)得到：

d=L(m)

(1)

式中L——正演算子，根据模型的不同可以采用不同的正演算子，如声波算子、弹性波算子、黏弹各向异性算子等，算子越复杂，得到的正演数据中的波场也就越复杂；

d——从模型得到的数据，也可以认为是野外地震采集到的数据，这个数据的大小和质量由设定的采集观测系统或者说采集参数所决定；

m——地质模型，也就是我们常说的地下构造，其中包含我们需要勘探的地质目标。

一般认为，采集到的地震数据资料越丰富，通过地震资料处理后得到的内容也就越完全，但采集成本也就越高。

Tarantola认为，可以通过地震数据恢复出地震模型，但我们知道，如果地震数据中包含的信息极少，就算再好的反演方法也不能得到满意的结果。目前，几乎所有全波形反演研究人员都把精力集中在全波形反演方法的本身，而没有考虑到数据采集系统对反演结果的影响。在已有地震数据的情况下，利用式(2)可以求出地质模型：

m=L-1(d)

(2)

式中L-1——作用在地震数据上的反演算子，表示利用数学算法恢复模型的过程。

模型速度恢复的准确度直接决定了叠前深度偏移成像的结果。

对于采集观测系统设计，不考虑模型的体积模量和密度信息，使用简化的标量声波方程(3)进行全波形反演，可以节省大量计算时间：

(3)

式中v(x,y,z)——模型的纵波速度；

u(x,y,z,t)——空间不同时刻的声波压力场；

f(x,y,z,t)——地震震源。

将方程加入边界条件和初值条件，离散后会得到相应的差分方程。利用有限差分法或者有限元法能够对方程进行求解，并得到不同位置、不同时刻的空间位置的波场。用准确模型计算得出的声波压力场u(x,y,z,t)，可以近似地认为它就是野外观测到的地震记录，记为dobs，用初始模型或者逐次迭代得到的模型进行计算得出的u(x,y,z,t)可以记为ucal。一般来说，常规反射地震勘探方法的检波点都位于地表或者接近于地表，所以u(x,y,z,t)变化为u(x,y,0,t)。由于只记录了地表位置的波场，而缺少大部分其他空间位置的波场，这给地下地质构造的反演造成了很大的困难。

定义初始模型波场记录与观测到地震之间的差值为：

δu=ucal-dobs

(4)

对于得到的地震记录可以通过最小二乘方法进行求解。定义全波形反演的目标函数为：

(5)

式中S——炮点位置；

R——检波点位置；

T——总的记录时间；

EN——可视为残差函数的L2范数。

方程(5)的最优化解受炮点位置、检波点位置、记录时间等参数的约束。要使目标函数达到最小，就要用初始模型模拟得到的ucal与实际数据dobs进行对比，并修改初始模型，用更新后的模型做正演模拟得新一轮的ucal，再对比，再更新速度模型，直至方程(5)的值达到最小，也就是速度模型和实际模型最为接近。理论上，这是一个迭代收敛的过程，但实际上很多原因会使上述迭代过程不能收敛。例如数据量不足就会使反演不收敛，故一个好的观测系统参数是地震全波形反演能够收敛到好的结果的必要前提。另外，实际数据中缺少低频信息也会造成迭代陷入局部极小，不会收敛到真实的模型。很多学者开展了缺少低频信息数据的全波形反演方法的研究。Wu等提出了利用包络反演的方法进行初始模型的建立，罗静蕊等[14-15]研究了地震包络反演对局部极小值的抑制特性，陈生昌等[16]使用时间积分波场进行全波形反演期望获得良好的初始速度场，胡勇等[17]进行了基于精确震源函数的解调包络多尺度全波形反演。此部分内容与地震采集观测系统的关系不大，这里不再赘述。

反演速度逐次迭代过程中的一个非常重要的步骤是速度模型的更新，这个步骤的核心是梯度场的计算，为了提高计算效率，这里避免Ferchet矩阵导数的计算，而采用伴随方法计算梯度。这种方法利用炮点正传波场与检波点残差的逆时传播波场进行零延迟互相关，逐步累加形成残差函数的梯度。伴随方法计算梯度的计算公式如下：

δv(x,y,z)=∑T[uf(x,y,z,t)ub(x,y,z,t)]

(6)

式中δv(x,y,z)——速度更新量；

uf(x,y,z,t)——不同时刻炮点正向传播波场；

ub(x,y,z,t)——不同时刻检波点数据残差的逆时传播波场。

伴随方法计算梯度的过程中使用了炮点正传波场与检波点残差的逆时传播波场的互相关，这个步骤把地震波入射信息与反射信息结合起来，形成了反射全波形反演。同样的技术也在叠前深度偏移成像的过程中进行类似的应用，但叠前深度偏移主要利用波形的相位信息，达到同相叠加的目的，最终形成以反射同相轴为反射信息的地震偏移成像结果，其成像分辨率主要受偏移速度模型、地震子波和最大偏移距的影响。如图1所示，KS为炮点在地下某空间位置的炮点波数，KR为同一点处的检波点波数，其合成矢量为KSR，其展布范围体现此炮检对的分辨率信息，垂向展布范围体现垂向分辨率，横向展布范围表示横向分辨率，而对应于观测系统参数中的最大偏移距主宰了局部反射的张角信息，其张角越大，成像分辨率越高，但仍受地震子波波长的影响。

图1 反射点处的局部波数矢量Fig.1 Local wave number vector at the reflection point

KSR=KS-KR

(7)

对于全波形反演方法，由于最终反演得到的是速度模型，而不是反射波同相轴的地震剖面，所以它不受地震波长的限制，其分辨率要比地震偏移成像的分辨率高得多，在此基础上，要获得良好的速度模型分辨率，就要求有比叠前偏移更长的最大偏移距。这也是地震全波形反演比常规叠前偏移成像对采集观测系统要求更高的原因之一。为此，这里发展了面向全波形反演的地震采集观测系统设计优化研究。

2 模型试算

由于三维波动正演与全波形反演的计算量巨大，为了阐明道理，我们这里只使用二维模型进行数值试验。图2为国内某工区模型，模型深度为3 500 m，模型长度为15 000 m，模型速度范围从1 800 m/s到3 500 m/s。模型中包含5个层位，并被多个断层所截断，由于断距较大，因此可以认为是横向速度变化较大的地质结构。后面的数值计算结果都基于此模型。

图2 国内某工区速度模型Fig.2 Velocity model of a certain area in China

由于全波形反演对最大偏移距信息比较敏感，因此首先分析不同的最大偏移距对反演模型的影响。为了提高计算速度，将模型离散化成25 m×25 m，在模型地表共布设81炮，炮点范围从2 500 m 至 12 500 m，炮点间距为125 m，检波点位于炮点两侧接收，检波点间距为25 m，设计最大偏移距分别为500 m、1 000 m、1 500 m、2 000 m和2 500 m五种观测系统方案，由于模型横向范围只有15 000 m，考虑到满覆盖区域的问题，因此没有设置更大最大偏移距的滚动接收方案。

使用以上5种观测方案进行正演模拟得到合成记录，我们粗略地认为它们就是野外观测到的地震记录。在此我们使用标量声波方程，空间10阶差分精度，时间2阶差分精度，2 ms采样间隔，6 s记录长度。由于模型网格尺寸的限制，为了避免数值频散，震源子波采用9 Hz的雷克子波，同时为保证记录范围内波形的完整性，子波的延迟时采用222 ms。

图3中分别为模型8 000 m处不同偏移距的有限差分正演模拟记录，从地震记录中可以看出，随着最大偏移距的增加，单炮记录内的信息逐渐增多，特别是复杂构造断面波和绕射波都可以在单炮记录中出现。在500 m最大偏移距的情况下，深层的断面波没有出现，一方面是由于接收排列太短没有接收到，另一方面是由于正演模拟算法使用的模型范围过小。

图3 不同最大偏移距单炮有限差分模拟数据Fig.3 Finite difference data of single shots with different maximum offsets

全波形反演过程中，初始模型的选择会影响反演的收敛速度，很多学者已经进行了相关的研究，本文不做重点分析。我们对原始模型使用1 000 m的平滑半径进行平滑，得到的模型作为初始模型(图4)进行反演。可以看到，相对于原始模型(图2)，初始模型中高波数成分已经丢失，只含有低波数成分，对于这样的初始模型，只要观测系统中含有足够的信息就可以反演出较为精细的速度模型。

图4 使用原始模型平滑后的反演初始模型Fig.4 Inversion initial model after smoothing with the original model

图5是使用不同最大偏移距观测方案与初始模型(图4)进行全波形反演迭代39次得到的结果。反演过程中约束了最小速度1 800 m/s与最大速度3 500 m/s。能够看出，当最大偏移距为500 m时，经过全波形反演，不但没有收敛到精细的模型，反而破坏了原有低波数模型的准确度，尤其是在模型的中深部产生了类似于波场传播的噪音。随着最大偏移距的增加，反演对模型的破坏越来越小，在最大偏移距为2 000 m时，深部构造逐渐显现，看似噪音，但与正确模型相比，其构造是正确的。

当偏移距达到2 500 m时，其结果已经相当令人满意，也没有破坏初始模型，整个过程一直处于收敛状态，这说明此模型在此种条件下，2 500 m应该是保证全波形反演能够收敛的最小的“最大偏移距”，低于此值，将不能保证全波形反演的收敛。这里的数值例子也从侧面证明了很多实际资料很难完成全波形反演任务的原因之一就是观测系统不能满足全波形反演对最大偏移距的要求。

由于模型横向长度的限制，同时考虑到覆盖次数的均匀，我们只算到了2 500 m最大偏移距的滚动排列方式。为了说明最大偏移距的重要性，我们采用全模型范围激发和接收，其中炮点间隔为125 m，检波点间隔为25 m，共计121炮，其他正演模拟参数与反演参数同上。图6为全模型范围激发与接收39次迭代的全波形反演结果。其反演结果质量要远远好于2 500 m最大偏移距方案，图中断层界面清晰可见，与准确模型之间的差别已经很小，可以用作叠前时间偏移的初始模型。

图6 全模型接收时的全波形反演结果Fig.6 FWI results with the offset covers the whole velocity model

为了说明本文提出方法的有效性，我们与目前工业界使用的方法进行对比。在面向复杂模型的观测系统设计中，现在最流行的方法就是波动照明分析及逆时偏移成像的观测系统分析。对于上述使用的5种观测系统，我们分别从波动照明分析和逆时偏移的分析思路进行计算。

图7中展示的是不同方案的照明分析结果。从结果可以看出，随着最大偏移距的增加，不同观测方案对应的照明能量也逐渐增加，对应的观测系统能够探测的深度也在增加，但我们无法判断出哪种观测系统能够满足我们的勘探要求。

图7 不同最大偏移距观测系统波动照明分析结果Fig.7 Illumination results of surveys with different maximum offset

图8是图2所示模型中最下层(L6)对应的不同观测方案下的照明曲线，通过该曲线也能看出，随着最大偏移距的增加，照明能量也在逐渐增强，但同样无法说明哪个方案的观测系统能够满足要求。

图8 不同观测方案L6层照明能量对比Fig.8 Comparison of illumination energy of layer L6 in different surveys

图9中展示了不同最大偏移距方案的逆时偏移结果，其中偏移速度模型使用准确速度模型，随着最大偏移距的逐步增加，逆时偏移结果逐渐变好，由于500 m 的最大偏移距过小，在偏移剖面中存在大量的噪音，同时断点位置不是很清晰，并且由于逆掩推覆构造的阻挡作用，弱照明区域的成像也不是很好。但所有构造都归位到正确的位置。当最大偏移距增加到1 500 m时，逆时偏移成像结果已经很好，断点和断面都比较清晰，下伏岩层成像准确。当继续增加最大偏移距时，成像结果继续变好，但变好程度有限。这说明，当偏移速度模型准确时，叠前深度偏移成像对数据的要求并不高，利用叠前偏移成像来分析观测系统不是一个有效的手段。而实际资料很难获得准确成像的一个主要原因就是很难获得准确的速度模型，所以速度模型的准确建立要比常规叠前偏移成像的要求更高，而全波形反演恰恰是一种很有潜力的速度模型建立手段。若想获得良好的地震资料，需要这种面向全波形反演的地震采集观测系统设计与优化的手段。

图9 不同最大偏移距方案的逆时偏移结果Fig.9 Reverse time migration results for different maximum offsets

3 结束语

地震勘探采集设计中，需要对采集参数进行论证，在对最大偏移距的论证过程中，常常采用分偏移距叠加的思路判断最大偏移距是否能够满足要求，但在对地质目标精度逐渐提高的今天，水平叠加的思路已经远远不能满足要求，应运而生的面向复杂模型的波动照明分析在参数论证上发挥了巨大的作用，但照明分析也存在其自身特有的弱点，它只体现了能流密度，没有充分考虑到反射波动相位信息，所以不能真正地体现偏移成像，这使得波动照明分析很难在精细勘探中发挥作用。而叠前深度偏移成像分析没有考虑到速度模型的建立，在准确速度模型的条件下，会使最大偏移距、道距、炮密度等参数的敏感性较差。

本文提出将全波形反演应用于采集观测系统设计，由于全波形反演的高分辨率特性，使得全波形反演对采集观测系统的要求有了进一步的提高，尤其在复杂高陡构造区，基于全波形反演的观测系统设计将会给采集观测系统设计打开一扇大门，为复杂区观测系统设计提供有力的支撑。

由于篇幅限制，本文只论证了最大偏移距的影响，其完全可以用于论证分析道距、炮距、炮道密度等与观测系统相关的其他各项参数；在三维情况下，还可以分析观测系统的纵横比、观测方位等参数。