CN115358311B - 地表变形监测多源数据融合处理方法 - Google Patents

Abstract

本发明公开了一种地表变形监测多源数据融合处理方法，涉及数据处理方法技术领域。所述方法包括如下步骤：利用测量获取的高精度点状监测数据，基于近邻距离加权平均法的同名点InSAS监测数据提取；在工点区域监测重合点位上进行数据融合计算；利用融合后数据对其他点位监测值进行改正计算；获取工点区域高精度面状竖向沉降及水平位移信息。所述方法能够获取高精度、高时间分辨率、高空间分辨率的地表变形监测成果，准确反映目标工点或大区域变形现状和规律，为防灾减灾决策提供数据支撑。

Description

地表变形监测多源数据融合处理方法

技术领域

本发明涉及数据处理方法技术领域，尤其涉及一种地表变形监测多源数据融合处理方法。

背景技术

地表变形是地表发生高程下降及水平位移的一种地质复合运动现象，一般是由自然因素和人为因素的共同作用导致的。目前地表变形监测的方法很多，各有优劣，针对较大范围地表变形监测的技术方法主要有全站仪或水准测量、GNSS测量和合成孔径雷达干涉(InSAR)技术等。这三种方法在测量精度、时间分辨率和空间分辨率各不相同，全站仪或精密水准测量技术测量精度最高，GNSS连续测量技术时间分辨率最高，InSAR测量技术空间分辨率最高。同时水准与GNSS测量方法仅能监测点状空间的地表沉降及水平位移信息，且需花费大量人力、物力，InSAR测量技术则可以覆盖全区域的地表变形信息。

尤其近年随着GNSS、InSAR等技术的快速进步发展，将这些方法同时应用于地表变形监测领域已逐步普及。当目标对象工点或大区域利用了三种方法同时展开地表变形监测，不同方法获取的监测数据在精度和时空上均不一致，会造成分析结果有差异，甚至判断失误风险。考虑到当前地表变形监测领域多集中在上述单种方法的工程应用上，对多源监测数据综合利用研究应用较少。

发明内容

本发明所要解决的技术问题是如何提供一种能够获取高精度、高时间分辨率、高空间分辨率的地表变形监测成果，准确反映目标工点或大区域变形现状和规律，为防灾减灾决策提供数据支撑的地表变形监测多源数据融合处理方法。

为解决上述技术问题，本发明所采取的技术方案是：一种地表变形监测多源数据融合处理方法，其特征在于包括如下步骤：

利用测量获取的高精度点状监测数据，基于近邻距离加权平均法的同名点InSAS监测数据提取；

在工点区域监测重合点位上进行数据融合计算；

利用融合后数据对其他点位监测值进行改正计算；

获取工点区域高精度面状竖向沉降及水平位移信息。

进一步的技术方案在于，所述基于近邻距离加权平均法的同名点InSAS监测数据提取的具体步骤如下:

采用近邻距离加权平均法进行同名点提取，选择与全站仪测点或水准点位置相距一定距离d以内的邻近点，将其加权平均值视为该全站仪测点或水准点、GNSS点对应的InSAS结果；在估计提取同名点的InSAS监测值时,考虑离全站仪测点或水准点最近的N个InSAS点,这N个InSAS点对同名点的影响与距离成反比；二维场景下，设某全站仪或水准及GNSS重合点为s(x₀,y₀)到周围InSAS点(x_i,y_i)距离D_ik，则：

求出距离全站仪或水准及GNSS重合点s最近的N个InSAS点的距离D_ik(i＝1,2,…,N)，则该目标范围的估计值为：

其中，Z_i为InSAS点i上的测量值,γ为提取同名点InSAS监测估计值；

若以水准、GPS观测成果为地表变形的最或然值，选择标准差作为衡量精度的指标，求解公式如下：

其中，y为水准监测地表变形值；

为InSAR监测地表变形值。

进一步的技术方案在于，所述重合点位上进行数据融合计算的方法具体包括如下步骤：

设有n种监测方法来监测某位置的地表变形向量为Y，设X为待估参数向量，则有观测方程：

Y＝KX+ε

其中，K为观测矩阵，ε为监测数据噪声；

利用最小二乘法对观测方程参数向量X进行最优估计，则有：

