阅读说明：本技术 一种地壳三维结构模型融合方法及装置 (Method and device for fusing three-dimensional structure models of earth crust ) 是由 李卫东 赵晨曦 单新建 王亚兵 张定文 董前林 刘甲 吴峥嵘 时春波 于 2020-06-05 设计创作，主要内容包括：本发明提供一种地壳三维结构模型融合方法及装置,属于地球物理领域。该方法包括：在参与融合的多个地壳三维结构模型的重合区域上统一各模型的结构；在深度方向上将重合区域分成多个地层,使每个地层上结构统一后的各模型的物性参数均呈线性关系；结合协同克里金方法,以及每个地层上各模型内部的自协方差函数、各模型间的互协方差函数和各模型的权重比,得到相应地层的融合矩阵；求解每个地层的融合矩阵得到相应地层上结构统一后的各模型对应的最优加权系数矩阵；结合最优加权系数矩阵计算每个地层上所有点融合后的物性参数,得到重合区域上所有点融合后的物性参数,完成模型融合。本发明能消除各模型之间的结构差异和保留各模型原有的地质特征。(The invention provides a method and a device for fusing a three-dimensional structure model of a ground shell, and belongs to the field of geophysical. The method comprises the following steps: unifying the structure of each model on the superposition area of the plurality of crustal three-dimensional structure models participating in fusion; dividing the overlapping area into a plurality of stratums in the depth direction, so that the physical property parameters of the models with unified structures on each stratum are in a linear relation; combining a collaborative kriging method, and an autocovariance function inside each model on each stratum, a cross covariance function among the models and a weight ratio of each model to obtain a fusion matrix of the corresponding stratum; solving the fusion matrix of each stratum to obtain an optimal weighting coefficient matrix corresponding to each model with unified structure on the corresponding stratum; and calculating the physical parameters after all the points on each stratum are fused by combining the optimal weighting coefficient matrix to obtain the physical parameters after all the points on the overlapped area are fused, thereby completing model fusion. The invention can eliminate the structural difference among the models and keep the original geological characteristics of the models.)

一种地壳三维结构模型融合方法及装置

技术领域

本发明涉及一种地壳三维结构模型融合方法及装置，属于地球物理技术领域。

背景技术

随着全球不同区域地壳三维结构模型研究的增加，研究统一地壳结构模型的构建，以实现对区域地壳结构的拼接和整合，进一步提高地质解译准确性，为解译深部地壳结构提供一种具有区域完整性的三维物性结构模型具有重要意义。但是由于每种地壳三维结构模型是基于有限约束条件得到的最优解，模型本身存在多解性问题，而且，利用不同方法获得的不同地壳三维结构模型之间存在结构差异问题(模型之间分辨率、深度及属性的差异)，这些因素均限制了地壳结构的研究以及统一地壳结构模型的构建。

目前，多源异构重力、地震模型数据融合研究方法的提出主要依托于不同物性经验关系研究和重力、地震数据联合反演方法的发展。

