CN111766627A - 基于模型分辨度的自适应平滑面波成像方法 - Google Patents

Abstract

本发明属于地球物理领域，公开基于模型分辨度的自适应平滑面波成像方法，在m*n个网格点的网格上，除了边界上的格点，对于某一路径i，网格中第j个点在G矩阵的某条路径下的敏感度参数值为S，其周围四个点分别为S1，S2，S3，S4，其中对应点位置的1改为P_SS1，其它四个点依次类推，其中Mid为当前路径下所有点的敏感度参数的中位数，构成新的敏感度差分矩阵L；通过计算每个格点附近经过的射线路径的数量，自适应的将每个格点都与相邻格点都赋予更为合理的平滑参数，避免了因射线路径分布不均，各格点间均使用相同平滑参数而带来的反演结果假象。

Description

基于模型分辨度的自适应平滑面波成像方法

技术领域

本发明属于地球物理领域，具体涉及基于模型分辨度的自适应平滑面波成像方法。

背景技术

目前，基于射线理论的面波层析成像方法被广泛应用于地球物理领域。速度模型的平滑方法对于基于射线理论的面波成像方法有着非常重要的作用。(韩复兴等，2019；Yuan et al.,2017)。由于射线理论基于高频近似假设，只能对射线路径上的格点具有敏感性，因此在射线密度有限，地震波波场不清楚的情况下，为满足射线追踪的理论要求，需要对速度模型进行一定的平滑处理。类似的平滑处理虽然会使成像结果出现一定程度的偏差，但也是继续进行相关研究的必要步骤。

(韩复兴等，2009)

为实现高频近似假设下的射线追踪，前人做了大量的研究工作(Lailly andSinoquet,1996；Gray,2000；Pacheco and Larner2000；Gao et al.,2017)为克服在使用射线理论的研究中出现数据有限和结构复杂的情况下，压制结果中出现的假象，即由于面波射线路径分布不均匀而反演中仍使用统一的平滑参数而导致的反演结果假象，为了解决这个问题，提出基于模型分辨度的自适应平滑面波成像方法。

发明内容

针对现有技术的不足，本发明的目的在于提供基于模型分辨度的自适应平滑面波成像方法，解决了由于面波射线路径分布不均匀而反演中仍使用统一的平滑参数而导致的反演结果假象的问题。

本发明的目的可以通过以下技术方案实现：

基于模型分辨度的自适应平滑面波成像方法，包括以下步骤：

在m*n个网格点的网格上，除了边界上的格点，对于某一路径i，网格中第j个点在G矩阵的某条路径下的敏感度参数值为S，其周围四个点分别为S1，S2，S3，S4，关于定义S与S1点的平滑关系为：

其中对应点位置的1改为P_SS1，其它四个点依次类推，其中Mid为当前路径下所有点的敏感度参数的中位数，构成新的敏感度差分矩阵L；

将L带入至目标函数中，得到最终目标函数：

C_D和C_M分别表示数据协方差矩阵和模型协方差矩阵；Δd＝d_obs-d_cal表示数据残差，d_obs和d_cal分别为观测数据和模型为m时的模拟数据；m_p为先验模型；ε表示阻尼参数；η为平滑参数；L为敏感度差分矩阵。

进一步地，所述阻尼参数ε的选取包括以下步骤：每次迭代前都计算

和

的比值作为阻尼参数ε，其中Δd＝d_obs-d_cal，Δm＝m-m_p，σ_d和σ_Δm分别表示数据误差和模型误差。

进一步地，所述数据误差作为先验信息给出，而模型误差通过：

其中，c是相速度，l为路径长度，t为走时，然后再求σ_Δm；

式(8)在每一次迭代后重新估算。

进一步地，所述平滑参数η选取包括以下步骤：选定阻尼参数后，在每一次迭代中再取不同的值作为平滑参数，进对应的走时残差和模型平滑度做归一化后画出L曲线，并选取L曲线拐点处的平滑参数作为最合适的η值。