其中，V为残差，P为加权权阵；

对上式估值

求一阶偏导，并令一阶偏导为零，可得最优估计：

由于二阶导数：

所以，估值

时，存在满足方程达到最小；

参数X估值误差：

进一步的技术方案在于，所述对其他点位监测值进行改正计算的方法包括如下步骤：

设有m个已知点位O(x_i,y_i)融合后的地表变形估值，可用n个核函数的总和去拟合或逼近其他任意点位的地表变形值

若待估点位观测变形量为Z_i，则该点位的检测值改正数或残差为：/>

以区域内待求点的残差或改正数为拟合对象，由多面函数的一般式可建立任意点位的地表变形改正模型为：

式中，α_i为模型待估计参数；n为选定的已知变形量的点，即结点个数；Q(x,y；x_i,y_i)为核函数，核函数的选择包括正双曲面函数、倒双曲面函数和三次曲面函数；

①正双面曲函数

Q(x,y；x_i,y_i)＝[(x-x_i)²+(y-y_i)²+δ²]^1/2

②倒双面曲函数

Q(x,y；x_i,y_i)＝[(x-x_i)²+(y-y_i)²+δ²]^-1/2

③三次曲面函数

Q(x,y；x_i,y_i)＝[(x-x_i)²+(y-y_i)²]^3/2++δ²

式中，δ为光滑因子；k为指数因子，决定了核函数的形状。

平滑因子δ的作用是改变核函数的形状，δ值越大，核函数所表达的曲面越平缓，反之曲面越陡峭，不同的核函数对δ的变化敏感程度是不同；

此处，残差方程矩阵形式为：

当m＞n时，依据最小二乘法原则，并按等权求解得待估参数最优估值：

将求得的系数

代入到多面函数方程中，可获得多面函数方程：

再利用点位的坐标，求得任意点位的改正数或残差，利用求得的改正数或残差对点位监测值进行改正，从而获取点位的拟合估值为：

采用上述技术方案所产生的有益效果在于：本申请通过多种方法对地表变形监测数据的融合分析处理，保证了数据的正确性，提高了数据测量精度，及监测点位的时空分辨率，可以通过有效地数据组合方式来解决监测中的费用过高的问题，能够利用低成本地多种监测方法进行监测，采用数据融合技术达到既满足监测要求，提高监测效率，又降低费用的要求，为地表变形变化规律的分析，提供更可靠的技术支持。

附图说明

下面结合附图和具体实施方式对本发明作进一步详细的说明。

图1是本发明实施例所述方法的主流程图；

图2是本发明实施例所述方法中地表变形修正计算流程图。

具体实施方式

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

在下面的描述中阐述了很多具体细节以便于充分理解本发明，但是本发明还可以采用其他不同于在此描述的其它方式来实施，本领域技术人员可以在不违背本发明内涵的情况下做类似推广，因此本发明不受下面公开的具体实施例的限制。

如图1所示，本发明实施例公开了一种地表变形监测多源数据融合处理方法，包括如下步骤：

首先，利用全站仪或精密水准测量和GNSS测量获取的高精度点状监测数据，基于近邻距离加权平均法的同名点InSAS监测数据提取；

其次，在工点区域监测重合点位上进行数据融合计算；

最后，利用融合后数据对其他点位监测值进行改正计算。

通过上述步骤，即可获取工点区域高精度面状竖向沉降及水平位移信息，最后为预测分析提供可靠有效数据支持。

下面结合具体方法对上述步骤进行说明：

一、基于近邻距离加权平均法的同名点InSAS监测数据提取

根据全站仪测点或水准点、GNSS点坐标位置提取对应的InSAS沉降及水平位移变形速率值。由于InSAS测量所获取的是地表相关目标的变形值，而全站仪或水准、GNSS测量获取的是地表特定点位置的变形值，二者并非同一位置，需要根据一定的准则提取全站仪测点或水准点、GNSS点对应位置的InSAS监测结果。

