一种多层近地表地震地质复杂性的定量分析方法

摘要

本发明提供了一种多层近地表地震地质复杂性定量分析的方法，该方法利用边界元对近地表结构几何特征的精确表征，研究起伏地表和非规则地质分界面对地震波传播的影响，由基本解在边界上的外法向导数经过高斯数值积分表征边界单元对节点处的波场影响，不仅描述了任意两点之间的相互影响，同时还刻画了边界形状特征的影响，能够较为准确地对多层近地表地震地质复杂性进行定量分析。

权利要求

1.一种多层近地表地震地质复杂性定量分析的方法，其特征在于，该方法包括以下步骤：(1)获取地形线和地质分界线，对获得的地形线及地质分界线进行离散，得到边界单元，确定各边界单元端点的空间位置坐标，并确定每个边界单元上相应节点的空间位置坐标；

(2)任意选择一个步骤(1)确定的节点作为当前节点，任意选择一个步骤(1)确定的边界单元作为当前边界单元，用格林函数在当前边界单元上的外法向导数表征当前节点与当前边界单元上任意点之间的相互影响，所述格林函数的关系式如下式(5)所示，其中，为0阶第一类汉克尔函数，

k0＝ω/v0为背景波数，

ω为圆频率，v0为背景波速，

r为当前节点位置，

r'为边界单元上任意点位置，

i和j表示边界单元序号，即i,j＝1,2,…,s；

(3)用积分表征当前边界单元对当前节点处波场的影响，所述积分的表达式如下式(6)所示，其中，

Γj为序号为j的边界单元，

n是Γj的外法向；

(4)重复步骤(2)和步骤(3)，直至用积分表征了所有节点与所有边界单元之间的相互影响，形成第一矩阵H；

(5)以步骤(4)中得到的第一矩阵H为基础，形成相应的散射振幅系数矩阵 并将散射振幅系数矩阵 的元素总和与矩阵维数的比值表示为多层近地表地震地质复杂性系数，实现多层近地表地震地质复杂性定量分析。

2.根据权利要求1所述的方法，其特征在于，在步骤(1)中，对获得的地形线和地质分界线用常单元进行离散。

3.根据权利要求1所述的方法，其特征在于，步骤(4)中，所述的第一矩阵H为，其中，矩阵的每个元素Hij由下式(13)所示高斯数值积分得到，其中，

lij为节点ri到任意点rj'的距离，为 的一阶导数；

x1，y1和x2，y2是边界单元Γj两个端点的实际坐标；

K表示高斯积分点数；

ξk表示积分点位置；

wk表示积分点权系数；

lik表示边界单元Γi上节点ri与边界单元Γj上高斯积分点ξk之间的距离。

4.根据权利要求1所述的方法，其特征在于，步骤(5)中，所述 矩阵如下式(14)所示：

5.根据权利要求1所述的方法，其特征在于，步骤(5)中，多层近地表地震地质复杂性系数如下式(15)所示，其中，

comp表示多层近地表地震地质复杂性系数；

i和j表示边界单元序号，即i,j＝1,2,…,s；

|Hij|表示散射振幅系数矩阵 的元素；

s表示散射振幅系数矩阵 的维数。

6.根据权利要求1所述的方法，其特征在于，步骤(5)得到的多层近地表地震地质复杂性系数的数值与多层近地表地震地质复杂性呈正相关关系。

7.根据权利要求1所述的方法，其特征在于，还包括以下步骤：依据步骤(5)中得到的多层近地表地震地质复杂性系数，来选择期望的地震成像方法、地震波数据处理方法、地表起伏复杂地区的数据采集方法。

说明书

一种多层近地表地震地质复杂性的定量分析方法

技术领域

[0001] 本发明涉及地球物理勘探技术领域，特别涉及一种地震勘探采集处理中的地震资料质量分析技术。

背景技术

[0002] 地表复杂性对地震成像方法的选用、地震波数据的处理和地表起伏复杂地的数据采集均具有指导意义，然而，目前地表复杂性大多是定性评估，难以用定量的方法进行评估。

[0003] 中国西部和南方地区近地表复杂性大，采集的地震数据信噪比极低，严重影响地震成像的品质，具体表现为：山体陡峭、沟壑纵横，高差悬殊，采集的原始地震数据走时畸变，导致严重的静校正问题；老地层强烈变形，逆冲推覆出露地表，导致剧烈的近地表速度变化和高陡构；出露岩层风化剥蚀严重，表层疏松，导致严重的频散效应和高频吸收。