进一步地，所述模型误差中的参考模型为观测数据的相速度分布的平均值。

进一步地，所述参考模型求解问题时需要进行迭代，具体如下：

m_n+1＝m_n+δm_n， (9)

每次迭代后的模型改变量为

每次迭代后都把当前的计算模型作为下一次迭代的参考模型，并更新相应的射线路径和G矩阵，最后通过多次迭代获得最后需要的速度结构。

本发明的有益效果：

本发明通过计算每个格点附近经过的射线路径的数量，自适应的将每个格点都与相邻格点都赋予更为合理的平滑参数，避免了因射线路径分布不均，各格点间均使用相同平滑参数而带来的反演结果假象。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1是本发明实施例的在m*n的网格中，对某一路径i，第j个的点在G矩阵中的敏感度参数值为S分布示意图；

图2是本发明实施例的1s周期的路径分布图示意图；

图3是本发明实施例的1.4s周期的路径分布图示意图；

图4展示了1.4s周期确定阻尼参数后在第一次和最后一次迭代时选取平滑参数η使用的L曲线和选取的平滑参数数值示意图；

图5是本发明实施例的展示了1.4s周期最后一次迭代后G矩阵对应到每个网格点上的敏感度总和的分布示意图；

图6是本发明实施例的构造敏感度差分矩阵的方法所体现的格点间平滑关系的示意图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其它实施例，都属于本发明保护的范围。

通常高斯概率分布下的面波走时反演问题需要通过使用地震波走时场的信息，通过一定的优化方法估计地下介质的地球物理模型，使计算的地震波走时数据与观测数据达到最佳组合，即最小化一个目标函数

但在处理实际数据时，由于测量误差，噪声的存在以及测量数据的不完整，在反演过程中我们对式(1)表示的目标函数进行适当的修改，加入阻尼项和平滑项(即正则化)

式(4a)(4b)中，C_D和C_M分别表示数据协方差矩阵和模型协方差矩阵；Δd＝d_obs-d_cal表示数据残差，d_obs和d_cal分别为观测数据和模型为m时的模拟数据；m_p为先验模型；ε表示阻尼参数；η为平滑参数；L为敏感度差分矩阵。

为了在满足数据和抑制相邻速度结构间过大的变化中寻求平衡，通常加入的平滑项(Sambridge,1990；Tarantora,1987)

Ω(m)＝m^TD^T (3)

其中，Dm表示有限差分形式的空间离散矩阵。例如，在一个有梯度的界面上Dm可以表示为

其中，m_i表示界面网格上的每个点，i＝1,2,3…M，M表示网格点的总数。但这种方法只能使当前格点与下一格点有平滑效应。

本公开中，在做2D面波成像时，我们把Lm改为另外一种形式，即在一个有m*n个网格点的网格上，除了边界上的格点，我们让每个格点都与它相邻的4个格点上的值做平滑，并且把D矩阵中的正负1的数值关系用每个格点上的敏感度参数关系，即G矩阵中在该点位置的参数关系的来代替，构成新的敏感度差分矩阵。原目标函数变化为

C_D和C_M分别表示数据协方差矩阵和模型协方差矩阵；Δd＝d_obs-d_cal表示数据残差，d_obs和d_cal分别为观测数据和模型为m时的模拟数据；m_p为先验模型；ε表示阻尼参数；η为平滑参数；L为敏感度差分矩阵。

具体来说，对于某一路径i，网格中第j个点在G矩阵的某条路径下的敏感度参数值为S，其周围四个点分别为S1，S2，S3，S4，比如定义S与S1点的平滑关系为(如图1.)

即把L矩阵中对应点位置的1改为P_SS1，其它四个点一次类推，其中Mid为当前路径下所有点的敏感度参数的中位数，构成新的敏感度差分矩阵L。这样就使得在路径分布密集的地方有较少的平滑，在边界和路径分布稀疏的地方有较大的平滑。