此处采用近邻距离加权平均法进行同名点提取，选择与全站仪测点或水准点位置相距一定距离d以内的邻近点，将其加权平均值视为该全站仪测点或水准点、GNSS点对应的InSAS结果。在估计提取同名点的InSAS监测值时,考虑离全站仪测点或水准点最近的N个InSAS点,这N个InSAS点对同名点的影响与距离成反比。二维场景下，不妨假设某全站仪或水准及GNSS重合点为s(x₀,y₀)到周围InSAS点(x_i,y_i)距离D_ik，则：

求出距离全站仪或水准及GNSS重合点s最近的N个InSAS点的距离D_ik(i＝1,2,…,N)，则该目标范围的估计值为：

其中，Z_i为InSAS点i上的测量值,γ为提取同名点InSAS监测估计值。

若以水准、GPS观测成果为地表变形的最或然值，选择标准差作为衡量精度的指标，求解公式如下：

其中，y为水准监测地表变形值；

为InSAR监测地表变形值。

二、重合点位数据融合计算

设有n种监测方法来监测某位置的地表变形向量为Y，设X为待估参数向量，则有观测方程：

Y＝KX+ε

其中，K为观测矩阵，ε为监测数据噪声。

利用最小二乘法对观测方程参数向量X进行最优估计，则有：

其中，V为残差，P为加权权阵。

对上式估值

求一阶偏导，并令一阶偏导为零，可得最优估计：

由于二阶导数：

所以，估值

时，存在满足方程达到最小。

参数X估值误差：

三、基于改进多面函数拟合法的其他点位监测值改正算法

考虑到地表变形监测场景包括竖向沉降及水平位移变形内容，本申请提出了一种多面函数拟合地表变形的改进算法，也就是对点位竖向沉降及水平位移分别采用不同的核函数进行建模分析计算的思想。具体来讲，即通过变换光滑因子和指数因子，利用量值残差比较，分别求取两者最佳拟合核函数，从而更加精确有效地建立工点区域地表变形模型，有效改善模型拟合效果，提高建模精度。

设有m个已知点位O(x_i,y_i)融合后的地表变形估值，可用n个核函数的总和去拟合或逼近其他任意点位的地表变形值

若待估点位观测变形量为Z_i，则该点位的检测值改正数(或残差)为：/>

此处以区域内待求点的残差或改正数为拟合对象，由多面函数的一般式可建立任意点位的地表变形改正模型为：

式中，α_i为模型待估计参数；n为选定的已知变形量的点，即结点个数；Q(x,y；x_i,y_i)为核函数，核函数的选择应结合实际情况进行，不同形式的核函数有不同的拟合效果，目前效果较好的核函数有正双曲面、倒双曲面、三次曲面等。

①正双面曲函数

Q(x,y；x_i,y_i)＝[(x-x_i)²+(y-y_i)²+δ²]^1/2

②倒双面曲函数

Q(x,y；x_i,y_i)＝[(x-x_i)²+(y-y_i)²+δ²]^-1/2

③三次曲面函数

Q(x,y；x_i,y_i)＝[(x-x_i)²+(y-y_i)²]^3/2++δ²

式中，δ为光滑因子；k为指数因子，决定了核函数的形状。

多面函数拟合法的关键是确定平差结点、指数因子和光滑因子。选择平差结点的原则是尽量分布均匀，覆盖整个研究区域，防止出现过拟合现象。平滑因子δ的作用是改变核函数的形状，δ值越大，核函数所表达的曲面越平缓，反之曲面越陡峭。不同的核函数对δ的变化敏感程度是不同。

此处，残差方程矩阵形式为：

当m＞n时，依据最小二乘法原则，并按等权求解得待估参数最优估值：

将求得的系数

代入到多面函数方程中，可获得多面函数方程：

再利用点位的坐标，求得任意点位的改正数(或残差)，利用求得的改正数(或残差)对点位监测值进行改正，从而获取点位的拟合估值为：

综上，本发明通过多种方法对地表变形监测数据的融合分析处理，保证了数据的正确性，提高了数据测量精度，及监测点位的时空分辨率，可以通过有效地数据组合方式来解决监测中的费用过高的问题，甚至利用低成本地多种监测方法进行监测，采用数据融合技术达到既满足监测要求，提高监测效率，又降低费用的要求，为地表变形变化规律的分析，提供更可靠的技术支持。