在密度-速度的物性关系研究中，Birch、Ludwing、Gardner、冯锐和Christensen等通过理论研究和实验测量提出多种适用于不同条件的密度速度公式(详见以下文献：Birch在1961年Journal of Geophysical Research第66卷第7期，公开的一篇名为The velocityof compressional waves in rocks to 10kilobars:2的文献；Ludwig在1970年Deep SeaResearch&Oceanographic Abstracts第17卷第3期，公开的一篇名为The Manila Trenchand West Luzon Trough—III.Seismic-refraction measurements的文献；Gardner在1974年Geophysics第39卷第6期，公开的一篇名为Formation velocity and density-thediagnostic basics for stratigraphic traps的文献；冯锐在1985年地震学报第2期，公开的一篇名为中国地壳厚度及上地幔密度分布(三维重力反演结果)的文献；Christensen在1996年Journal of Geophysical Research Solid Earth第101卷第B2期，公开的一篇名为Poisson's ratio and crustal seismology的文献)。Brocher等对比分析多种密度-弹性波速度关系的适用条件，并给出了适用于多种条件的地壳纵横波转换公式(详见以下文献：Brocher在2005年Bulletin of the Seismological Society of America第95卷第6期，公开的一篇名为Empirical Relations between Elastic Wavespeeds and Densityin the Earth’s Crust的文献)，这为模型融合提供了统一的属性基础。

在联合反演方面，以地质统计学方法为约束的地球物理反演将先验知识融入协同克里金矩阵中，能够有效提高反演结果准确性(例如以下文献中构建的反演矩阵：Asli在2000年Proceedings of the 6th International Geostatistics Congress-Geostat2000.Marshalltown:Geostatistical Association of South Africa，公开的一篇名为Direct Inversion of Gravity Data by Cokriging的文献；Shamsipour在2010年Geophysics第75卷第1期，公开的一篇名3D Stochastic Inversion of Gravity DataUsing Cokriging and Cosimulation为的文献；Shamsipour在2011年Geophysics第73卷第4期，公开的一篇名为3D Stochastic Inversion of MagneticData的文献；Geng作者在2014年Geophysics第79卷第4期，公开的一篇名为3D Inversion of Airborne Gravity-Gradiometry Data Using Cokriging的文献；耿美霞在2016年地球物理学报第59卷第5期，公开的一篇名为基于协同克里金方法的重力梯度全张量三维约束反演的文献；高秀鹤在2017年吉林大学学报(地球科学版)第47卷第2期，公开的一篇名为重力梯度数据协克里金三维反演确定岩脉倾向的文献；高秀鹤在2019年地球物理学报第62卷第3期，公开的一篇名为基于阈值约束的协克里金法联合反演重力与重力梯度数据的文献)。但这些反演矩阵的构建是以期望存在且已知为假设，在解决多源异构地壳三维结构模型重合区域的统一表达问题时，这种反演矩阵无法附加不同模型相关性的约束条件，也无法融合不同地壳结构模型的地质特征。

综上，在融合多源异构地壳三维结构模型时，如何在消除模型间结构差异的同时保留各模型原有的地质特征仍是需要解决的问题。

发明内容

本发明的目的是提供一种地壳三维结构模型融合方法及装置，用以解决目前融合地壳三维结构模型时，无法同时消除模型间结构差异和保留各模型原有的地质特征的问题。

为实现上述目的，本发明提供了一种地壳三维结构模型融合方法，该方法包括以下步骤：

(1)在参与融合的多个地壳三维结构模型的重合区域上对各模型的结构进行统一，包括对各模型的物性参数和各模型的分辨率进行统一；

(2)在深度方向上将所述重合区域分成多个地层，使每个地层上结构统一后的各模型的物性参数均呈线性关系；

(3)基于协同克里金方法分层构建融合矩阵，其中，通过对某个地层上结构统一后的各模型进行处理得到该地层的融合矩阵的过程为：

将该地层上各模型的三维变异函数作为该地层上各模型内部的自协方差函数；

根据该地层上各模型的物性参数的线性关系和各模型内部的自协方差函数构建该地层上各模型间的互协方差函数；

根据该地层上各模型的物性参数的线性关系确定该地层上各模型的权重比；

结合协同克里金方法，以及该地层上各模型内部的自协方差函数、各模型间的互协方差函数和各模型的权重比，得到该地层的融合矩阵；

(4)对每个地层的融合矩阵进行求解，得到相应地层上结构统一后的各模型对应的最优加权系数矩阵；针对每个地层上的任一点，根据该点以及该点邻域内的各点对应的不同模型的物性参数和该点邻域内的各点对应的不同模型的最优加权系数，得到该点融合后的物性参数，从而得到相应地层上所有点融合后的物性参数，进而得到重合区域上所有点融合后的物性参数，完成模型融合。

本发明还提供了一种地壳三维结构模型融合装置，该装置包括处理器和存储器，所述处理器执行由所述存储器存储的计算机程序，以实现上述的地壳三维结构模型融合方法。

本发明的地壳三维结构模型融合方法及装置的有益效果是：首先对重合区域上各模型的结构进行统一，以消除模型间的结构差异；其次，在深度方向上将重合区域分成多个地层，分层原则是保证在每个地层上结构统一后的各模型的物性参数均呈线性关系，以满足二阶平稳假设提高模型融合的可靠性；接着，基于协同克里金方法分层构建融合矩阵，由于各模型内部的自协方差函数体现了各模型内部数据间的自相关性，则利用自协方差函数和获得的线性关系对互协方差函数加以约束，能在保留各模型原有先验知识的同时提高互协方差函数的精度，从而将各模型内部的自协方差函数和各模型间的互协方差函数加入融合矩阵能够保留各模型原有的地质特征；各模型的权重比体现了各模型间的相关性，从而将各模型的权重比加入融合矩阵能够体现不同模型地质特征之间的相关性；最后，通过求解重合区域上所有点融合后的物性参数完成模型融合。综上，本发明不仅能消除各模型之间的结构差异，还能保留各模型原有的地质特征，并能体现不同模型地质特征之间的相关性。

进一步地，在上述地壳三维结构模型融合方法及装置中，某个地层的融合矩阵为：

式中，λ₁，λ₂分别为模型U和模型V对应的最优加权系数矩阵，γ表示模型U和模型V间的线性关系，为模型U的权重比，

为模型V的权重比，cov(u_k,u_j)为模型U的自协方差函数，cov(v_k,v_j)为模型V的自协方差函数，cov(u_k,v_j)为模型U和模型V间的互协方差函数，cov(v_k,u_j)为模型V和模型U间的互协方差函数，

u_k表示点k对应的模型U的物性参数，u_j表示点k邻域内的第j个点对应的模型U的物性参数，n表示点k邻域内点的总个数，v_k表示点k对应的模型V的物性参数，v_j表示点k邻域内的第j个点对应的模型V的物性参数，u_o表示点o对应的模型U的物性参数，v_o表示点o对应的模型V的物性参数，μ₁、μ₂是构建融合矩阵的拉格朗日参数。

进一步地，在上述地壳三维结构模型融合方法及装置中，在某个地层上以幂函数模型分别拟合结构统一后的各模型的x轴、y轴和z轴三个轴向的变异函数，然后采用各向异性套合方法对各模型三个轴向的变异函数进行套合，得到该地层上各模型的三维变异函数。

进一步地，在上述地壳三维结构模型融合方法及装置中，所述三维变异函数为：γ^*(h)＝γ₁(h)+γ₀(|Ah|)-γ₁(|Ah|)，式中，h为样本点间距，γ₁(h)为水平方向变异函数，γ₀(|Ah|)-γ₁(|Ah|)为垂直方向变异函数，A为距离变换矩阵，γ^*(h)为所述三维变异函数。

进一步地，在上述地壳三维结构模型融合方法及装置中，所述结构统一后的各模型的物性参数为P波速度或S波速度或密度。

进一步地，在上述地壳三维结构模型融合方法及装置中，当所述结构统一后的各模型的物性参数为P波速度时，利用横纵波经验关系式和速度-密度经验关系式将各模型的物性参数统一为P波速度，所述横纵波经验关系式为：

所述速度-密度经验关系式为：

式中，v_p为P波速度，v_s为S波速度，ρ为密度。

附图说明

图1是本发明方法实施例中的地壳三维结构模型融合方法流程图；

图2是本发明方法实施例中的华北克拉通中部山西断陷带区域图；

图3-a是本发明方法实施例中在深度0km-1.4km开展线性关系拟合得到的东、西部模型线性关系图；

图3-b是本发明方法实施例中在深度1.4km-15km开展线性关系拟合得到的东、西部模型线性关系图；

图3-c是本发明方法实施例中在深度15km-39km开展线性关系拟合得到的东、西部模型线性关系图；

图3-d是本发明方法实施例中在深度39km-50km开展线性关系拟合得到的东、西部模型线性关系图；

图4-a是本发明方法实施例中东部模型在重合区域的地壳三维P波速度结构模型图；

图4-b是本发明方法实施例中西部模型在重合区域的地壳三维P波速度结构模型图；

图4-c是本发明方法实施例中Krige模型在重合区域的地壳三维P波速度结构模型图；

图4-d是本发明方法实施例中G模型在重合区域的地壳三维P波速度结构模型图；

图5-a是本发明方法实施例中NCISP4剖面位置的西部模型图；

图5-b是本发明方法实施例中NCISP4剖面位置的东部模型图；

图5-c是本发明方法实施例中NCISP4剖面位置的Krige模型图；

图5-d是本发明方法实施例中NCISP4剖面位置的G模型图；

图6-a是本发明方法实施例中文登-阿拉善剖面位置的东部模型图；

图6-b是本发明方法实施例中文登-阿拉善剖面位置的西部模型图；

图6-c是本发明方法实施例中文登-阿拉善剖面位置的Krige模型图；

图6-d是本发明方法实施例中文登-阿拉善剖面位置的G模型图；

图7-a是本发明方法实施例中盐城-包头剖面位置的东部模型图；

图7-b是本发明方法实施例中盐城-包头剖面位置的西部模型图；

图7-c是本发明方法实施例中盐城-包头剖面位置的Krige模型图；

图7-d是本发明方法实施例中盐城-包头剖面位置的G模型图；

图8-a是本发明方法实施例中东部模型截取剖面BB`位置的东部模型图；

图8-b是本发明方法实施例中东部模型截取剖面BB`位置的西部模型图；

图8-c是本发明方法实施例中东部模型截取剖面BB`位置的Krige模型图；

图8-d是本发明方法实施例中东部模型截取剖面BB`位置的G模型图；

图9是本发明装置实施例中的地壳三维结构模型融合装置结构图。

具体实施方式

方法实施例

本实施例的地壳三维结构模型融合方法如图1所示，下面以参与融合的地壳三维结构模型模型有2个(以下称模型U和模型V)为例，详细介绍利用该方法进行地壳三维结构模型融合的过程，具体如下：

步骤1、在模型U和模型V的重合区域上对这两个模型的结构进行统一，包括对模型U和模型V的物性参数和分辨率进行统一。

由于利用不同方法获得的不同地壳三维结构模型之间存在结构差异问题，利用该步骤能消除模型间的结构差异，为之后的模型融合奠定基础。

当前基于正、反演方法获取并构建的地壳三维结构模型主要包含P波速度、S波速度和密度的结构模型，一些模型在反演结果中已给出经过转换的P波速度结构模型，因此，本实施例中，选用P波速度作为统一后的物性参数，并利用横纵波经验关系式和速度-密度经验关系式将模型U和模型V的物性参数统一为P波速度，利用插值方法将模型U和模型V的分辨率进行统一。

横纵波经验关系式为：

速度-密度经验关系式为：

式中，v_p为P波速度，v_s为S波速度，ρ为密度。

作为其他实施方式，还可以选用S波速度或密度作为统一后的物性参数，并选择相应的物性关系经验公式将模型的物性参数统一为S波速度或密度。

步骤2、在深度方向上将模型U和模型V的重合区域分成多个地层，使每个地层上模型U和模型V的P波速度均呈线性关系。

由于重合区域地球物理场存在同源性，为满足二阶平稳假设，对不同模型的重合区域进行深度分层处理，使每个地层上各模型均满足正态分布假设，能够提高模型融合的可靠性。

理论上地壳结构内相同地层具有空间一致性，因此对于不同地壳三维结构模型在相同位置存在的差异可利用线性关系表达：

u＝γ·v (3)

式中，u为模型U的P波速度，v为模型V的P波速度，γ表示模型U和模型V间的线性关系。

步骤3、基于协同克里金方法分层构建融合矩阵。

其中，通过对某个地层上结构统一后的模型U和模型V进行处理得到该地层的融合矩阵的过程为：

(1)分别构建该地层上模型U和模型V内部的自协方差函数；

本实施例中，考虑到地壳三维结构模型存在各向异性，首先基于该地层走向确定该地层的三个主轴方向(即x轴、y轴和z轴三个轴向)，然后以幂函数模型分别拟合模型U的x轴、y轴和z轴三个轴向的变异函数，最后采用各向异性套合方法对模型U三个轴向的变异函数进行套合，得到该地层上模型U的三维变异函数，并将该地层上模型U的三维变异函数作为模型U内部的自协方差函数；同理，该地层上模型V的三维变异函数可类似得到，并将该地层上模型V的三维变异函数作为模型V内部的自协方差函数。

其中，各轴向的变异函数以及三维变异函数的计算方法如下：

由于地壳结构模型所占区域尺度较大，而地壳融合以地震波P波速度作为属性，地层速度差变化一般小于0.2km/s，因此本实施例中采用幂函数模型(式A1)作为计算变异函数的统一模型，采用最小二乘法拟合变异函数模型，并以变异系数检验验证拟合函数精度。

γ(h)＝scale*h^exponet+nugget (A1)

式中，scale是变异函数比例系数，exponet是变异函数指数系数，nugget是变异函数块金值，h为两样本点间距。

x、y、z三轴向变异函数套合方法采用各向异性套合，该方法将水平方向上变异函数特征作为三维空间中各项同性特征γ₁(|h|)，整个三维空间中套合变异函数γ^*(h)是在垂直方向以γ₂(h_z)对各水平方向的套合。套合后得到的三维变异函数如式A2所示：

γ^*(h)＝γ₁(h)+γ₀(|Ah|)-γ₁(|Ah|) (A2)

式中，h为样本点间距，γ₁(h)为水平方向变异函数，γ₀(|Ah|)-γ₁(|Ah|)为垂直方向变异函数，A为距离变换矩阵，γ^*(h)为套合后得到的三维变异函数。

由于变异函数可以定量表达和分析样本数据的空间自相关性，因此将三维变异函数作为模型内部的自协方差函数能够体现模型内部数据间的自相关性，从而能够体现模型原有的地质特征。

(2)构建该地层上模型U和模型V间的互协方差函数；

其中，根据该地层上模型U和模型V的P波速度的线性关系，以及模型U和模型V内部的自协方差函数构建模型U和模型V间的互协方差函数，具体如下：

若模型U和模型V的P波速度存在如公式(3)所示的线性关系，则自协方差函数和互协方差函数存在如下关系：

cov(u_k,v_j)＝γ·cov(u_k,u_j) (4)

式中，cov(u_k,u_j)为模型U的自协方差函数，cov(u_k,v_j)为模型U和模型V间的互协方差函数。

本实施例中，在计算该地层上点o处模型U和模型V间的互协方差函数时，对点o处的互协方差函数cov(u_o,v_k)和cov(v_o,u_k)进行变形，加入模型U、模型V在点o处的自协方差函数cov(u₀,u_k)和cov(v₀,v_k)，以此利用原模型具有的先验知识对模型间的互协方差函数加以约束，以保留各模型原有的地质特征，同时提高互协方差函数精度。变形如下：

(3)确定该地层上模型U和模型V的权重比；

本实施例中，根据该地层上模型U和模型V的P波速度的线性关系确定模型U和模型V的权重比；若该地层上模型U和模型V的P波速度存在如公式(3)所示的线性关系，则模型U的权重比为

模型V的权重比为

(4)结合协同克里金方法，以及该地层上模型U和模型V内部的自协方差函数、模型U和模型V的互协方差函数和模型U和模型V的权重比，得到该地层的融合矩阵。

本实施例中，考虑到不同正、反演方法获取的地壳结构模型在地层深度和地质特征三维空间展布上存在差异，为了使模型融合后能够保留模型U和模型V原有的地质特征，将模型U和模型V内部的自协方差函数和模型U和模型V间的互协方差函数加入融合矩阵；并且，为了使模型融合后能够体现模型U和模型V间的相关性，将模型U和模型V的权重比作为约束条件加入融合矩阵，在此基础上结合协同克里金方法最终得到的某个地层的融合矩阵为：

式中，λ₁，λ₂分别为模型U和模型V对应的最优加权系数矩阵，γ表示模型U和模型V间的线性关系，为模型U的权重比，为模型V的权重比，cov(u_k,u_j)为模型U的自协方差函数，cov(v_k,v_j)为模型V的自协方差函数，cov(u_k,v_j)为模型U和模型V间的互协方差函数，cov(v_k,u_j)为模型V和模型U间的互协方差函数，

u_k表示点k对应的模型U的P波速度，u_j表示点k邻域内的第j个点对应的模型U的P波速度，n表示点k邻域内点的总个数，v_k表示点k对应的模型V的P波速度，v_j表示点k邻域内的第j个点对应的模型V的P波速度，u_o表示点o对应的模型U的P波速度，v_o表示点o对应的模型V的P波速度，μ₁、μ₂是构建融合矩阵的拉格朗日参数。

步骤4、对每个地层的融合矩阵进行求解，得到相应地层上结构统一后的各模型对应的最优加权系数矩阵；针对每个地层上的任一点，根据该点以及该点邻域内的各点对应的不同模型的P波速度和该点邻域内的各点对应的不同模型的最优加权系数，得到该点融合后的P波速度，从而得到相应地层上所有点融合后的P波速度，进而得到重合区域上所有点融合后的P波速度，完成模型融合(即得到最终的融合模型)。

例如某个地层上点o融合后的P波速度的计算过程如下：

对该地层的融合矩阵进行求解，得到该地层上结构统一后的各模型对应的最优加权系数矩阵，即模型U、模型V对应的最优加权系数矩阵λ₁、λ₂；

根据点o以及点o邻域内的各点对应的不同模型(即模型U和模型V)的P波速度，以及点o邻域内的各点对应的不同模型(即模型U和模型V)的最优加权系数，得到点o融合后的P波速度。

当参与融合的模型有3个时，先随机选择两个模型按照本实施例的方法进行融合，得到一个新模型，然后再按照本实施例的方法将得到的新模型和剩下的模型进行融合，得到最终的融合模型。当参与融合的模型>3个时，融合过程与此类似，不再赘述。

下面基于实际数据进行实验，并分别利用统计学评价方法和实测剖面评价方法对本实施例的地壳三维结构模型融合方法(以下简称本实施例方法)的有效性进行验证。

基于统计学评价方法的验证实验：

以如图2所示的华北克拉通中部山西断陷带区域为例，利用本实施例方法构建研究区统一地壳三维P波速度结构模型，并引入加权平均法构建简单的融合模型作为比对模型。验证方法采用均方根误差和pearson相关系数。

对华北克拉通中部山西断陷带区域已获取的两个地壳三维结构模型开展模型融合实验，参与模型融合的三维地壳速度结构模型有：利用Krige插值方法构建的华北克拉通中东部地壳三维P波速度结构模型；利用接收函数方法获取的华北克拉通西部地壳三维S波速度结构模型，S波模型已给出转换后的P波速度。在实际检验中，由于三维地壳结构模型尺度较大且包含大量地质构造特征，直接提取并融合重合区域数据具有难以符合二阶平稳假设和变异函数难以拟合的难题，所以需要依照重合区域主要地质构造对重合区域进行分块，本实施例以山西断陷带作为该区域主要地质构造特征提取对应模型的待融合数据。

华北克拉通中东部地壳三维P波速度结构模型(以下简称东部模型)在深度10km的分辨率为0.25°×0.25°×1km，在深度10km-50km的分辨率为0.25°×0.25°×2km。华北克拉通西部地壳三维S波速度结构模型(以下简称西部模型)在深度7km的分辨率为0.2°×0.2°×0.2km，在深度7km-50km的分辨率为0.2°×0.2°×1km。东部模型和西部模型给出的P波速度结构分层如表1所示，表中East Model为东部模型，West Model为西部模型，主要地层包括沉积层(Sediment)、上地壳(Upper Crust)、中地壳(Middle Crust)、下地壳(LowerCrust)、壳幔过渡带(Crust-mantle transition zone)和岩石圈地幔(Lithosphericmantle)。

表1东部模型和西部模型各主要地层P波速度

本实施例中选取经度112°-113°、纬度36°-42°范围作为待融合区(见图2中虚线框区域)，该区域主要地质构造特征以山西断陷带为主。以分辨率较高的西部模型网格作为待融合模型网格，利用IDW方法(即反距离加权插值方法)对东部模型进行插值，统一各模型分辨率，最终得到东部模型与西部模型分层速度结构。

为满足二阶平稳假设，分别以0km-1.4km深度、1.4km-15km深度、15km-39km深度和39km-50km深度对模型进行分层。选取各层主要P波速度特征进行正态检验获取最大样本点数量，并均匀抽取各层包含P波速度属性样本点坐标。由于区域内盆地地质特征呈北东向，现定义区域变量特征方向为[1,1,1]，区域坐标系三主轴方向分别为[1,0,0]，[0,1,0]，[0,0,1]，并以幂函数模型拟合x、y、z三轴向变异函数，利用套合方法获取最终三维变异函数。

在定义融合模型关系时，分层拟合各模型对应地层线性关系，计算互协方差函数，最终构建基于协同克里金方法的融合矩阵，对各点进行邻域点搜索和矩阵求解，获取研究区统一地壳三维P波速度结构模型(即利用本实施例方法得到的融合模型)。其中，依据深度0km-1.4km、1.4km-15km、15km-39km和39km-50km开展线性关系拟合得到的东、西部模型线性关系分别如图3-a、图3-b、图3-c和图3-d所示。

在分析验证时，加入G模型作对比分析。G模型以东部模型和西部模型的线性关系作为权重系数，对各点速度进行加权求和，东部模型和西部模型所占比重为

和

图4-a至图4-d为参与融合的模型与融合模型在重合区域的地壳三维P波速度结构模型对比图，其中，东部模型和西部模型为参与融合的模型，Krige模型表示利用本实施例方法得到的融合模型，G模型表示利用加权方法得到的融合模型。结合图4-a至图4-d可知，东、西部模型P波速度在5.5km/s以下、6km/s、6.5km/s及7.5km/s存在峰值，由此可以将东、西部模型划分为沉积层、上地壳、下地壳和上地幔四个符合二阶平稳假设的地层结构，而最终通过本实施例方法和加权方法得到的Krige模型和G模型也符合四层结构。

模型精度评估采用均方根误差(RMSE))和Pearson相关系数作为评价指标，表2和表3是对比分析结果和各模型间RMSE和Pearson系数值。