[0004] 如此复杂的近地表地震地质条件给地震资料采集处理带来了极大困难，导致激发接收条件和地震子波(波形、能量、频率)一致性极差。

[0005] 对地面地震勘探而言，复杂近地表最大的影响是对所有深层反射波的强散射作用，加上近地表严重的吸收和频散，形成半随机半相干的近地表强散射噪音背景，弥漫整个炮集，淹没了深层反射信号，导致有效波波形破碎，波散化严重，几乎看不到同相性较好的有效波组。

[0006] 因此，近地表的复杂性是导致地震资料极低信噪比的主要原因，严重制约深部构造的地震成像效果，是复杂地区石油勘探的瓶颈问题，也是国内外石油勘探长期攻关的难题。

[0007] 复杂近地表散射衰减所有的地面观测波场，主要表现为聚焦、散焦和波型转换，与起伏地表的粗糙度和近地表随机介质的相关长度密切相关。如何研究分析和评价近地表复杂性及其对地震波传播的影响，长期以来一直是地震学的研究热点，现有技术的大部分研究工作主要集中在复杂近地表地震波场数值模拟，大多采用大尺度不规则边界结构+内部随机介质的速度模型，侧重模拟算法研究，而对波场复杂性的分析较少。

[0008] 现有技术中存在采用地震全波形数值模拟方法来研究近地表散射机制，但其并不能直观了解散射与地表粗糙度和近地表随机介质参数的相关性，而且涉及海量模拟计算，虽然大量的解析/半解析方法应运而生，但由于逼近精度的局限性，这些方法不能完全适用于复杂近地表散射表征。

[0009] 总之，近地表强散射对勘探地震而言是噪音，如何量化地震炮集的散射噪音强度，一直是地面地震勘探，即，石油勘探未解决的问题，这与近地表速度横向变化和结构倾角分布密切相关。

[0010] 由于不规则起伏地表和内部非均质极强的随机性，目前难以找到一套行之有效的近地表复杂性量化表征方法，如，Ashford等较为系统地给出了地形倾角和相对高程与地震波类型、波长、入射角之间的相互关系；Chou以起伏地表数据特征为基础来定义地形复杂度；中国专利CN200810105063公开了一种地下非均匀介质地震探测复杂性定量分析方法，其基于复杂介质地质非均质谱(速度横向变化和地层倾角变化非均质谱)与地震传播算子角谱的相互作用；Fu提出了一种复杂介质复杂度概念及其计算方法，并将其应用于库车坳陷高陡构造的复杂性定量分析和复杂海底地貌特征描述，然而，该定量分析方法并不能适用于地表复杂性的定量分析。

[0011] 在现有技术中，对于地表复杂性的定量分析方法尚未建立起来，因此，需要开发一种对于多层近地表地震地质复杂性的定量方法，并基于此可选择适当的地震成像方法、地震波数据的处理方法、和地表起伏复杂地区的数据采集方法等，从而指导地质勘探。

发明内容

[0012] 为了解决上述问题，本发明人进行了锐意研究，结果发现：以地形线和地质分界线为起点，以边界波动积分方程为基础，用格林函数表征波场的复杂性，再结合波场的复杂性，用高斯数值积分表征不同产状的边界单元对不同节点处波场的影响，形成矩阵，然后以所形成的矩阵为基础得到散射振幅矩阵，将散射振幅矩阵元素的总和与矩阵维数之比的比值作为表征多层近地表地震地质复杂性的系数，用其来定量分析多层近地表地震地质复杂性，并基于此可选择适当的地震成像方法、地震波数据的处理方法和地表起伏复杂地区的数据采集方法等，从而指导地质勘探，从而完成了本发明。

[0013] 本发明的目的在于提供一种用多层近地表散射振幅矩阵来定量分析多层近地表地震地质复杂性的方法，该方法包括以下步骤：

[0014] (1)获取地形线和地质分界线，对获得的地形线及地质分界线进行离散，得到边界单元，确定各边界单元端点的空间位置坐标，并确定每个边界单元上相应节点的空间位置坐标；