本方法适于水利水电、公路铁路、市政及新能源等多行业工程地面或地表沉降与水平变形病害或灾害安全防治领域，为工程防灾减灾决策提供可靠技术支持。

Claims (4)

1.一种地表变形监测多源数据融合处理方法，其特征在于包括如下步骤：

利用测量获取的高精度点状监测数据，基于近邻距离加权平均法的同名点InSAR监测数据提取；

在工点区域监测重合点位上进行数据融合计算；

利用融合后数据对其他点位监测值进行改正计算；

获取工点区域高精度面状竖向沉降及水平位移信息；

所述基于近邻距离加权平均法的同名点InSAR监测数据提取的具体步骤如下：

采用近邻距离加权平均法进行同名点提取，选择与全站仪测点或水准点位置相距一定距离d以内的邻近点，将其加权平均值视为该全站仪测点或水准点、GNSS点对应的InSAR结果；在估计提取同名点的InSAR监测值时，考虑离全站仪测点或水准点最近的N个InSAR点，这N个InSAR点对同名点的影响与距离成反比；二维场景下，设某全站仪或水准及GNSS重合点为s(x₀，y₀)到周围InSAR点(x_i，y_i)距离D_ik，则：

求出距离全站仪或水准及GNSS重合点s最近的N个InSAR点的距离D_ik(i＝1，2，...，N)，则目标范围的估计值为：

其中，Z_i为InSAR点i上的测量值，γ为提取同名点InSAR监测估计值；

若以水准、GNSS观测成果为地表变形的最或然值，选择标准差作为衡量精度的指标，求解公式如下：

其中，y为水准监测地表变形值；

为InSAR监测地表变形值。

2.如权利要求1所述的地表变形监测多源数据融合处理方法，其特征在于：利用全站仪或精密水准测量和GNSS测量获取高精度点状监测数据。

3.如权利要求1所述的地表变形监测多源数据融合处理方法，其特征在于，所述重合点位上进行数据融合计算的方法具体包括如下步骤：

设有n种监测方法来监测某位置的地表变形向量为Y，设X为待估参数向量，则有观测方程：

Y＝KX+ε

其中，K为观测矩阵，ε为监测数据噪声；

利用最小二乘法对观测方程待估参数向量X进行最优估计，则有：

其中，V为残差，P为加权权阵；

对上式估值

求一阶偏导，并令一阶偏导为零，可得最优估计：

由于二阶导数：

所以，估值

时，存在满足方程达到最小；

待估参数向量X估值误差：

4.如权利要求1所述的地表变形监测多源数据融合处理方法，其特征在于，所述对其他点位监测值进行改正计算的方法包括如下步骤：

设有m个已知点位O(x_i，y_i)融合后的地表变形估值，可用n个核函数的总和去拟合或逼近其他任意点位的地表变形值

若待估点位观测变形量为Z_i，则该点位的监测值改正数或残差为：/>

以区域内待估点位的残差或改正数为拟合对象，由多面函数的一般式可建立任意点位的地表变形改正模型为：

式中，α_i为模型待估参数；n为选定的已知变形量的点，即结点个数；Q(x，y；x_i，y_i)为核函数，核函数的选择包括正双曲面函数、倒双曲面函数和三次曲面函数；

①正双面曲函数

Q(x，y；x_i，y_i)＝[(x-x_i)²+(y-y_i)²+δ²]^1/2

②倒双面曲函数

Q(x，y；x_i，y_i)＝[(x-x_i)²+(y-y_i)²+δ²]^-1/2

③三次曲面函数

Q(x，y；x_i，y_i)＝[(x-x_i)²+(y-y_i)²]^3/2++δ²

式中，δ为光滑因子；k为指数因子，决定了核函数的形状；

平滑因子δ的作用是改变核函数的形状，δ值越大，核函数所表达的曲面越平缓，反之曲面越陡峭，不同的核函数对δ的变化敏感程度是不同；

此处，残差方程矩阵形式为：

当m＞n时，依据最小二乘法原则，并按等权求解得待估参数最优估值：

将求得的系数

代入到多面函数方程中，可获得多面函数方程：

再利用点位的坐标，求得任意点位的改正数或残差，利用求得的改正数或残差对点位的监测值进行改正，从而获取点位的拟合估值为：

Images