表2各地层RMSE

东部模型和西部模型为不同方法获取的重合区域地壳速度结构模型，由表2可以看出，模型均方根误差在0.342，误差较大。而融合后的Krige模型对于东部模型和西部模型的均方根误差均小于0.3，仅在16km-39km深度误差上升，证明Krige能够近似满足对东、西部模型的表达，而G模型的均方根误差小于0.2，精确程度高于Krige模型。

表3各地层pearson系数

表3中Pearson相关系数为分层后各模型相关性分析。未融合前，东部模型和西部模型的相关性仅在39km-50km处大于0.7。融合后，Krige模型与东部模型的相关性显著提高，而G模型作为东部模型和西部模型利用加权平均方法得到的融合模型，其相关性高于Krige模型。

综上，基于统计学评价方法得到的评价结果为：利用本实施例方法构建的融合模型(即Krige模型)在均方根误差和Pearson相关系数上的精度低于加权平均法构建的融合模型(即G模型)。

但是由于进行多源异构地壳模型融合时，希望融合后的模型能保留多个不同模型原有的地壳结构特征及地质体特征，因此选用的评价指标应该能反映出融合后的模型是否能体现参与融合的各模型原有的地质特征，而均方根误差和Pearson相关系数是统计学中通用的评价指标，这些指标可以普遍应用在任何领域，无法体现多源异构地壳模型融合的特殊性，依据这些指标并不能真实地反映出融合后的模型是否实现了对多源异构地壳模型地质特征的融合表达。因此，需要选用新的评价指标进行进一步验证。

