中国西部和南方地区因近地表复杂性，采集的地震数据信噪比极低，严重影响地震成像的品质.山体陡峭、沟壑纵横，高差悬殊，采集的原始地震数据走时畸变，导致严重的静校正问题；老地层强烈变形，逆冲推覆出露地表，导致剧烈的近地表速度变化和高陡构造；出露岩层风化剥蚀严重，表层疏松，导致严重的频散效应和高频吸收.如此复杂的近地表地震地质条件给地震资料采集处理带来极大困难，激发接收条件和地震子波 (波形、能量、频率) 一致性极差.总之，对地面地震勘探而言，复杂近地表最大的影响是对所有深层反射波的强散射作用，加上近地表严重的吸收和频散，形成半随机半相干的近地表强散射噪音背景，弥漫整个炮集，淹没深层反射信号.因此，近地表的复杂性是导致地震资料极低信噪比的主要原因，严重制约深部构造的地震成像效果，是复杂地区石油勘探的瓶颈问题，也是国内外石油勘探长期攻关的难题.极低性噪比资料通常具有波形破碎、散化严重、同相轴无法追踪、噪声强度大于有效信号强度的特征.

复杂近地表散射衰减所有的地面观测波场，主要表现为聚焦、散焦和波型转换，与起伏地表的粗糙度和近地表随机介质的相关长度密切相关.如何研究和评价近地表复杂性及其对地震波传播的影响，长期以来一直是地震学的研究热点，大部分研究工作主要集中在复杂近地表地震波场数值模拟 (Bouchon, 1973; Chen, 1990; Carcione, 1994; Hestholm and Ruud, 1994; Komatitsch et al., 1996; Robertsson, 1996; Fu, 2003; 周红和陈晓非, 2007; 管西竹等, 2011)，大多采用大尺度不规则边界结构+内部随机介质的速度模型，侧重模拟算法研究，波场复杂性分析较少.采用地震全波形数值模拟方法来研究近地表散射机制，并不能直观了解散射与地表粗糙度和近地表随机介质参数的相关性，而且涉及海量模拟计算.大量的解析/半解析方法应运而生 (Gilbert and Knopoff, 1960; Hudson et al., 1973; Snieder, 1986; Fu et al., 2002; Fu, 2005; Hu et al., 2009; Yu and Fu, 2012; Rao and Wang, 2015)，由于逼近精度的局限性，这些方法不能完全适用于复杂近地表散射表征.总之，近地表强散射对勘探地震而言是噪音，如何量化地震炮集的散射噪音强度一直是石油勘探未解决的问题，这与近地表速度横向变化和结构倾角分布密切相关.

由于不规则起伏地表和内部非均质极强的随机性，难以找到一套行之有效的近地表复杂性量化表征方法.Ashford等 (Ashford et al., 1997) 较为系统地给出了地形倾角和相对高程与地震波类型、波长、入射角之间的相互关系，Chou (Chou et al., 1999) 以起伏地表数据特征为基础来定义地形复杂度.基于复杂介质地质非均质谱 (速度横向变化和地层倾角变化非均质谱) 与地震传播算子角谱的相互作用，Fu (Fu, 2010) 提出了一种复杂介质复杂度概念及其计算方法，并将其应用于库车坳陷高陡构造的复杂性定量分析 (董文等, 2011; 符力耘等, 2013) 和复杂海底地貌特征描述 (李绪宣等, 2011).

基于复杂近地表边界元法波动积分方程数值模拟研究，本文提出一种复杂近地表散射振幅矩阵方法来分析近地表散射强度.首先对复杂近地表结构进行边界元配置方法离散，根据边界积分方程生成矩阵方程HU=GQ，其中H和G是散射振幅系数矩阵，U和Q是位移向量及其法向导数.我们不求解该矩阵方程 (涉及海量计算)，只是利用矩阵分析技术来解析散射振幅系数矩阵，研究复杂近地表结构对不同频率波场的散射强度.H矩阵的求取是整个方法的关键，它不仅需要计算边界上任意两点之间的基本解，同时还需求取该基本解在边界上的法向导数.事实上基本解包含了任意两个单元的距离信息，而基本解的法向导数则反映了这两个单元的相对位置关系，H矩阵的这种特殊性为我们提供了评价不同地质结构复杂度的可能性.通过积分算子计算任意两单元之间的格林函数，并进行高斯数值积分计算得到H矩阵，对H矩阵进行矩阵分析来研究近地表散射特征.作为初步评价手段，我们采用矩阵元素总和与矩阵维数之比作为表征散射振幅系数矩阵散射特征的标量复杂系数.通过理论和实际模型测试，形成了一套行之有效、计算快速的近地表复杂性分析方法.