[0015] (2)任意选择一个步骤1确定的节点作为当前节点，任意选择一个步骤1确定的边界单元作为当前边界单元，用格林函数在当前边界单元上的外法向导数表征当前节点与当前边界单元上任意点之间的相互影响；

[0016] (3)用积分表征当前整个边界单元对当前节点处波场的影响；

[0017] (4)重复步骤2和步骤3，直至用积分表征了所有节点与所有边界单元之间的相互影响，形成第一矩阵H；

[0018] (5)以步骤4中得到的第一矩阵H为基础，形成相应的散射振幅系数矩阵 并将散射振幅系数矩阵 的元素总和与矩阵维数的比值表示为多层近地表地震地质复杂性系数，实现多层近地表地震地质复杂性定量分析。

附图说明

[0019] 图1示出本发明多层近地表地震地质复杂性定量分析方法流程框图；

[0020] 图2示出包含源点的边界波动积分方程模型示意图；

[0021] 图3示出山峰高度为160m的模型；

[0022] 图4示出山峰高度为240m的模型；

[0023] 图5示出山峰高度为320m的模型；

[0024] 图6示出山峰高度为400m的模型；

[0025] 图7示出用发明方法得到的图3示出模型相应的 矩阵；

[0026] 图8示出用发明方法得到的图4示出模型相应的 矩阵；

[0027] 图9示出用发明方法得到的图5示出模型相应的 矩阵；

[0028] 图10示出用发明方法得到的图6示出模型相应的 矩阵；

[0029] 图11示出山峰坡度为25°的模型；

[0030] 图12示出山峰坡度为35°的模型；

[0031] 图13示出山峰坡度为45°的模型；

[0032] 图14示出山峰坡度为65°的模型；

[0033] 图15示出用发明方法得到的图11示出模型相应的 矩阵；

[0034] 图16示出用发明方法得到的图12示出模型相应的 矩阵；

[0035] 图17示出用发明方法得到的图13示出模型相应的 矩阵；

[0036] 图18示出用发明方法得到的图14示出模型相应的 矩阵；

[0037] 图19示出水平跨度为2000m的模型；

[0038] 图20示出水平跨度为1500m的模型；

[0039] 图21示出水平跨度为1000m的模型；

[0040] 图22示出水平跨度为500m的模型；

[0041] 图23示出用发明方法得到的图19示出模型相应的 矩阵；

[0042] 图24示出用发明方法得到的图20示出模型相应的 矩阵；

[0043] 图25示出用发明方法得到的图21示出模型相应的 矩阵；

[0044] 图26示出用发明方法得到的图22示出模型相应的 矩阵；

[0045] 图27示出山峰数为1个的模型；

[0046] 图28示出山峰数为2个的模型；

[0047] 图29示出山峰数为3个的模型；

[0048] 图30示出山峰数为4个的模型；

[0049] 图31示出用发明方法得到的图27示出模型相应的 矩阵；

[0050] 图32示出用发明方法得到的图28示出模型相应的 矩阵；

[0051] 图33示出用发明方法得到的图29示出模型相应的 矩阵；

[0052] 图34示出用发明方法得到的图30示出模型相应的 矩阵；

[0053] 图35示出地层数为1层的模型；

[0054] 图36示出地层数为2层的模型；

[0055] 图37示出地层数为3层的模型；

[0056] 图38示出地层数为4层的模型；

[0057] 图39示出用发明方法得到的图35示出模型相应的 矩阵；

[0058] 图40示出用发明方法得到的图36示出模型相应的 矩阵；

[0059] 图41示出用发明方法得到的图37示出模型相应的 矩阵；

[0060] 图42示出用发明方法得到的图38示出模型相应的 矩阵；

[0061] 图43示出地震波频率为10Hz的模型；

[0062] 图44示出地震波频率为35Hz的模型；

[0063] 图45示出地震波频率为50Hz的模型；

[0064] 图46示出地震波频率为65Hz的模型；

[0065] 图47示出用发明方法得到的图43示出模型相应的 矩阵；

[0066] 图48示出用发明方法得到的图44示出模型相应的 矩阵；

[0067] 图49示出用发明方法得到的图45示出模型相应的 矩阵；

[0068] 图50示出用发明方法得到的图46示出模型相应的 矩阵；

[0069] 图51示出剖分单元数为80的模型；