基于实测剖面评价方法的验证实验：

现采用三条实测剖面(分别见图2中NN`、YY`、WW`)和一条东部地壳结构模型截取剖面(见图2中BB`)对Krige模型地质特征进行分析，以加权方法获取的G模型作为对比验证模型，确定本实施例方法的最终可行性。

图2中，BB`表示东部地壳结构模型截取剖面，NN`表示被动源接收函数方法获取的NCISP4剖面，YY`表示深地震方法获取盐城-包头剖面，WW`表示深地震方法获取的文登-阿拉善剖面。

(1)被动源接收函数方法获取的NCISP4剖面比对验证结果

NCISP4是利用接收函数方法获取的地壳S波速度结构成像。该成像区域地层起伏，地壳内有倾斜及水平展布的低速层，Moho面最深达到45km，保留华北克拉通演化留下的构造痕迹，揭示华北克拉通自西向东地壳减薄构造特征。该区域最主要特征为下地壳中低速层L2，它是东部古洋壳带入岩石圈顶部的古陆壳折返到下地壳的残留体。

图5-a至图5-d是NCISP4剖面位置的模型对比图，结合图5-a至图5-d可知，西部模型和东部模型在沉积层、康拉德面的深度及形态上存在差异。西部模型在剖面西侧30km深度附近存在P波速度小于6.2km/s的低速层L2，东部模型在康拉德面以下存在P波速度为6.45km/s的凹陷。