反演中各参数的选取方法

为了在计算过程中尽量使数据项和正则项在量级上保持稳定的比例关系，在每次迭代时都根据当时的结果变化选取合适的阻尼参数ε和平滑参数η。具体做法是每次迭代前都计算

和

的比值作为阻尼参数ε，其中Δd＝d_obs-d_cal，Δm＝m-m_p，σ_d和σ_Δm分别表示数据误差和模型误差。这里的数据误差作为先验信息给出，而模型误差通过下面等式在每一次迭代后重新估算

其中，c是相速度，l为路径长度，t为走时，然后再求σ_Δm。

选定阻尼参数后，在每一次迭代中再取不同的值作为平滑参数，进对应的走时残差和模型平滑度做归一化后画出L曲线，并选取L曲线拐点处的平滑参数作为最合适的η值。

参考模型通常使用平均速度模型，即观测数据的相速度分布的平均值。

求解问题时需要进行迭代，

m_n+1＝m_n+δm_n， (9)

每次迭代后的模型改变量为

每次迭代后都把当前的计算模型作为下一次迭代的参考模型，并更新相应的射线路径和G矩阵，最后通过多次迭代获得最后需要的速度结构。至此，目标函数中的可变参数均可在计算过程中做到自适应的调节，以获得最优结果。

通过计算每个格点附近经过的射线路径的数量，自适应的将每个格点都与相邻格点都赋予更为合理的平滑参数，避免了因射线路径分布不均，各格点间均使用相同平滑参数而带来的反演结果假象。

试验验证：

为说明该方法在实例中的应用，本文使用了台北盆地的频散数据(Benson etal.,2007；Huang et al.，2010)来计算当地的速度结构。该频散数据从0.5s至3s，间隔为0.1s。使用Yao et al.，(2006)和Huang et al.(2010)的从背景噪音互相关中提取经验格林函数来测算Rayleigh波相速度的方法获得。通过基于快速行进的面波射线追踪方法获得各周期的路径分布图，如图2和图3，其中图示黑色三角形表示台站，通过路径的弯曲程度可以看出地下速度结构比较复杂。

本文在计算时将区域在南北向分为15个网格点，东西向19个网格点，网格间距为0.02；参考模型和初始模型都使用平均模型，数据协方差矩阵C_D中的值都设置为0.01；平滑参数η的从10^-4开始，每取值一次往上乘以5，一共乘20次，即一共取值21个；迭代一般在20-30次可以达到收敛。速度结构成像的结果主要呈现出中部低速，周围高速的特征，这与台北盆地中部有较厚的沉积层速度较低，周围为山脉为高速体的地形特征也比较吻合。

为了说明平滑参数η的选取过程，我们做了L曲线的图像。图4展示了1.4s周期确定阻尼参数后在第一次和最后一次迭代时选取平滑参数η使用的L曲线和选取的平滑参数数值。图中彩色的圆圈的颜色代表当前平滑参数η的大小，位置代表当前η值对应的归一化后的走时残差和模型平滑度的关系，红色星号代表使用L曲线后所选取的最合适的平滑参数。

图4通过L曲线选取平滑参数η。圆圈颜色对应当前平滑参数的大小，黑色星号代表选取的拐点处的平滑参数。

图5展示了1.4s周期最后一次迭代后G矩阵对应到每个网格点上的敏感度总和的分布(图5)以及按照公开中所述的构造敏感度差分矩阵的方法所体现的格点间平滑关系(图6)。上图中的横纵坐标分别表示E-W和N-S方向上的15和19个格点，灰度表示当前格点处的敏感度大小。