[0070] 图52示出剖分单元数为150的模型；

[0071] 图53示出剖分单元数为230的模型；

[0072] 图54示出剖分单元数为535的模型；

[0073] 图55示出用发明方法得到的图51示出模型相应的 矩阵；

[0074] 图56示出用发明方法得到的图52示出模型相应的 矩阵；

[0075] 图57示出用发明方法得到的图53示出模型相应的 矩阵；

[0076] 图58示出用发明方法得到的图54示出模型相应的 矩阵；

[0077] 图59示出实际地下速度结构模型；

[0078] 图60示出实际地下速度结构模型经过地表平滑后的模型；

[0079] 图61示出用发明方法得到的图60示出模型相应的 矩阵；

[0080] 图62示出用发明方法得到的图61示出模型相应的 矩阵；

[0081] 图63示出用本发明方法得到的图59和图60示出模型的复杂性系数曲线，其中，[0082] 曲线1对应图59示出的模型；

[0083] 曲线2对应图60示出的模型；

[0084] 图64示出实验例1～实验例5和实验例7得到的复杂性系数曲线；

[0085] 图65示出实验例6得到的复杂性系数曲线。

具体实施方式

[0086] 下面通过对本发明进行详细说明，本发明的特点和优点将随着这些说明而变得更为清楚、明确。

[0087] 在这里专用的词“示例性”意为“用作例子、实施例或说明性”。这里作为“示例性”所说明的任何实施例不必解释为优于或好于其它实施例。

[0088] 以下详述本发明。

[0089] 根据本发明提供的多层近地表地震地质复杂性的定量分析方法，如图1所示，包括以下步骤：

[0090] 步骤1，获取地形线和地质分界线，对获得的地形线及地质分界线进行离散，得到边界单元，确定各边界单元端点的空间位置坐标，并确定每个边界单元上相应节点的空间位置坐标。

[0091] 在本发明中，所采集的信息不仅包括了地表地形信息，同时也包括地下地质分界面信息，即，本发明提供的分析方法考虑全面，不仅考虑了地表地形对多层近地表地震地质复杂性的影响，也包含了地下地质分界面对全局的影响，分析方法较为全面、客观、科学。

[0092] 在本发明中，对获取地形线和地质分界面信息的方法不做特别限定，可以使用现有技术中任意一种获取地形线和地质分界面信息的方法，如地质推断法、测井法和/或地震探测法等。

[0093] 优选地，对地形线和地质分界线进行编号，便于计算机的识别。在本发明中，对地形线和地质分界线的编号顺序不做特别限定，只需要将每条地形线和地质分界线赋予唯一的编号即可。

[0094] 在本发明中，对获得的地形线和地质分界线用常单元进行离散，即节点取在边界单元的中点，并且每个边界单元上任意点的函数值都等于节点处的函数值，从而得到s个相互独立的边界单元，并确定了上述边界单元端点的空间位置坐标，进而确定了每个边界单元上相应地节点的空间位置坐标。

[0095] 本发明以包含源点的边界波动积分方程为基础，定义了多层近地表地震地质复杂性系数，来评价和分析多层近地表地震地质复杂性。

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

[0097]

[0098] 其中，

[0099] 等号左边第一项和第二项依次表示边界散射场和源波场，即总波场u(r)由这两部分共同决定。

[0100] 所述包含源点的边界波动积分方程能够由如图2所示的模型描述：

[0101] 均匀介质 由不规则边界Γ和内部区域Ω两部分组成，

[0102] r、r'和r0为 上的任意点，

[0103] 震源S(ω)位于r0处，

[0104] ω＝2πf为圆频率，

[0105] f为地震波频率，

[0106] 为二维背景介质中的格林函数，

[0107] 为汉克尔函数，

[0108] k0＝ω/v0为背景波数，

[0109] v0为背景波速，

[0110] C(r)＝θ(r)/2π为边界点r处与边界几何形状有关的常数，由r点处的边界切线面向 的夹角θ(r)决定，对于光滑边界，C(r)＝0.5。

[0111] 由于包含源点的边界波动积分方程(1)在整个计算区域 上成立，因此能够将r点和r'点移至边界上，将等式左边第一项拆分后移至等式右边，从而使边界上的函数u(r')及其法向导数 分离，便于单独处理，则上式转化为下式(2)

