说明书地震观测系统优化设计的层状介质双聚焦方法及其应用
技术领域
本发明涉及地震勘探技术领域，特别面向复杂地区特定勘探目标的三维观测系统优化设计，具体是一种地震观测系统优化设计的层状介质双聚焦方法及其应用。
背景技术
激发、接收和观测系统是地震资料采集的3个主要组成部分。三维观测系统设计水平，直接影响地下构造的成像质量，同时还直接决定野外施工的工程总量，因此它是三维地震资料采集技术设计中的第一重点。
目前在常规三维地震勘探中普遍采用的优化设计方法是CMP面元属性分析法。这种方法是以地下共中心点(CMP)水平叠加理论为基础，通过优化CMP面元的纵横向覆盖次数、炮检距分布和炮检线方位角分布等特性参数，以求在数据处理中取得最佳的CMP叠加和叠后偏移成像结果。这种优化方法是宏观的、全局性的，所以它在以往的常规三维地震勘探中发挥了极其重要的作用。但面对的是复杂地区复杂地质目标的油气勘探，相应地要求将采集设计重点集中到局部区域和局部目标上来，研究开发另一类以地下共反射点(CRP)为对象的，适合叠前深度偏移和CRP偏移叠加处理的方法。近年来，主要研究集中在三维射线追踪法和双聚焦特征函数分析法。
三维射线追踪法是一种首选的方法，因为二维射线追踪法的具体方案很多，算法也比较成熟，多年应用积累了许多成功的经验。但是在研究和应用过程中发现，三维的情况比二维复杂得多，由于三维射线追踪本身就比较困难，如果再遇到复杂地质模型，往往会造成追踪失败，所以至今尚未有成功的工业化应用。
双聚焦理论是荷兰Delft科技大学Berkhout在20世纪80年代提出的。起初主要用于地震偏移成像技术上，20世纪末，他与Volker开始将均匀介质模型的点散射双聚焦成像原理应用到三维观测系统优化设计中去，通过估算、分析3个双聚焦特征函数(分辨率函数、AVP函数和照明强度分布切片)，完成三维观测系统优化设计这一目标。发表了4篇论文(A.J.Berkhout，L.Ongkiechong.Analysis of seismic acquisition geometries byfocal beams.SEG 1998 Expanded Abstracts。A.W.E.Volker，G.Blaequiere.Analysisof 3-D seismic acquisition geometries by focal beams，SEG 1999 Expanded Abstracts。A.W.E.Volker etc.Analysis and Optimization of 3-D seismic acquisition geometriesby focal beams.SEG 2000 Expanded Abstracts。A.W.E.Volker etc.Acquisition designfor optimum ampiitude accuracy.SEG 2001 Expanded Abstracts.)。
理论上三维观测系统优化设计的双聚焦分析法可适用于层状介质和复杂介质模型，但目前只完成了均匀介质双聚焦分析法的研究，开发了均匀介质双聚焦的工业化软件(狄帮让等.三维观测系统优化设计的双聚焦方法.石油地球物理勘探，2003，38(5)：463～469)。
上述三种分析法的主要特点汇集于下表：面元属性分析法三维射线追踪法双聚焦分析法
共中心点CMP共反射点CRP共反射点CRP
全局性的面向复杂地质目标地面向复杂地质目标的
均匀介质水平层状介质均匀介质层状介质和复杂介质均匀介质层状介质和复杂介质
最佳水平叠加最佳叠前偏移成像最佳叠前和叠后偏移成像
定性的定量的定量的
均匀介质双聚焦特征函数计算方法可以简述如下：
在空间频率域，均匀介质一次反射波的WRW正演模型可用图1和(1)式表示，
P ( Z 0 ) = D ( Z 0 ) W D ( Z 0 , Z m ) R ( Z m ) W S ( Z m , Z 0 ) S ( Z 0 ) - - - ( 1 ) ]]>
其中：P—数据矩阵(地震记录)，WD—上行传播算子矩阵，WS—下行传播算子，R—反射系数矩阵，D—检波点矩阵，S—震源矩阵，Z—深度坐标，Z0—地表，Zm—地下目标层。
设FD(Zm，Z0)表示检波聚焦算子，FS(Z0，Zm)表示震源聚焦算子，定义聚焦检波束BD(Zm，Z0)和聚焦震源束BS(Z0，Zm)如下：
B D ( Z m , Z 0 ) = F D ( Z m , Z 0 ) D ( Z 0 ) W D ( Z 0 , Z m ) - - - ( 2 ) ]]>
B S ( Z 0 , Z m ) = W S ( Z m , Z 0 ) S ( Z 0 ) F S ( Z 0 , Z m ) - - - ( 3 ) ]]>
将聚焦算子作用于地震记录P(Z0)，可以得到FPF双聚焦成像结果ΔP(Zm)，即
ΔP ( Z m ) = F D ( Z M , Z 0 ) P ( Z 0 ) F S ( Z 0 , Z M ) ]]>
ΔP ( Z m ) = B D ( Z m , Z 0 ) R ( Z m ) B S ( Z 0 , Z m ) ]]>
如果FD(Zm，Z0)和FS(Z0，Zm)均为优化的聚焦算子，则散射点i(正演目标点)与聚焦点j(反向成像点)重合，聚焦束BD(Zm，Z0)和BS(Z0，Zm)就皆为单位矩阵，于是 ΔP ( Z m ) = R ( Z m ) , ]]>实现了理想的反射系数成像。如果聚焦算子有偏差，则聚焦点就会偏离散射点，成像结果也会发生偏差。聚焦成像的分辨率函数定义为下式：
RF ( Z m ) = Σ B D ( Z m ) B S ( Z m ) - - - ( 4 ) ]]>
式中∑表示对聚焦孔径内所有观测系统模板求和。分辨率函数定量地描述了在地面特定的观测系统下，成像点与真实目标点位置偏差对偏移成像结果的影响。
AVP(Amplitude versus ray-parameter)函数描述了振幅与射线关系的参数信息。在空间频率域，AVP函数可由聚焦检波束与聚焦震源束的褶积得到；在波数频率域，则可由他们的乘积得到。
AVP ( Z m ) = Σ B D ( Z m ) * B S ( Z m ) - - - ( 5 ) ]]>
沿地下目标层，求取每个目标面元的分辨率函数极大值，即可得到照明强度分布切片，也即采集痕迹。
综上所述，双聚焦的3个特征函数都取决于聚焦束，所以，如何求取聚焦检波束BD(Zm，Z0)和聚焦震源束BS(Z0，Zm)，这是双聚焦算法研究中的关键。
但是，均匀介质是理想的模型，它与实际地质构造相差甚远，实际构造通常是层状介质模型，或者是复杂介质模型。实践表明，均匀介质双聚焦分析法不能应用于岩性勘探、储层勘探等复杂地质目标勘探的三维地震采集设计。
发明内容
本发明的目的是提供一种层状介质双聚焦方法，以及利用上述方法在复杂地区特定勘探目标的三维地震观测系统优化设计。
本发明由以下步骤建立层状介质双聚焦模型：
1)从地表震源到地下目标点的正向传播效应WS(Zm，Z0)用一系列的层间传播算子(矩阵)相乘WS(zm，zm-1)WS(zm-1，zm-2)…WS(z1，z0)来实现；
上式中Z0为地表，地下各层为Z1、Z2…Zm，
Zm为目标层；
2)从地下目标点到地表检波排列片的正向传播效应WD(z0，zm)用一系列的层间传播算子(矩阵)相乘WD(z0，z1)WD(z1，z2)…WD(zm-1，zm)来实现；
3)同理，从地下目标点到地表震源的反向聚焦震源效应及从地表检波排列片到地下目标点的反向聚焦检波效应也用一系列的层间聚焦算子(矩阵)相乘来实现；
4)由步骤1)-3)即可套入均匀介质双聚焦计算公式求得聚焦检波束和聚焦震源束，紧接着可求得层状介质双聚焦的3个特征函数，
由以下步骤实现的层状介质双聚焦计算方法：
1)从地表每一个观测系统模板的震源排列出发，先模拟正向传播过程，将压力波场从地表逐层外推到地下目标层的目标反射点；
2)以地下目标点为中心，确定聚焦成像平面，随后模拟反向聚焦过程，将压力波场由成像平面逐层外推到地表震源排列，于是得到单模板的聚焦震源束；
3)将地下目标反射点视为散射点，从散射点出发，模拟正向传播过程，将压力波场由地下逐层外推到地表检波器排列；
4)随后模拟反向聚焦，将压力波场从地表逐层外推到地下成像平面，得到单模板的聚焦检波束；
5)将单模板的聚焦震源束和聚焦检波束相乘，得到单模板的分辨率函数；
6)对所有的单模板分辨率函数求和得到双聚焦成像分辨率函数；
7)利用二维Fourier变换，将单模板的聚焦震源束和聚焦检波束从空间域(x，y，z)变换到波束域(kx，ky，z)，随后两者相乘即可得到单模板的AVP函数；
8)对所有的单模板AVP函数求和，得到双聚焦成像AVP函数；
9)沿地下目标层，求取每个反射面元的分辨率函数极大值，即可得到照明强度分布切片。
本发明层状介质双聚焦算法模型实现正向传播效应和反向聚焦效应为：
1)从地表震源到地下目标点的正向传播效应为
W S ( z m , z 0 ) = w S ( z m , z m - 1 ) W S ( z m - 1 , z m - 2 ) · · · W S ( z 1 , z 0 ) ]]>
上式中Z0为地表，地下各层为Z1、Z2…Zm，Zm为目标层；
2)从地下目标点到地表检波排列片的正向传播效应为
W D ( z 0 , z m ) = W D ( z 0 , z 1 ) W D ( z 1 , z 2 ) · · · W D ( z m - 1 , z m ) ]]>
3)从地下目标点到地表震源的反向聚焦震源效应FS(z0，zm)及从地表检波排列片到地下目标点的反向聚焦检波效应FD(zm，z0)为
F S ( z 0 , z m ) = F S ( z 0 , z 1 ) F S ( z 1 , z 2 ) · · · F S ( z m - 1 , z m ) ]]>
F D ( z m , z 0 ) = F D ( z m , z m - 1 ) F D ( z m - 1 , z m - 2 ) · · · F D ( z 1 , z 0 ) ]]>
本发明将模拟波场逐层正向传播和模拟波场逐层反向聚焦的递推分解为两个聚焦束实现，其中聚焦震源束为：
下行波正向传播：从地表(z0)观测系统单模板的震源排列出发，每个震源激发的压力波场PS(x0，y0，z0)皆设为1，由上向下将压力波场逐层递推到地下目标层zm的目标点i，得到叠加波场PSi(zm)，
P S ( x 1 , y 1 , z 1 ) = ∫ ∫ W S ( z 1 , z 0 ) P S ( x 0 , y 0 , z 0 ) S S ( x 0 , y 0 , z 0 ) dxdy ]]>
P S ( x 2 , y 2 , z 2 ) = ∫ ∫ W S ( z 2 , z 2 ) P S ( x 1 , y 1 , z 1 ) S ( x 1 , y 1 , z 1 ) dxdy ]]>
P Si ( z m ) = ∫ ∫ W Si ( z m , z m - 1 ) P S ( x m - 1 , y m - 1 , z m - 1 ) S ( x m - 1 , y m - 1 , z m - 1 ) dxdy ]]>
反向聚焦：成像平面网格以地下目标点i为中心，每个成像网格单元j的压力波场PSj(zm)皆设为1，从网格j出发，由下向上将压力波场逐层聚焦到地表的震源排列，得到聚焦震源束BSij(zm)，
B Sj ( x m - 1 , y m - 1 , z m - 1 ) = P Sj ( z m ) S j ( z m ) F S ( z m - 1 , z m ) ]]>
B Sj ( x m - 2 , y m - 2 , z m - 2 ) = ∫ ∫ P Sj ( x m - 1 , y m - 1 , z m - 1 ) S ( x m - 1 , y m - 1 , z m - 1 ) F S ( z m - 2 , z m - 1 ) dxdy ]]>
B Sj ( x 0 , y 0 , z 0 ) = ∫ ∫ P Sj ( x 1 , y 1 , z 1 ) S ( x 1 , y 1 , z 1 ) F S ( z 0 , z 1 ) dxdy ]]>
B Sij ( z m ) = ∫ ∫ P Si ( z m ) B Sj ( x 0 , y 0 , z 0 ) S S ( x 0 , y 0 , z o ) dxdy ]]>
以上SS(x0，y0，z0)为地表震源算子，S(x，y，z)为各层的网格单元算子，均为1；聚焦检波束为：上行波正向传播：从地下zm层的目标散射点i出发，压力波场PDi(Zm)设为1，由下向上逐层递推到地面观测系统模板的每个检波点，得到波场PD(x0，y0，z0)，
P D ( x m - 1 , y m - 1 , z m - 1 ) = P Di ( Z m ) S i ( z m ) W D ( z m - 1 , z m ) ]]>
P D ( x m - 2 , y m - 2 , z m - 2 ) = ∫ ∫ P D ( x m - 1 , y m - 1 , z m - 1 ) S ( x m - 1 , y m - 1 , z m - 1 ) W D ( z m - 2 , z m - 1 ) dxdy ]]>
P D ( x 0 , y 0 , z 0 ) = ∫ ∫ P D ( x 1 , y 1 , z 1 ) S ( x 1 , y 1 , z 1 ) W D ( z 0 , z 1 ) dxdy ]]>
反向聚焦：从单模板的每个检波点SD(x0，y0，z0)出发，由上向下将压力波场PD(x0，y0，z0)逐层聚焦到zm层成像网格单元j，得到聚焦检波束BSij(zm)，
B D ( x 1 , y 1 , z 1 ) = ∫ ∫ F D ( z 1 , z 0 ) P D ( x 0 , y 0 , z 0 ) S D ( x 0 , y 0 , z 0 ) dxdy ]]>
B D ( x 2 , y 2 , z 2 ) = ∫ ∫ F D ( z 2 , z 1 ) B D ( x 1 , y 1 , z 1 ) S ( x 1 , y 1 , z 1 ) dxdy ]]>
B Dij ( z m ) = ∫ ∫ F Dj ( z m , z m - 1 ) B D ( x m _ 1 , y m - 1 , z m - 1 ) S ( x m - 1 , y m - 1 , z m - 1 ) dxdy ]]>
以上SD(x0，y0，z0)为地表检波算子。
以上SD(x0，y0，z0)为地表检波算子，BSij(zm)为聚焦震源束，S(x，y，z)为各层的网格单元算子，均为1。
利用本发明地震观测系统优化设计的层状介质双聚焦方法，结合常规CMP面元属性分析法的三维观测系统优化设计应用，步骤如下：
1)输入若干个三维观测系统；
2)根据输入的三维观测系统，计算其CMP面元的3个属性参数，即纵横向覆盖次数、炮检距分布和炮检线方位角分布，不同的观测系统得到不同的CMP面元属性参数；
3)分析对比不同CMP面元属性参数对CMP叠加质量的影响，以此优选有利于水平叠加的最佳观测系统；
4)输入层状地质模型及主要地质目标；
5)选择第一个三维观测系统，从地表震源排列出发，先模拟正向传播过程，将压力波场从地表逐层外推到地下目标层的目标反射点；
从地表(z0)观测系统单模板的震源排列出发，每个震源激发的压力波场PS(x0，y0，z0)皆设为1，由上向下将压力波场逐层递推到地下目标层zm的目标点i，得到叠加波场PSi(Zm)，
P S ( x 1 , y 1 , z 1 ) = ∫ ∫ W S ( z 1 , z 0 ) P S ( x 0 , y 0 , z 0 ) S S ( x 0 , y o , z 0 ) dxdy ]]>
P S ( x 2 , y 2 , z 2 ) = ∫ ∫ W S ( z 2 , z 2 ) P S ( x 1 , y 1 , z 1 ) S ( x 1 , y 1 , z 1 ) dxdy ]]>
P Si ( z m ) = ∫ ∫ W Si ( z m , z m - 1 ) P S ( x m - 1 , y m - 1 , z m - 1 ) S ( x m - 1 , y m - 1 , z m - 1 ) dxdy ]]>
6)以地下目标点为中心，确定聚焦成像平面，随后模拟反向聚焦过程，将压力波场由成像平面逐层外推到地表震源排列，于是得到聚焦震源束；
从地下zm层的目标散射点i出发，压力波场PDi(Zm)设为1，由下向上逐层递推到地面观测系统模板的每个检波点，得到波场PD(x0，y0，z0)，
P D ( x m - 1 , y m - 1 , z m - 1 ) = P Di ( Z m ) S i ( z m ) W D ( z m - 1 , z m ) ]]>
P D ( x 0 , y 0 , z 0 ) = ∫ ∫ P D ( x 1 , y 1 , z 1 ) S ( x 1 , y 1 , z 1 ) W D ( z 0 , z 1 ) dxdy ]]>
7)从目标反射点(散射点)出发，模拟正向传播过程，将压力波场由地下逐层外推到地表检波器排列；
由地下zm层的目标散射点i出发，压力波场PDi(Zm)为1，由下向上逐层递推到地面检波点PD(x0，y0，z0)
P D ( x m - 1 , y m - 1 , z m - 1 ) = P Di ( Z m ) W D ( z m - 1 , z m ) S ( x m - 1 , y m - 1 , z m - 1 ) ]]>
P D ( x m - 2 , y m - 2 , z m - 2 ) = ∫ ∫ P D ( x m - 1 , y m - 1 , z m - 1 ) W D ( z m - 2 , z m - 1 ) S ( x m - 2 , y m - 2 , z m - 2 ) dxdy ]]>
P D ( x 0 , y 0 , z 0 ) = ∫ ∫ P D ( x 1 , y 1 , z 1 ) W D ( z 0 , z 1 ) S D ( x 0 , y 0 , z 0 ) dxdy ]]>
8)随后模拟反向聚焦，将压力波场从地表逐层外推到地下成像平面，于是得到聚焦检波束；
从单模板的每个检波点SD(x0，y0，z0)出发，由上向下将压力波场PD(x0，y0，z0)逐层聚焦到zm层成像网格单元j，得到聚焦检波束BSij(zm)，
B D ( x 1 , y 1 , z 1 ) = ∫ ∫ F D ( z 1 , z 0 ) P D ( x 0 , y 0 , z 0 ) S D ( x 0 , y 0 , z 0 ) dxdy ]]>
B D ( x 2 , y 2 , z 2 ) = ∫ ∫ F D ( z 2 , z 1 ) B D ( x 1 , y 1 , z 1 ) S ( x 1 , y 1 , z 1 ) dxdy ]]>
B Dij ( z m ) = ∫ ∫ F Dj ( z m , z m - 1 ) B D ( x m _ 1 , y m - 1 , z m - 1 ) S ( x m - 1 , y m - 1 , z m - 1 ) dxdy ]]>
以上SD(x0，y0，z0)为地表检波算子，BSij(zm)为聚焦震源束，S(x，y，z)为各层的网格单元算子，均为1；
9)由聚焦震源束和聚焦检波束可以直接求得该观测系统对应的分辨率函数、AVP函数和照明强度分布切片；
10)选择第二个、第三个三维观测系统，重复上述第5)-9)步计算过程，不同的观测系统得到不同的双聚焦特征函数；
11)以指定地质目标的最佳偏移成像为评价标准，对计算得到的几组双聚焦特征函数进行分析对比，优选三维观测系统；
12)以双聚焦分析法为主，参考CMP面元分析法，最后确定最佳三维观测系统。
附图说明
图1是均匀介质一次反射波的WRW正演模型；
图2显示了本发明的层状介质聚焦震源束和聚焦检波束的推算方法；
图3更具体地显示了本发明的聚焦检波束计算步骤和计算方法；
图4是本发明双聚焦法三维观测系统优化设计流程图；
图5显示了三个等速模型；
图6显示了三个模型分辨率函数的水平切片；
图7显示了三个模型分辨率函数的垂直切片；
图8显示了一个具有两个倾斜界面的三维地质模型；
图9显示了不同倾角成像面的分辨率函数；
图10显示了不同倾角成像面的空间分辨率垂直切片；
图11是宽窄两种三维观测系统CMP面元属性对比图；
图12是宽窄两种三维观测系统双聚焦特征函数对比图；
下面结合附图及具体实例对本发明进行详细说明。
具体实施方式
本发明由以下两部分组成：
1)层状介质双聚焦特征函数计算方法；
2)以层状介质双聚焦特征函数分析法为主结合常规CMP面元属性分析法为辅的三维观测系统优化设计方法。
层状介质双聚焦算法的基础是均匀介质双聚焦算法，其模型如下：
地表为Z0，地下各层依次为Z1、Z2…Zm，其中Zm为目标层；层与层之间的下行传播算子用WS(Z1，Z0)、WS(Z2，Z1)…WS(Zm，Zm-1)表示，上行传播算子用WD(Zm-1，Zm)…WD(Z1，Z2)、WD(Z0，Z1)表示，显然层与层之间的两种聚焦算子也可以此类推。
层状介质双聚焦算法的基础是均匀介质双聚焦算法，其模型如下：
1)从地表震源到地下目标点的正向传播效应WS(Zm，Z0)可用一系列的层间传播算子(矩阵)相乘WS(zm，zm-1)WS(zm-1，zm-2)…WS(z1，z0)来实现，即
W S ( z m , z 0 ) = W S ( z m , z m - 1 ) W S ( z m - 1 , z m - 2 ) · · · W S ( z 1 , z 0 ) - - - ( 6 ) ]]>
2)从地下目标点到地表检波排列片的正向传播效应WD(Z0，Zm)也可用一系列的层间传播算子(矩阵)相乘WD(z0，z1)WD(z1，z2)…WD(zm-1，zm)来实现，即
W D ( z 0 , z m ) = W D ( z 0 , z 1 ) W D ( z 1 , z 2 ) · · · W D ( z m - 1 , z m ) - - - ( 7 ) ]]>
3)同理，从地下目标点到地表震源的反向聚焦震源效应及从地表检波排列片到地下目标点的反向聚焦检波效应也可用一系列的层间聚焦算子(矩阵)相乘来实现，即
F D ( z m , z 0 ) = F D ( z m , z m - 1 ) F D ( z m - 1 , z m - 2 ) · · · F D ( z 1 , z 0 ) - - - - ( 8 ) ]]>
F S ( z 0 , z m ) = F S ( z 0 , z 1 ) F S ( z 1 , z 2 ) · · · F S ( z m - 1 , z m ) - - - ( 9 ) ]]>
4)由(6)-(9)式即可套入均匀介质双聚焦公式(2)和(3)，求得聚焦检波束和聚焦震源束如下，
B D ( z m , z 0 ) = ∫ ∫ F D ( z m , z 0 ) S D ( z 0 ) W D ( z 0 , z m ) d x D d y D - - - - ( 10 ) ]]>
B S ( z 0 , z m ) = ∫ ∫ W S ( z m , z 0 ) S S ( z 0 ) F S ( z 0 , z m ) d x S d y S - - - - ( 11 ) ]]>
其中：SD为地表检波采样算子，SS为地表震源采样算子，其值在检波点或震源点处为1，其他位置均为0。另外，聚焦算子可取传播算子的共轭矩阵，即F＝W*；
5)紧接着可求得3个特征函数。
上述为理论模型，目的为了证明层状介质算法与均匀介质算法的关联性。实际可以实现的算法如下：将每一层内部视为均匀介质，于是可以采用均匀介质的传播算子和聚焦算子进行波场外推；在层与层之间，则是采用逐层递推的递归方法，应用波动方程积分解，完成正反向波场外推(参考图2)。
下面以聚焦检波束为例说明具体的计算公式及计算步骤(参考图3)。
正向传播：由地下zm层的目标散射点i出发，将压力波场PDi(zm)(可设为1)由下向上逐层递推到地面检波点PD(x0，y0，z0)。对于平界面，选择Rayleigh积分公式；对于弯曲界面，采用Kirchhoff积分解进行外推。
P D ( x m - 1 , y m - 1 , z m - 1 ) = P Di ( Z m ) S i ( z m ) W D ( z m - 1 , z m ) ]]>
P D ( x m - 2 , y m - 2 , z m - 2 ) = ∫ ∫ P D ( x m - 1 , y m - 1 , z m - 1 ) S ( x m - 1 , y m - 1 , z m - 1 ) W D ( z m - 2 , z m - 1 ) dxdy ]]>
P D ( x 0 , y 0 , z 0 ) = ∫ ∫ P D ( x 1 , y 1 , z 1 ) S ( x 1 , y 1 , z 1 ) W D ( z 0 , z 1 ) dxdy ]]>
反向聚焦：从地面检波点出发，将压力波场PD(x0，y0，z0)逐层聚焦到zm层成像平面网格j，j可以在成像平面上围绕i点移动。
B D ( x 1 , y 1 , z 1 ) = ∫ ∫ F D ( z 1 , z 0 ) P D ( x 0 , y 0 , z 0 ) S D ( x 0 , y 0 , z 0 ) dxdy ]]>
B D ( x 2 , y 2 , z 2 ) = ∫ ∫ F D ( z 2 , z 1 ) B D ( x 1 , y 1 , z 1 ) S ( x 1 , y 1 , z 1 ) dxdy ]]>
B Dij ( z m ) = ∫ ∫ F Dj ( z m , z m - 1 ) B D ( x m _ 1 , y m - 1 , z m - 1 ) S ( x m - 1 , y m - 1 , z m - 1 ) dxdy ]]>
BSij(zm)为聚焦震源束，S(x，y，z)为各层的网格单元算子，均为1。
上两式中，积分区域可选择有限的孔径范围进行计算，也可扩大到整个模型边界；SD(x0，y0，z0)为地面检波算子，网格大小b应满足空间采样定理，即b≤Vi/2fd，其中：Vi为层速度，fd为压力波主频。
同理可以计算聚焦震源束BSij(zm)，随之可以求得分辨率函数，AVP函数和照明强度分布切片。
本发明的层状介质双聚焦特征函数计算方法包括以下步骤(参考图4)：
1)从地表每一个观测系统模板的震源排列出发，先模拟正向传播过程，将压力波场从地表逐层外推到地下目标层的目标反射点；
2)以地下目标点为中心，确定聚焦成像平面，随后模拟反向聚焦过程，将压力波场由成像平面逐层外推到地表震源排列，于是得到单模板的聚焦震源束；
3)将地下目标反射点视为散射点，从散射点出发，模拟正向传播过程，将压力波场由地下逐层外推到地表检波器排列；
4)随后模拟反向聚焦，将压力波场从地表逐层外推到地下成像平面，于是得到单模板的聚焦检波束；
5)将单模板的聚焦震源束和聚焦检波束相乘，得到单模板的分辨率函数；
6)对所有的单模板分辨率函数求和得到双聚焦成像分辨率函数；
7)利用二维Fourier变换，将单模板的聚焦震源束和聚焦检波束从空间域(x，y，z)变换到波束域(kx，ky，z)，随后两者相乘即可得到单模板的AVP函数；
8)对所有的单模板AVP函数求和，得到双聚焦成像AVP函数；
9)沿地下目标层，求取每个反射面元的分辨率函数极大值，即可得到照明强度分布切片。
本发明第二项核心技术是利用双聚焦特征函数来分析优化三维观测系统设计，其基本特点如下：先采用常规的CMP面元属性分析法，进行宏观的全局性观测系统设计；随后针对主要地质目标，采用层状介质双聚焦分析法，进行局部的优化设计。
本发明的三维观测系统优化设计层状介质双聚焦方法包括以下步骤：
1)输入若干个三维观测系统；
2)根据输入的三维观测系统，计算其CMP面元的3个属性参数，即纵横向覆盖次数、炮检距分布和炮检线方位角分布，不同的观测系统得到不同的CMP面元属性参数；
3)分析对比不同CMP面元属性参数对CMP叠加质量的影响，以此优选有利于水平叠加的最佳观测系统；
4)输入层状地质模型及主要地质目标；
5)选择第一个三维观测系统，从地表震源排列出发，先模拟正向传播过程，将压力波场从地表逐层外推到地下目标层的目标反射点；
6)以地下目标点为中心，确定聚焦成像平面，随后模拟反向聚焦过程，将压力波场由成像平面逐层外推到地表震源排列，于是得到聚焦震源束；
7)从目标反射点(散射点)出发，模拟正向传播过程，将压力波场由地下逐层外推到地表检波器排列；
8)随后模拟反向聚焦，将压力波场从地表逐层外推到地下成像平面，于是得到聚焦检波束；
9)由聚焦震源束和聚焦检波束可以直接求得该观测系统对应的分辨率函数、AVP函数和照明强度分布切片；
10)选择第二个、第三个三维观测系统，重复上述第5)-9)步计算过程，不同的观测系统得到不同的双聚焦特征函数；
11)以指定地质目标的最佳偏移成像为评价标准，对计算得到的几组双聚焦特征函数进行分析对比，优选三维观测系统；
12)以双聚焦分析法为主，参考CMP面元分析法，最后确定最佳三维观测系统；
实例1：层状介质双聚焦法验证
三维观测系统选择8线8炮线束型模板；纵向排列形式为0-200-2500m，道距100m，滚动距离200m，滚动20次；横向线间距200m，炮点距100m，滚动距离800m，滚动5次；施工面积3.3km×4.6km，模型大小定为3300×4600×3000m3，目标点指定为模型底中心部位，压力波主频选择30Hz。
设计3个等速模型(参见图5)，各层速度均为3000m/s，模型1为均匀介质，模型2为两层的倾斜介质，模型3为三层的水平介质。预测实验结果三者应该具有相同的分辨率，以此验证算法的正确性。
图6、图7是3个模型的分辨率函数图像，它们几乎完全相同。进一步定量分析表明：纵向x方向的分辨率都为68m，横向y方向的分辨率都为62m，由此证明算法的正确性。
实例2：空间分辨率计算
三维模型大小为4000×4000×3000m3，从上到下三层介质速度依次为2500m/s、2900m/s和3200m/s，中间两层界面倾角分别为10°和5°，倾向相反，走向均为横测线方向，界面中点深度分别为800m和1500m，散射点位于模型底部中间位置，深度为3000m(参见图8)。
观测系统位于地表水平面上，8线8炮，炮点距100m，检波点距50m，检波线距200m，模板横向滚动间距800m，纵向滚动间距200m，滚动次数为横向4次，纵向5次。
成像平面与水平面夹角分别取0°，10°，30°和60°，走向与倾斜层面走向一致。地震波选取30Hz的频率分量。
图9是不同倾角成像面的分辨率函数，从图中可以看出，成像平面与水平面夹角分别为0°和10°时的分辨率函数没有明显变化，而倾角为30°和60°时分辨率函数在倾向方向明显变宽。
图10为不同倾角的空间分辨率垂直切片，左图为成像面倾向方向剖面，右图为成像面走向方向剖面，图中由离散的点，实线，点线和虚线分别表示成像面与水平面夹角为0°、10°、30°和60°时得到的分辨率函数垂直切片。从图中可以看出，四个分辨率函数在成像面走向方向上几乎完全重合，在成像面倾向方向上随倾角增大有明显变化。
进一步的定量分析如下表2：
表2 不同倾角成像面的空间分辨率倾角0°10°30°60°
x-z切片71m72m82m142m
y-z切片109m109m109m109m
图9、图10与表2表明：成像面倾角越大，空间分辨率越差。
本发明的应用，即双聚焦法三维观测系统优化设计流程，见图4。常规的CMP面元属性分析法是宏观的，全局性的，它强调获得最佳的水平叠加效果；而双聚焦特征函数分析法是面向指定勘探目标的，局部性的，强调获得最佳的偏移成像精度。因此，两者恰好互为补充，相辅相成，将两种方法结合在一起，可以极大地完善和提高三维观测系统的设计水平。
下面结合实例3说明双聚焦法三维观测系统优化设计方法。
实例3：宽窄两种三维观测系统分析对比
选择窄方位6线6炮和宽方位12线9炮两种观测系统进行对比。排列总道数均为96道，道距50m，中心放炮，炮点距100m，覆盖次数都是48次，目标点埋深2500m，介质均方根速度2300m/s。
图11是CMP面元属性对比图。对比结果为：6L6S观测系统的柱状图分布比较均匀，纵向信息丰富而完整，有利于保证速度分析的精度；但是缺少横向的信息，不利于三维参数分析。12L9S观测系统的炮检方位线分布比较合理，纵横向信息比较均衡，有利于三维AVO分析；但是它的超级面元柱状图分布不如6L6S。
图12是双聚焦特征函数对比图，两者之间的差别是非常明显的。表3是定量分析结果，分析表明：1)在纵向(y)主测线方向上，6L6S窄方位观测系统的偏移成像分辨率略高于12L9S炮宽方位的，而在横向(x)上两者相同。2)AVP函数图象是在波数平面上，显然12L9S的平面波分布比6L6S丰富，因此有利于三维的物性参数和岩性参数分析；3)12L9S炮的照明强度分布略次于6L6S炮，但仍达到94％以上，因此这两种观测系统都不会产生采集痕迹。
表3 双聚焦特征函数的定量分析 观测系统 6L6S 12L9S
XZ切片 55m 65m
YZ切片 60m 60m
检波采集痕迹 99.3％ 97.7％
震源采集痕迹 99.7％ 96.7％
综合采集痕迹 99.0％ 94.4％
综合上面CMP面元属性分析和双聚焦特征函数分析结果，可以认为这两种观测系统的总体技术性能都是优良的，它们可供不同地质目标的采集设计做不同的选择。
本发明的层状介质双聚焦算法可以选择逐层外推的递归方法，也可选择三维网格化方法。采用大套地层外推方法的主要优点如下：可以用二维网格取代三维网格以减少对内存空间的占有量；以大步长取代小步长可以几倍几十倍的减小运算量，因此本发明所选择的逐层递推方法是一种优选的方案。