Krige模型的沉积层近似于东部模型，康拉德面则介于于东、西部模型康拉德面之间，在深度25km-40km存在低速异常L2`，P波速度低于6.2km/s，可与西部模型中低速带L2对应。

G模型沉积层近似于西部模型，康拉德面为P波速度为6.3km/s的等值线，近似于东部模型，不存在对低速层L2的描述。对于莫霍面形态及深度的表达，东、西模型相近，而Krige、G两个融合模型莫霍面则介于它们之间。

(2)深地震方法获取文登-阿拉善剖面比对验证结果

东部模型中的文登-阿拉善P波速度结构剖面是采用二维射线追踪方法对折射波和反射波建模得到的二维深部地震测深剖面。介于山西断陷带与太行山***交界处，由G、C1、C3和M四个界面对沉积层、上地壳、中地壳、下地壳和上地幔，而C2界面则对中地壳再次进行分层。沉积层厚度为0.5km–2km，地壳厚度为39km，下地壳的厚度为15km–19km。Moho深度自西向东增加至40km以上。地壳内构造异常复杂多变，中上地壳有明显高、低速异常结构。

图6-a至图6-d是文登-阿拉善剖面位置的模型对比图，结合图6-a至图6-d可知，东、西部模型仅在沉积层与莫霍面界面存在深度及形态相似特征，沉积层均位于深度1.4km附近且界面平稳，莫霍面均位于深度40km附近，自西向东逐渐上升。东部模型中、上地壳分界面C1位于深度15km，中、下地壳分界面C2位于深度23km，下地壳内部C3界面位于深度30km。上地壳中存在P波速度大于6.2km/s的高速异常，中地壳中存在P波速度低于6.15km/s的低速异常。西部模型中、上地壳分界面C1位于深度10km处，中、下地壳分界面C2位于深度20km处，C3界面深度位于30km-40km之间，剖面东侧下地壳中存在P波速度低于6.49km/s的低速层，C3界面是由西侧P波速度为6.5km/s的低速层顶部与东侧P波速度为6.7km/s的等值线组成。

Krige模型沉积层界面高于东、西部模型沉积层界面，莫霍面近似东部模型莫霍面。C1、C2界面与东部模型结构近似，C2、C3界面深度略低于东部模型。C1、C2界面为6.15km/s的P波速度等值面，因此介于C1和C2界面之间，深度5km-15km的地层为高速层，而深度在15km-30km的地层为低速层，与东部模型文登-阿拉善剖面的中、上地壳高、低速异常结构交错对应。且在C3界面以下，存在与西部模型对应的低速异常。

G模型的沉积层和莫霍面与东、西部模型近似。C1界面低于东、西部模型C1界面深度，C2界面与Krige模型C2界面近似，C3界面在形态与深度上与西部模型P波速度为6.7km/s的C3界面相同。但G模型C1、C2界面以上P波速度自上而下呈递增变化，仅在P波速度为6.2km/s等值线以下存在低速异常，未保留东部模型中、上地壳的高、低速异常特征。剖面西侧P波速度低于6.5km/s的低速层位于C3界面以上，与西部模型低速层位置不对应。

(3)深地震方法获取盐城-包头剖面比对验证结果

东部模型中的盐城-包头P波速度结构剖面是采用射线追踪方法反演高分辨率折射和宽角反射/折射联合探测剖面资料得到的岩石圈二维P波速度结构剖面。该剖面分层结构同文登阿拉善一致，剖面中山西断陷带遭受强烈地质构造变形,地壳结构破碎，上地壳存在明显低速凹陷区。在中地壳上部C1界面与C2界面之间存在低速体。

图7-a至图7-d是盐城-包头剖面位置的模型对比图，结合图7-a至图7-d可知，东、西部模型莫霍面深度位于45km且形态差异不大，剖面中部沉积层深度轻微增加。东部模型沉积层深度位于1.4km左右，C1界面深度位于15km-20km之间，C2界面位于深度20km处，C3界面位于深度30km处。上地壳存在明显低速凹陷区，中地壳为P波速度低于6.15km/s的低速带，C3界面以下存在P波速度低于6.6km/s的低速层。西部模型沉积层深度处于0km-5km之间，C1界面位于5km-15km之间，C2界面位于20km-30km之间，上下界面速度差大于0.2km/s。C3界面深度位于35km处，C3界面以下地壳结构为P波速度小于6.55km/s的低速带。

Krige模型沉积层仅中部存在凹陷，深度小于5km，由剖面中部开始向两端减薄，存在上地壳出露，与东、西部模型差异较大，莫霍面深度位于45km处，形态与东部模型近似。C1界面深度位于10km-15km之间，C2界面自西向东逐渐下降，深度位于20km-30km之间，C3界面深度位于30km处，东侧接近C2界面。位于C1、C2界面之间的低速带与东部模型对应，20km-30km的高速带使C2界面抬高，但与西部模型C2、C3界面之间的地壳结构相对应，且在C3界面以下40km深度存在低速层，能够与东部模型的低速凹陷与西部模型低速带相对应。

G模型主要间断面(沉积层、C1、C2、C3、莫霍面)的深度和形态与西部模型相近，在C3界面以下存在P波速度小于6.53km/s的低速层，但中上地壳结构缺少对东部模型中低速带的描述。

(4)华北克拉通中东部地壳三维速度结构模型截取BB`剖面比对验证结果