[0112]

[0113] 为叙述简便，令上式等号左边等于M，则得到下式(3)

[0114]

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

[0117] 其中，

[0118] u(rj')是边界单元Γj上任意点rj'处的函数值；

[0119] i,j＝1,2,…,s为边界单元序号；

[0120] s为边界单元总数。

[0121] 步骤2，任意选择一个步骤1确定的节点作为当前节点，任意选择一个步骤1确定的边界单元作为当前边界单元，用格林函数在当前边界单元上的外法向导数表征当前节点与当前边界单元上任意点之间的相互影响。

[0122] 在本发明步骤2中，基于步骤1得到的每个节点的空间位置坐标和每个边界单元端点的空间位置坐标，任意选择一个步骤1确定的节点作为当前节点，任意选择一个步骤1确定的边界单元作为当前边界单元，从选定的节点开始计算其与当前边界单元上任意点的格林函数，并求得到的格林函数在当前边界单元上的外法向导数，以此来表征当前节点与当前边界单元上任意点之间的相互影响。

[0123] 所述的格林函数如下式(5)所示，

[0124]

[0125] 其中，

[0126] 为0阶第一类汉克尔函数，

[0127] k0＝ω/v0为背景波数，

[0128] ω为圆频率，v0为背景波速，

[0129] r为当前节点位置，

[0130] r'为边界单元上任意点位置，

[0131] i和j表示边界单元序号，即i,j＝1,2,…,s。

[0132] 步骤3，用积分表征当前整个边界单元对当前节点处波场的影响。

[0133] 在本发明中，对步骤2得到的外法向导数进行积分，用得到的积分表征当前整个边界单元对当前节点处波场的影响，其中，所述的积分的表达式如下式(6)所示，[0134]

[0135] 其中，

[0136] Γj为序号为j的边界单元，

[0137] n是Γj的外法向。

[0138] 在本发明中，由于是用常单元对地形线和地质分界线进行的离散，因此有u(rj')＝u(rj)，则可将上式(6)进一步化为

[0139]

[0140] 由式(6)可知，上式(7)等号右边的积分项即为积分 则有

[0141]

[0142] 上式(8)表示对于当前计算节点i，需将其与所有边界单元(包括节点i所在单元本身)进行一次积分，对于有s个节点的边界，共需计算s次积分，对于当前节点i计算结束后依次移向下一个节点，直至所有节点计算完成，总积分次数为s2。

[0143] 本发明人注意到在上式(8)中，对所有i,j＝1,2,…,s，在j循环过程中会出现一次i＝j的情况，此时有u(ri)＝u(rj)，可以将其合并为 因此，重新定义

[0144]

[0145] 其中，

[0146] 为Kronecker函数，则式(8)中Mi可进一步化简为

[0147]

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

[0150] 从而表征了整个模型中不同产状的边界单元对不同节点处波场的影响。

[0151] 本发明提供的多层近地表地震地质复杂性的定量分析方法，主要利用上述矩阵方程的系数矩阵，即第一矩阵H进行散射复杂度分析。其中，

[0152] 所述的第一矩阵H为，

[0153]

[0154] 其中，

[0155] 矩阵的每个元素Hij由下式(13)计算得到，

[0156]

[0157] 其中，

[0158] lij为节点ri到任意点rj'的距离，

[0159] 为 的一阶导数；

[0160]

[0161] x1，y1和x1，y2是边界单元Γj两个端点的实际坐标；

[0162] K表示高斯积分点数；

[0163] ξk表示积分点位置；

[0164] wk表示积分点权系数；

[0165] lik表示边界单元Γi上节点ri与边界单元Γj上高斯积分点ξk之间的距离。

[0166] 步骤4，重复步骤2和步骤3，直至用积分表征了所有节点与所有边界单元之间的相互影响，形成第一矩阵H。

[0167] 在本发明中，重复步骤2和步骤3，从而计算了所有节点与所有边界单元之间的格林函数，并相应地求得格林函数在其边界单元上的外法向导数，进而用相应地高斯数值积分表征了每一个节点产生的波场之间的相互影响，即，表征了波场的复杂性，此时形成第一矩阵H。

[0168] 在本发明步骤4中，每一个节点和所有边界单元之间的相互影响都会一一被计算到，得到的信息全面而细致。

[0169] 步骤5，以步骤4中得到的第一矩阵H为基础，形成相应地散射振幅系数矩阵 并将散射振幅系数矩阵 的元素总和与矩阵维数的比值表示为多层近地表地震地质复杂性系数，实现多层近地表地震地质复杂性定量分析。

[0170] 由于H矩阵是一个复数矩阵，其维度表示边界单元的总数，矩阵中每个元素可以认为是从源节点传播到场节点的波场在场节点边界外法向上的分量，其与节点之间的距离、相对位置和节点所在边界的几何形态密切相关，而在复数H矩阵基础上对每个元素求模，将矩阵H转化为振幅矩阵 能够便于在实数域说明问题，从而实现对多层近地表地震地质复杂性定量分析。

[0171] 步骤5得到的多层近地表地震地质复杂性系数的数值与多层近地表地震地质复杂性呈正相关关系，即，其数值越大说明多层近地表地震地质复杂性越强。

[0172] 在本发明中，依据步骤5中得到的多层近地表地震地质复杂性系数选择地震成像方法、指导地震波数据的处理方法和/或选择地表起伏复杂地区的数据采集方法等，如利用频率域偏移方，当主频在30Hz时，在一般情况下，当多层近地表地震地质复杂性系数小于0.04时，采用列步傅里叶地震成像方法；当多层近地表地震地质复杂性系数介于0.04～
0.08时，采用一阶分离变量地震成像方法；当多层近地表地震地质复杂性系数大于0.08时，采用傅里叶有限差分地震成像方法。

[0173] 因此，在本发明步骤5中，在步骤4中得到的第一矩阵H的基础上计算散射振幅系数矩阵 从而将复数域波场转化为实数域振幅系数，以便得到标量复杂度系数，其中， 矩阵的计算方式如下式(14)所示：

[0174]

[0175] 在本发明中，由于整个模型的复杂性由所有边界单元之间的相互关系共同决定，因此根据上述矩阵元素的分布特征，本发明人定义了一个标量系数comp，来直接分析多层近地表地震地质复杂性，所述标量系数comp的定义公式如下式(15)所示，

[0176]

[0177] 其中，

[0178] comp表示多层近地表地震地质复杂性系数；

[0179] i和j表示边界单元序号，即i,j＝1,2,…,s；

[0180] |Hij|表示散射振幅系数矩阵 的元素；

[0181] s表示散射振幅系数矩阵 的维数。

[0182] 在本发明中，用散射振幅系数矩阵 的维数s来消除步骤1中边界单元划分数量对多层近地表地震地质复杂性系数产生的影响，进而，能够对不同规模模型的复杂性系数进行比较。

[0183] 根据本发明提供的多层近地表地震地质复杂性的定量分析方法，具有以下有益效果：

[0184] (1)整个模型的复杂性由所有边界单元之间的相互关系共同决定；

[0185] (2)对于同一个模型，可以针对不同频率的地震波给出不同的复杂度系数，实现了分频处理；

[0186] (3)该方法尤其适用于地表严重起伏和地下地质情况复杂区域，包括山地、高原、盆地等，可以较准确和细致地给出不同模型的构造复杂性评价参数；

[0187] (4)对地震数据的采集，地震数据处理中的去噪、偏移成像有着非常重要的实际意义；

[0188] (5)本发明提供的方法，以任意两点之间的基本解及其在边界上的外法向导数经过高斯数值积分计算得到的 矩阵，描述了边界上所有点对当前节点产生的散射影响，其优点在于包含了边界几何形态信息，为近地表构造复杂性分析提供了可能性。

[0189] 实验例

[0190] 实验例1不同高度的山峰地貌模型

[0191] 在本实验例中，设定山峰高度分别为160m、240m、320m和400m的四个模型，分别记为模型(a)、模型(b)、模型(c)和模型(d)，如图3～图6所示，四个模型的水平跨度均为2000m，边界单元数均为200，利用本发明提供的方法定量分析其复杂性，其相应 矩阵如图
7～图10所示，其灰度值随着山峰高度的增加而递增，其复杂性系数如图64中曲线1所示，也是随着山峰高度的增加而递增，递增的复杂性系数反映了复杂性随相对高差增加而增加。

[0192] 实验例2不同坡度的山峰地貌模型

[0193] 在本实验例中，设定山峰坡度逐渐增大的四个模型，其坡度分别为25°、35°、45°和65°，分别记为模型(a)、模型(b)、模型(c)和模型(d)，如图11～图14所示，利用本发明提供的方法定量分析其复杂性，其相应 矩阵如图15～图18所示，其灰度值随着山峰坡度的增加而递增，其复杂性系数如图64中曲线2所示，也是随着山峰坡度的增加而递增。

[0194] 实验例3不同水平跨度的地貌模型

[0195] 在本实验例中，设定四个模型不同水平跨度的地貌模型，分别为2000m、1500m、1000m和500m，分别记为模型(a)、模型(b)、模型(c)和模型(d)，如图19～图22所示，利用本发明提供的方法定量分析其复杂性，其相应 矩阵如图23～图26所示，其灰度值随着水平跨度的增加而递减，其复杂性系数如图64中曲线3所示，由图64中曲线3可知，水平跨度越小复杂性越强，即当山峰不变而水平部分无限延伸时，山峰的影响将逐渐减弱。

[0196] 实验例4不同山峰个数模型

[0197] 在本实验例中，设定四个模型不同山峰个数的地貌模型，分别为1个、2个、3个和4个，分别记为模型(a)、模型(b)、模型(c)和模型(d)，四个模型中每个山峰的高度、跨度和边界单元数分别为360m、2000m和200，如图27～图30所示，利用本发明提供的方法定量分析其复杂性，其相应 矩阵如图31～图34所示，其灰度值随着山峰个数的增加而剧烈增加，反映了多峰起伏地表具有很强的复杂性，这与实际地面地震观测中多峰起伏地表很强的散射噪音相符，其复杂性系数如图64中曲线4所示。

[0198] 实验例5不同地层数模型

[0199] 在本实验例中，设定四个不同地层数的地貌模型，地表及地层分界面均随机起伏，层内为均匀介质，如图35～图38所示，地层数依次为1层、2层、3层和4层，分别记为模型(a)、模型(b)、模型(c)和模型(d)，其中黑点表示边界单元的节点，其中，

[0200] 图35地下为均匀介质模型，图36、图37和图38的地层数分别为2层、3层和4层，各模型地形线、地质分界面形态及边界剖分方式完全相同。

[0201] 由图35～图38可知，随着非规则地层数的增加，模型越来越复杂。

[0202] 利用本发明提供的方法定量分析其复杂性，其相应 矩阵如图39～图42所示。

[0203] 在本实验例中，模型的边界形态呈无规律展布状态，因此对应的 矩阵分布表现得较为混乱，无法根据其特征评估模型构造特性。

[0204] 利用本发明提供的方法计算图35～图38所示模型的多层近地表地震地质复杂性系数，并将四个结果绘制成图64中曲线5，由图64中曲线5可知，前三种模型的复杂性系数分别为0.015144、0.031128、0.046847，三者基本呈线性增加趋势，这是因为地形线和两条非圈闭地质分界面的形态具有相似性，且地层平均厚度也不存在明显差异所致；而第四种模型的复杂性系数为0.073751，从图64中可以看出该段斜率明显增大，引起该变化的主要因素是模型中的圈闭层，尤其是底部和边界两侧的陡倾角位置是成像的难点，同时也是加剧复杂度的主要因素，这与理论预测相符。

[0205] 实验例6不同频率的地震波模型

[0206] 在本实验例中，设定的四个模型地形完全相同，均取自于实验例5中模型的第一层和第三层地质分界面，但地震波频率参数不同，分别为10Hz、35Hz、50Hz和65Hz，分别记为模型(a)、模型(b)、模型(c)和模型(d)，如图43～图46所示。

[0207] 相同的地质模型对于不同频率的地震波表现出不同的复杂度，理论推知，频率越高，波长越小，地震波越容易受构造影响，模型复杂度相对也就越大。

[0208] 利用本发明提供的方法定量分析其复杂性，其相应 矩阵如图47～图50所示，矩阵非零元素的分布完全相同，这是由于采用同一个模型的结果；而元素值的差异则体现了不同频率产生的影响，从图47～图50可以看出，矩阵的灰度值是逐渐变大的，可以初步判断整体复杂呈增加趋势。

[0209] 图65准确给出了利用本发明提供的方法得到的模型的多层近地表地震地质复杂性系数的结果，由图65可知，随着频率增加，复杂性明显上升，而且曲线的斜率与频率间隔近似为正相关。

[0210] 实验例7不同边界离散模型

[0211] 图51～图54所示的四个模型与图50所示模型完全相同，区别在于边界剖分精细程度(离散后所得边界单元个数)逐渐递增，剖分所得的边界单元数差别较大，分别为80、150、230和535，分别记为模型(a)、模型(b)、模型(c)和模型(d)，用于测试边界剖分对复杂系数计算的影响。

[0212] 模型离散要求边界单元与实际边界之间的误差尽量小，在起伏变化剧烈的地方用小而多的边界单元，而起伏平缓的地方则用长而少的边界单元，确保离散模型基本反映实际模型的主要特征。

[0213] 利用本发明提供的方法定量分析其复杂性，其相应 矩阵如图55～图58所示。

[0214] 利用本发明提供的方法得到的复杂性系数如图64中曲线6所示，复杂性系数基本维持在0.073左右，边界剖分对结果的影响十分微小。

[0215] 实验例8实际地下速度结构模型

[0216] 图59示出实际地下速度结构模型，其复杂性主要表现在两方面：严重起伏的地形和复杂的近地表结构。

[0217] 用本发明提供的方法分析和评价该地区的地质复杂性。

[0218] 首先根据速度分布勾勒出主要地层分界线，再利用常单元对地形线和地层分界线进行离散，节点和边界单元的标号遵循从左往右和从上到下的原则。

[0219] 离散结果如图59所示，共形成14条边界和445个边界单元，其中黑线表示边界，黑点表示边界单元端点。此外，为了对比分析，增加结果的可信度，在实际模型的基础上对地形线进行移动平滑，从而减小了地形对整个模型复杂度的影响，并在超出实际高程区域用邻域波速充填，以保证原速度结构的完整性，最后形成了图60所示的对比模型。

[0220] 显然，地形线平滑后的地质结构更加简单，用同样的方法划分和离散边界，边界数量仍为14，而边界单元随着地形线的去尖锐化减少为265。

[0221] 将本发明提供的方法应用于这两个模型，得到的 矩阵分别如图61和图62所示，由于本实验例中地形及分界面的随机性较强，矩阵缺少明确对比标志，因此仅凭感官难以准确区分两者的复杂度差异。选择频率10Hz、35Hz、50Hz和65Hz分别计算两者复杂性系数，结果如图63所示，

[0222] 图63中曲线1和曲线2分别对应图59和图60得到复杂性系数，从图63中可以看出，曲线1完全在曲线2的上方，由此说明在频率参数相同的情况下，用本发明提供的方法得出的复杂度系数，实际地下速度结构模型确实比表面平滑后的模型要复杂，而且，随着地震波频率的增大，两个模型的复杂度同时增加，这一点同样与理论相吻合。

[0223] 在本实验例中也能够看出，地表曲线的平滑对整个模型的整体影响并非很大，但本发明提供的方法从数值上准确反映了两者的差异，表现出对结构变化良好的响应。

[0224] 实验例9绘制复杂性系数曲线

[0225] 将实验例1～实验例8得到的复杂性系数曲线绘制于图64上，其中，[0226] 曲线1表示实验例1得到的复杂性系数；

[0227] 曲线2表示实验例2得到的复杂性系数；

[0228] 曲线3表示实验例3得到的复杂性系数；

[0229] 曲线4表示实验例4得到的复杂性系数；

[0230] 曲线5表示实验例5得到的复杂性系数；

[0231] 曲线6表示实验例7得到的复杂性系数。

[0232] 由图64可知，实验例1～实验例3所示模型地貌，其复杂性增加对应的复杂性系数递增斜率几乎相同，表明这些复杂性因素的影响属于同一个量级。

[0233] 此外，图64中曲线1到曲线3的初值逐渐升高，说明扁平山峰地貌复杂性弱于窄而高的山峰地貌。

[0234] 以上结合具体实施方式和范例性实例对本发明进行了详细说明，不过这些说明并不能理解为对本发明的限制。本领域技术人员理解，在不偏离本发明精神和范围的情况下，可以对本发明技术方案及其实施方式进行多种等价替换、修饰或改进，这些均落入本发明的范围内。本发明的保护范围以所附权利要求为准。