我们可以看到该周期在南北方向和东西方向上都在第10-12格点处的敏感度最强，其所代表的分辨率也应该最好。下图中的黑色星号表示在反演过程中所选取的南北向15个格点，东西向19个格点。从每个格点向其四个方向延伸出的黑色线段表示了当前各点与其周围四个格点的平滑关系。黑色线段越长，该点与黑色线段指向的格点间的平滑关系越剧烈。

对比图5，我们可以看出在敏感度较高的地方，相应格点间的平滑也就越小；对周围路径没有分辨率的地方，平滑也就达到最大值。在分辨率由高到低的边缘，我们也能看到递减的平滑关系。计算过程中，我们没有对网格最外围的点处的速度值做平滑。

1.4s周期最后一次迭代后G矩阵对应到每个网格点上的敏感度总和的分布(图5)，横纵坐标分别表示E-W和N-S方向上的格点数，图像的灰度值表示当前格点在G矩阵中的敏感度总和值，灰度越大，敏感度越高。

格点间平滑关系(图6)，黑色星号表示在反演过程中在经纬度上所选取的格点，从每个格点向其四个方向延伸出的黑色线段表示了当前各点与其周围四个格点的平滑关系大小，线段越长，表示当前格点与线段所指向方向的格点间平滑越大。

在本说明书的描述中，参考术语“一个实施例”、“示例”、“具体示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不一定指的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任何的一个或多个实施例或示例中以合适的方式结合。

以上显示和描述了本发明的基本原理、主要特征和本发明的优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。

Claims (6)

1.基于模型分辨度的自适应平滑面波成像方法，其特征在于，包括以下步骤：

在m*n个网格点的网格上，除了边界上的格点，对于某一路径i，网格中第j个点在G矩阵的某条路径下的敏感度参数值为S，其周围四个点分别为S1，S2，S3，S4，关于定义S与S1点的平滑关系为：

其中对应点位置的1改为P_SS1，其它四个点依次类推，其中Mid为当前路径下所有点的敏感度参数的中位数，构成新的敏感度差分矩阵L；

将L带入至目标函数中，得到最终目标函数：

C_D和C_M分别表示数据协方差矩阵和模型协方差矩阵；Δd＝d_obs-d_cal表示数据残差，d_obs和d_cal分别为观测数据和模型为m时的模拟数据；m_p为先验模型；ε表示阻尼参数；η为平滑参数；L为敏感度差分矩阵。

2.根据权利要求1所述的基于模型分辨度的自适应平滑面波成像方法，其特征在于，所述阻尼参数ε的选取包括以下步骤：每次迭代前都计算

和

的比值作为阻尼参数ε，其中Δd＝m_obs-d_cal，Δm＝m-m_p，σ_d和σ_Δm分别表示数据误差和模型误差。

3.根据权利要求2所述的基于模型分辨度的自适应平滑面波成像方法，其特征在于，所述数据误差作为先验信息给出，而模型误差通过：

其中，c是相速度，l为路径长度，t为走时，然后再求σ_Δm；

式(8)在每一次迭代后重新估算。

4.根据权利要求3所述的基于模型分辨度的自适应平滑面波成像方法，其特征在于，所述平滑参数η选取包括以下步骤：选定阻尼参数后，在每一次迭代中再取不同的值作为平滑参数，进对应的走时残差和模型平滑度做归一化后画出L曲线，并选取L曲线拐点处的平滑参数作为最合适的η值。

5.根据权利要求4所述的基于模型分辨度的自适应平滑面波成像方法，其特征在于，所述模型误差中的参考模型为观测数据的相速度分布的平均值。

6.根据权利要求5所述的基于模型分辨度的自适应平滑面波成像方法，其特征在于，所述参考模型求解问题时需要进行迭代，具体如下：

m_n+1＝m_n+δm_n， (9)

每次迭代后的模型改变量为

每次迭代后都把当前的计算模型作为下一次迭代的参考模型，并更新相应的射线路径和G矩阵，最后通过多次迭代获得最后需要的速度结构。

Images