式 (1) 是包含源点的边界波动积分方程，它将区域内的微分方程变换成边界上的积分方程，描述了源点扰动后所形成的综合波场，是地震波在均匀介质中传播的控制方程，

(1)

其中等号左边两项分别表示边界散射场和源波场，即总波场u(r) 由这两部分共同决定.其过程可由图 1所示模型描述，均匀介质Ω由不规则边界Γ和内部区域Ω两部分组成，r、r′、r₀为Ω上的任意点，震源S (ω)位于r₀处，ω=2πf为角频率，f为地震波频率，G(r, r′)=iH₀⁽¹⁾(k₀|r-r′|)/4为二维背景介质中的格林函数，H₀⁽¹⁾为汉克尔函数，k₀=ω/v₀为背景波数，v₀为背景波速，C(r)=θ(r)/2π为边界点r处与边界几何形状有关的系数，由r点处的边界切线面向Ω的夹角θ(r) 决定，对于光滑边界，C(r)=0.5.

图 1

Fig. 1

图 1 非规则边界下均匀介质模型 Fig. 1 A homogeneous medium with a free boundary

式 (1) 在整个计算区域Ω上成立，因此可以将r和r′点移至边界上，并将等式左边第一项拆分后移至等式右边，则式 (1) 转化为

(2)

之所以拆分第一项，是为了将边界上的函数u(r′) 及其法向导数∂u(r′)/∂n分离，便于单独处理.为叙述简便，令 (2) 式等号左边等于M，则有

(3)

将边界用常单元离散，即节点取在边界单元的中点，且每个边界单元上任意点的函数值都等于节点处的函数值，则式 (3) 可离散为 (对于常单元C(r)=0.5)

(4)

式 (4) 中u(r′_j) 是边界单元Γ_j上任意点r′_j处的函数值，i, j=1, 2, …, s为边界单元序号，s为边界单元总数.由于是常单元剖分，因此有u(r′_j)=u(r_j)，则式 (4) 可进一步化为

(5)

令 (5) 式右边的积分项等于，即

(6)

则式 (5) 化简为

(7)

式 (7) 表示对于当前计算节点i，需将其与所有边界单元 (包括节点i所在单元本身) 进行一次积分，对于有s个节点的边界，共需计算s次积分，i点结束后依次移向下一个节点，直至所有节点计算完成，总积分次数为s².注意到式 (7) 中，在j循环过程中会出现一次i=j的情况，此时有u(r_i)=u(r_j)，可以将其合并为.重新定义

(8)

其中为Kronecker函数，则式 (7) 可进一步化简为

(9)

将式 (9) 应用于所有边界单元，则可得到s个方程，写成矩阵形式为

(10)

其中等式左边系数剧矩阵便是本文研究的重点对象H矩阵.

由式 (8) 和 (6) 可知

(11)

将格林函数G(r_i, r′_j)=iH₀⁽¹⁾(k₀l_ij)/4代入式 (11) 并化简得

(12)