东部模型中BB`剖面为华北克拉通中东部地壳三维速度结构模型HBCrust1.0中的纵剖面。该剖面与NCISP4相交，由G、C和M三层间断面将模型分为沉积层、上地壳、下地壳和上地幔。剖面中部分地壳出露至地表，山西断陷带低速层从16km延续至30km深度。

图8-a至图8-d是东部模型截取剖面BB`位置的模型对比图，结合图8-a至图8-d可知，东、西部模型莫霍面与康拉德面形态近似，其康拉德面深度在20km-30km之间，莫霍面自西向东逐渐上升，深度在45km附近。东部模型BB`剖面沉积层位于深度5km附近，存在上地壳出露，下地壳自西向东逐渐减薄，存在于深度35km至39km之间。下地壳中P波速度为6.55km/s的等值线自西向东上升，对比康拉德面可以发现在下地壳西部区域存在凹陷。西部模型沉积层深度在5km左右，与东部模型差异较大。在深度20km-30km处，康拉德面以下的下地壳中存在低速层。

Krige模型沉积层在剖面东侧保持地壳出露特征，符合东部模型沉积层深度与形态特征，在上地壳中，存在P波速度为6.0km/s的低速异常，符合东部模型康拉德面在剖面西侧的凹陷特征。在深度25km至40km下地壳中存在低速异常近似表达同一位置东部模型P波速度为6.55km/s的凹陷形态和西部模型康拉德面以下P波速度为6.25km/s的低速层。

G模型康拉德面、莫霍面的深度与形态与东、西部模型相似，沉积层深度、形态符合西部模型特征。在剖面东侧上地壳中存在P波速度为6.15km/s的特征，仅符合对西部模型的描述。在P波速度为6.4km/s的等值线以下，深度在30km以下的G模型地壳中不存在任何速度异常。

综上，基于实测剖面评价方法的评价结果为：通过对比东部模型两条深地震测深剖面、一条截取剖面和西部模型一条接收函数成像剖面可以发现，在东、西部模型差异不大的莫霍面和康氏面处，Krige模型和G模型能够良好保留原有模型间断面特征。当东、西部模型存在差异时，Krige模型保留的主要间断面形态与东部模型近似，G模型保留的主要间断面形态与西部模型近似。在对东西部高速异常和低速异常的表达上，G模型仅能保留西部模型下地壳结构特征，而Krige模型能够通过校正主要间断面形态及P波速度差异同时保留东、西部模型的高速异常与低速异常，因此利用本实施例方法获得的融合模型具有更高的可信度。

综合统计学评价方法和实测剖面评价方法可以发现，虽然本实施例方法构建的融合模型在均方根误差和相关系数上的精度低于加权平均法构建的融合模型，但在比对实测剖面时，本实施例方法能够保留研究区多种地壳三维结构模型的地质特征，克服模型融合中地质特征差异难以融合的问题，因此本实施例方法是更具有可行性和可信度的。

本实施例的地壳三维结构模型融合方法，将地质统计学原理与地球物理方法相结合，解决了多源异构地壳三维结构模型在模型融合与地壳结构统一构建上的难题。针对不同地球物理方法获取的地壳三维结构模型在分辨率、深度及属性上的差异，该方法运用地震波速度与密度经验公式和相关插值方法得到具有相同分辨率、深度及属性的地壳三维P波速度结构模型，在构建保留多模型地质特征的统一地壳三维结构模型时，利用基于协同克里金方法的模型融合方法可以将地质特征融入融合矩阵中获取包含多模型地质特征的统一地壳结构。

该方法主要解决反演地壳结构模型的多解性问题，以及地震波、重力及二者联合反演获取的多种地壳三维结构模型差异整合问题。随着对地壳结构研究深入，利用不同方法获取的地壳三维结构模型越来越多，而各模型之间难以将特征融合开展数值模拟和地质解译。本实施例的地壳三维结构模型融合方法，旨在为多种地球物理方法交叉融通，以统一属性、融合建模、联合模拟的方法开展地壳结构研究提供一套模型融合解决方案。

装置实施例

本实施例的地壳三维结构模型融合装置，如图9所示，该装置包括处理器、存储器，存储器中存储有可在处理器上运行的计算机程序，所述处理器在执行所述计算机程序时实现上述方法实施例中的方法。

也就是说，以上方法实施例中的方法应理解为可由计算机程序指令实现地壳三维结构模型融合方法的流程。可提供这些计算机程序指令到处理器，使得通过处理器执行这些指令产生用于实现上述方法流程所指定的功能。

本实施例所指的处理器是指微处理器MCU或可编程逻辑器件FPGA等的处理装置。

本实施例所指的存储器包括用于存储信息的物理装置，通常是将信息数字化后再以利用电、磁或者光学等方式的媒体加以存储。例如：利用电能方式存储信息的各式存储器，RAM、ROM等；利用磁能方式存储信息的的各式存储器，硬盘、软盘、磁带、磁芯存储器、磁泡存储器、U盘；利用光学方式存储信息的各式存储器，CD或DVD。当然，还有其他方式的存储器，例如量子存储器、石墨烯存储器等等。

通过上述存储器、处理器以及计算机程序构成的装置，在计算机中由处理器执行相应的程序指令来实现，处理器可以搭载各种操作系统，如windows操作系统、linux系统、android、iOS系统等。