式中l_ij=|r_i-r_j′|为节点r_i到任意点r′_j的距离，H_-1⁽¹⁾(k₀l_ij) 为H₀⁽¹⁾的一阶导数 (根据汉克尔函数的导数=h_ij/l_ij为方向导数，h_ij是节点r_i到Γ_j的垂直距离，如图 2所示.

图 2

Fig. 2

图 2 H矩阵元素计算示意图 Fig. 2 Schematic diagram of element calculation in matrix H

式 (12) 的计算需分两种情况分别进行讨论：

①当i=j时，节点r_i和任意点r_j′在同一边界单元上，此时l_ii与n相互垂直，即，同时δ_ii=1，因此对于式 (12) 有H_ii=0.5，即H矩阵的主对角线元素都为0.5.

②当i≠j时，节点r_i和任意点r′_j在不同边界单元上，此时Γ_j上所有的点都参与积分运算，且有δ_ij=0.实际计算时，我们采用高精度的高斯积分进行，为此需要将边界进行无因次转换，转换关系由式 (13) 决定：

(13)

式中ξ∈[-1, 1]是无因次坐标，L_j=为边界单元Γ_j的实际长度，(x₁, y₁) 和 (x₂, y₂) 是边界单元Γ_j两个端点的实际坐标.将式 (13) 代入式 (11)，并离散成高斯积分有

(14)

其中K、ξ_k、w_k分别为高斯积分点数、积分点位置和积分点权系数，l_ik为边界单元Γ_i上节点r_i与Γ_j上高斯积分点ξ_k之间的距离.通过式 (14) 可以计算H中所有元素.

以任意两点之间的基本解及其在边界上的法向导数经过高斯数值积分计算得到的H矩阵，描述了所有边界上不同点之间产生的散射影响，同时包含了边界几何形态信息，为近地表构造复杂性分析提供了可能性.H矩阵是一个复数矩阵，其维度表示边界单元总数，矩阵中每个元素可以认为是从源节点传播到场节点的波场在场节点边界外法向上的分量，与节点之间的距离、相对位置和节点所在边界的几何形态密切相关.为了便于在实数域说明问题，我们对H矩阵每个元素求模，即将矩阵H转化为振幅矩阵H.元素值越大说明两个节点之间的散射影响越强，整个模型的复杂性则由所有边界单元之间的相互关系共同决定.

为了深入分析H矩阵的特性，我们设计了图 3上所示的三个2500 m×60 m水平层状模型.将每层边界划分成50个完全相同的常单元，并按从左到右的顺序依次编号，图中黑色小圆点和线段分别表示节点和单元端点，三个模型对应的H矩阵如图 3下所示.矩阵主对角线元素表示源节点和场积分点重合，根据前面推导，在水平层常单元时为常数0.5.主对角线元素以外区域表示源节点和场积分点不重合，此时若源节点和场节点在同一层边界上，由于任意点的外法向与所有边界单元相互垂直，根据积分式 (11)，元素值为0.以图 3b对应的矩阵H为例，我们用虚线将矩阵分为4个区域，分别用A、B、C、D标注.可见，区域A和D为边界1和边界2的自相互作用区，而B和C为边界1和边界2之间的相互作用区.从B和C区元素值大小来看，离主对角线位置越远，元素值越小，表示两个边界节点之间的相互影响越弱.例如图 3f中的C和G区对应距离较远的边界1和3之间的相互影响，浅灰色斜线反映元素值较小，各元素对应模型由图 3c中节点1和101、2和102、…、50和150等组合所形成，元素值大小与距离成反比.上述各层节点之间相互作用的清晰分布仅对应于水平层状模型，当层边界不规则时，H矩阵元素值分布变得不规则，但上述分析方法依然适用于所有情形.

图 3

Fig. 3

图 3 三个水平层模型 (上) 及其对应的H矩阵 (下) Fig. 3 Three horizontal-layered model and the corresponding matrices

根据上述H矩阵元素的分布特征，可以定义一个标量系数来直接分析近地表构造的复杂性.一般的矩阵特征分析多采用矩阵范数、条件数、特征值等，经过仔细研究，这些常规的H矩阵参数无法满足本文的研究目的.以矩阵1范数为例，它仅突出H矩阵中具有最大叠加影响的某个边界单元特征，不能反映H矩阵分布的复杂性.此外，反映图像纹理结构特征的二阶矩、对比度、相关、熵等特征量比矩阵分析参数要好一些，但依然不足以表征H矩阵来分析近地表构造的复杂性.考虑到H矩阵元素反映了不同节点之间相互关系，以及节点所在边界几何形态的影响，最简单的综合评价方法就是H矩阵所有元素求和，所有元素求和仅作为初步评价手段以突出矩阵整体特性，基本反映了H矩阵与近地表结构复杂性的正相关关系，但忽略了不同元素的作用权重及其之间的非线性叠加影响.H矩阵标量系数定义如下：

(15)

其中s为矩阵维数，用以消除边界单元划分数量对结果产生的影响.

如图 4-7所示，我们设计一系列从简单到复杂的模型来验证上述方法的有效性.每一组模型仅包含一种近地表结构复杂性，包括山峰个数、高度、坡度、地表延续度等.各例中的山峰均由不同参数控制的高斯曲线生成，背景速度取为2500 m·s^-1，子波主频为30 Hz，各速度模型对应的H矩阵用灰度表示元素值大小，计算的标量系数如图 10所示.

图 4

Fig. 4

图 4 不同高度山峰地貌模型 (上) 及其H矩阵 (下) Fig. 4 Models of different peak heights (up) and the corresponding matrices (down)

图 5

Fig. 5

图 5 不同坡度山峰地貌模型 (上) 及其H矩阵 (下) Fig. 5 Models of different slopes (up) and the corresponding matrices (down)

图 6

Fig. 6

图 6 不同跨度山峰地貌模型 (上) 及其H矩阵 (下) Fig. 6 Models of different spans (up) and the corresponding matrices (down)

图 7

Fig. 7

图 7 不同个数山峰地貌模型 (上) 及其H矩阵 (下) Fig. 7 Models of different peaks (up) and the corresponding matrices (down)

图 8

Fig. 8

图 8 单层到多层的起伏界面近地表模型 (上) 及其H矩阵 (下) Fig. 8 Models of different interfaces and layers (up) and the corresponding matrices (down)

图 9

Fig. 9

图 9 不同边界剖分模型 (上) 及其H矩阵 (下) Fig. 9 Models with different boundary discretization (up) and the corresponding matrices (down)

图 10

Fig. 10

图 10 图 4-9不同近地表模型对应的标量复杂系数 Fig. 10 Complexity coefficients of models from Fig. 4 to Fig. 9

图 4代表不同高度的山峰地貌模型，山峰高度分别为160 m、240 m、320 m和400 m，模型水平跨度均为2000 m，边界单元数均为200，其H矩阵递增的灰度值和图 10中曲线1递增的标量系数反映了复杂性随相对高差增加而增加.图 5中四个模型代表不同坡度的山峰地貌，其标量系数 (图 10中曲线2) 的递增表明坡度的增加对近地表复杂性的影响.图 6表示不同水平跨度地貌模型，分别为2000 m、1500 m、1000 m和500 m，其标量系数 (图 10中曲线3) 的递增表明跨度越小复杂性越强，即当山峰不变而水平部分无限延伸时，山峰的影响将逐渐减弱.上述三种类型的地貌其复杂性增加对应的标量系数递增斜率几乎相同，表明这些复杂因素的影响属于同一个量级.此外，图 10中曲线1到3的初值逐渐升高，说明扁平山峰地貌复杂性弱于窄而高的山峰地貌.图 7研究山峰个数增加的复杂性，四个模型中每个山峰的高度、跨度和边界单元数分别为360 m、2000 m和200，其H矩阵剧烈变化的灰度值及其标量系数 (图 10中曲线4) 的剧增反映了多峰起伏地表具有很强的复杂性，这与实际地面地震观测中多峰起伏地表很强的散射噪音相符.

图 8为单层到多层近地表模型，地层分界面随机起伏，层内为均匀介质，随着层数的增加，模型越来越复杂.对应的H矩阵灰度值剧烈变化及其标量系数 (图 10中曲线5) 的剧增表明近地表不规则多层结构具有更强的复杂性.其中，形态差不多的多层 (图 8前面三个模型) 其标量系数呈线性增加趋势，而形态迥异的多层 (图 8第四个模型) 其标量系数呈非线性增加，尤其是成像复杂的陡倾角结构加剧了近地表复杂性.一种有效的评价方法应该是稳定的，对相同或相似模型的评价应该获得相近的结果，不受离散方式的影响.图 9所示的四个模型与图 8d模型完全相同，区别在于边界剖分精细程度逐渐递增，剖分的边界单元数差别较大，分别为80、150、230和535，测试边界剖分对复杂系数计算的影响.计算的标量系数如图 10中近似水平的曲线6所示，复杂系数基本维持在0.073左右.模型离散要求边界单元与实际边界之间的误差尽量小，在起伏变化剧烈的地方用小而多的边界单元，而起伏平缓的地方则用长而少的边界单元，确保离散模型基本反映实际模型的主要特征.

图 11用于测试方法对不同频率的响应，四个模型完全相同，均取自于图 9中模型a的第一和第三层地质分界面，频率参数分别为10 Hz、35 Hz、50 Hz和65 Hz.相同的地质模型对于不同频率的地震波表现出不同的复杂度，频率越高，波长越小，地震波越容易受构造变化的影响，模型复杂度相对也就越大.图 11e-11h所示各频率H矩阵的非零元素分布模式完全相同，这是由于采用同一个地质模型的结果，而元素值的差异则体现了不同频率产生的影响，从图中可以看出，灰度值逐渐变大，整体复杂度呈明显增加趋势.图 12给出了不同频率H矩阵的标量复杂系数，可见，随着频率增加，近地表模型的标量复杂系数大幅度上升，而且上升曲线的斜率与频率间隔也表现出一定正相关性.

图 11

Fig. 11

图 11 相同模型对不同频率地震波的影响 (上) 及其H矩阵 (下) Fig. 11 The influence of the same model on seismic waves of different frequencies (up) and the matrices (down)

图 12

Fig. 12

图 12 图 11中不同频率对应标量复杂系数曲线 Fig. 12 Scalar complex coefficient curve of the four different frequencies in Fig. 11

图 13a是我国某地区实际近地表速度结构，其复杂性主要表现为强起伏地形和强横向变化速度.从图中可以看出，两个大规模山丘控制着整个地区的主要地形特征，隆起部位伴随着参差不齐的锯齿状波动，特别是在横坐标1600~2200 m和2900~3400 m范围内形成陡坡.近地表低速层较薄，成层性和连续性差.总体而言，该地区对地震勘探是一个极大的挑战.

图 13

Fig. 13

图 13 实际地表起伏模型 (上) 及其地表平滑模型 (中) 和对应的矩阵 (下，左对应原模型，右对应平滑模型) Fig. 13 A real model with undulating surface (up) and its duplicate that has the surface smoothed (middle) and the corresponding matrices (down)

我们用本文方法来分析该模型的地质复杂性.首先根据速度分布勾勒出主要地层分界面，再利用边界单元对地表线和分界面进行剖分，边界和边界单元的标号遵循从左往右和从上到下的原则.剖分结果如图 13a所示，共形成14条边界和445个边界单元.为了对比分析，增加结果的可信度，我们对原模型地形线进行平滑，减小地形起伏对整个模型复杂性的影响，平滑模型如图 13b所示，用同样的方法进行边界单元剖分，边界数量仍为14，边界单元数随着地形线的去尖锐化而减为265.两个模型对应的H矩阵分别如图 13c和13d所示.最后利用式 (16) 选择频率10 Hz、35 Hz、50 Hz和65 Hz分别计算两者标量复杂系数，结果如图 14所示，曲线1和2分别对应模型 (a) 和 (b).从图中可以看出，曲线1完全在曲线2的上方，由此说明在频率参数相同的情况下，本文得出的标量复杂系数确实比表面平滑后的模型要复杂.随着频率的增大，两个模型的复杂度同时增加，这一点同样与理论相吻合.实际上本例中，地表曲线的平滑对整个模型的整体影响并非很大，但我们的方法从数值上准确反映了两者的差异，表现出对结构变化良好的响应.

图 14

Fig. 14

图 14 不同频率下图 13中模型a和b对应复杂度系数曲线 Fig. 14 The corresponding complexity coefficient curve of the models a and b in Fig. 13 at different frequencies

为了增加上述结果的可信度，我们进一步从地震波场的角度来分析两个模型之间的差异.边界元法波动积分方程地震模拟采用的水平和深度方向网格间距为5 m，时间采样间隔为0.001 s，采样时间为1.2 s，Ricker子波主频为30 Hz，震源位于地表水平位置2500 m处.图 15a和15b分别为两个模型在0.5 s时刻的波场快照 (截取1000~4000 m之间的部分).可见，图 15a中S区域弥漫着各种相互干涉的反射和散射波场，可辨识度极低，这是因为ac地段地形起伏变化剧烈，形成极强的局部散射波分布，其中振幅尤为强烈的散射波前l1，是由位置较低的震源波场向上传播至c点附近时被散射回来的波.各个尖锐的山峰和山谷产生强散射噪音，特别是在陡峰b点附近，波场逐渐破碎，能量被限制在坡峰处直至耗散殆尽，该处的检波器接收到的几乎是强散射噪声.相比于图 15a，图 15b中波场相对简单，反射波源易于追踪，由于平滑作用削弱了地形的起伏变化，使散射波强度锐减，特别是陡峰b处变得相对平缓，发生局部振幅放大，形成连续性较好的明显波前l2，尽管也属于噪声，但可识别度高，相对于随机噪音更易消除，且不会完全淹没此处检波器接收的深部反射有效信号.图 15c和15d是对应的单炮记录，可见图 15c中的初至同相轴连续性差，走时起伏较大，波场较为破碎，近地表反射和散射占据整个剖面，特别是几处山峰强散射形成的信号杂乱区域E1-E4，各种波相互干涉，形成局部强散射噪声环境.而图 15d中的同相轴更为连续规则，波场可辨识度高，虽然剖面也受大量地形噪声影响，但整体显得更为规则和简单，几处山峰强散射形成的信号杂乱区域或消失或减弱，对有效波场的干扰范围明显减小，去噪更容易.

图 15

Fig. 15

图 15 原模型及其平滑模型地震模拟0.5 s时刻波场快照 (上) 及单炮记录 (下) Fig. 15 Snapshots at time 0.5 s (up) and the single shot gathers of the two models in Fig. 13 (down)

复杂近地表地区地震勘探实践表明，半随机半相干的近地表强散射噪音弥漫整个炮集，淹没了深层反射信号，是导致地震资料极低信噪比的主要原因.如何研究和评价近地表散射强度一直是石油勘探未解决的问题，与起伏地表的粗糙度、近地表速度横向变化和结构倾角分布密切相关.本文基于前期复杂近地表边界元法波动方程数值模拟研究，提出一种复杂近地表散射振幅矩阵方法来分析近地表散射强度，对于根据边界积分方程生成的矩阵方程，该方法不求解矩阵方程 (涉及海量计算)，只是利用矩阵分析技术来解析矩阵方程中的散射振幅系数矩阵，研究复杂近地表结构的地震响应复杂性及其对不同频率波场的散射强度.方法的优势在于利用边界元对近地表结构几何特征的精确描述，由基本解及其在边界上的法向导数经过高斯数值积分计算得到的散射振幅系数矩阵，表征了复杂近地表结构对地震波的散射特征.作为初步评价手段，我们采用矩阵元素总和与矩阵维数之比作为表征散射振幅矩阵散射特征的标量复杂系数.

通过理论模型试验表明，在起伏地貌中，山峰的高度、坡度和跨度对散射强度的影响属于同一个量级，而山峰个数的增加会极大提高散射强度，反映了多峰起伏地表具有很强的复杂性.引起标量复杂系数剧增的近地表模型是多层不规则界面近地表模型，特别是形态迥异、成像复杂的陡倾角结构加剧了近地表的复杂性，使标量系数呈非线性增加趋势.随着计算频率的增大，标量复杂系数大幅度上升，复杂近地表对高频波表现出更强的散射特征.实际近地表速度模型试验表明，在近地表结构中，界面几何结构和速度的局部剧烈变化会引起标量复杂系数剧增，特别是地形的局部剧烈变化导致极强散射波干涉，形成极其复杂的局部散射噪音分布区域，有效波可辨识度极低.总之，本文提出的近地表复杂性分析方法能够准确反映不同模型之间的结构复杂性差异，对结构变化反应灵敏，可作为地震勘探近地表复杂性的简单快速的评价方法.此外，边界剖分的精细程度对结果没有影响，只要边界单元尽可能逼近实际边界即可.

本文提出的方法还存在许多不足之处，例如标量复杂系数的物理含义还没有阐述清楚，其大小分布与起伏地表的粗糙度 (均方根高和相关长度) 和近地表随机介质的相关长度之间的直接关系是下一步将要开展的研究工作.