CN111798135A - 一种基于多源数据集成的高铁沉降危害性评估方法 - Google Patents

Abstract

本发明提供了一种基于多源数据集成的高铁沉降危害性评估方法，包括以下步骤：S1、利用InSAR技术采集地面形变数据，包括采集InSAR影像数据、光学遥感影像数据、水文地质资料数据、二级水准测量数据；S2、对步骤S1中采集的数据进行预处理；S3、采用PS‑InSAR+SBAS‑InSAR两种时序InSAR技术分析提取地面形变结果；S4、针对高铁线路背景区域沉降分离；S5、对影响地面形变的危害因子采用MIC算法赋权重；S6、地面沉降危害性评估。本发明所述的基于多源数据集成的高铁沉降危害性评估方法评估结果趋于大图幅、快速、客观、经济高效的地面沉降危害性评估体系，做到长效监测，保障城市的安全发展，减少其对生命财产经济安全的威胁。

Description

一种基于多源数据集成的高铁沉降危害性评估方法

技术领域

本发明属于地面沉降评估技术领域，尤其是涉及一种基于多源数据集成的高铁沉降危害性评估方法。

背景技术

地面沉降是一种普遍存在的一种不可逆的地质灾害，且对城市建筑群和连接城市之间的交通线路存在很大的安全隐患。随着我国城市化的进展，城区面积在不断扩大，且高速铁路承担着城际之间主要的人员流动，对地面的形变要求就更加严格。地面沉降对城市化的危害也不断增加，且地面沉降对于城市建筑物的影响是缓慢的发育过程，因此需要对其进行长期监测，并且对于地面沉降的进行危害性评价，建立长效的监测评价机制。InSAR地面形变监测已经广泛应用在地面形变监测上，相较于GPS控制网和水准控制网相比，可得到大幅图的形变监测结果，且结果形式多样，目前InSAR对地监测的结果可靠性已经很高，但目前对于地面监测的危害性评估多基于水准监测和以往的观测记录，时效性较差，且监测结果为离散数据点，在危害性评价上，得到的数据多为离散数据点，不能直观反映大幅图连续的地面沉降危害，且在进行危害性评估时，对于各个危害因子的权重赋值多基于主观的经验值，在不同的地区地形会无法形成统一有效的评估体系，且地域地区地形不同，主观经验值借鉴意义不大。

发明内容

有鉴于此，为克服上述缺陷，本发明旨在提出一种集成InSAR对地观测、光学遥感、水文地质资料、机器学习等多源数据多处理方法的地面沉降危害性评估体系，针对危害因子的赋值，引入统计学的最大信息系数法(Maximum Information Coefficient)，建立更客观的赋值方法，结合遥感技术，使得评估结果趋于大图幅、快速、客观、经济高效的地面沉降危害性评估体系，做到长效监测，保障城市的安全发展，减少其对生命财产经济安全的威胁。

为达到上述目的，本发明的技术方案是这样实现的：

一种基于多源数据集成的高铁沉降危害性评估方法，包括以下步骤：

S1、利用InSAR技术采集地面形变数据，包括采集InSAR影像数据、光学遥感影像数据、水文地质资料数据、二级水准测量数据；

S2、对步骤S1中采集的数据进行预处理；

S3、采用PS-InSAR+SBAS-InSAR两种时序InSAR技术分析提取地面形变结果；

S4、针对高铁线路背景区域沉降分离；

S5、对影响地面形变的危害因子采用MIC算法赋权重；

S6、地面沉降危害性评估。

进一步的，所述步骤S2中，对采集的数据进行预处理的方法包括：

对光学遥感影像数据的辐射定标、大气校正、几何校正；

对水文地质资料数据的插值化、矢量化栅格化处理；

对InSAR影像数据的去噪、多视、选取主影像。

进一步的，所述步骤S3中，针对有大量为稳定散射体的高铁沿线的处理方法为：进行PS-InSAR流程处理，在多时序SAR影像中选取时间空间基线小的影像为主影像，通过相干系数阈值和振幅离差阈值两种方法对影像进行选点处理，选取稳定的PS点；

针对有大量荒地或逐步建成的区域的处理方法为：进行SBAS-InSAR流程处理，除去必要的去噪之外，不必进行主影像选取，在影像基线优化处理时，引入最小二乘法优化干涉对组，在保障影像处理质量的同时，减少冗余干涉对，减轻计算负担。

进一步的，进行PS-InSAR流程处理的具体方法如下：

将InSAR影像数据进行PS-InSAR处理，需要在N+1景影像中选取1景主影像m，其余N景InSAR影像将与主影像匹配，生成N幅干涉图，基于PS处理流程的特殊性，一般要求N≧20，且主影像的时间空间基线在N+1景影像中最小，在选取主影像时，要建立相关函数模型：

上述公式(1)、(2)中，

是SAR影像N和M的空间基线，T^n,m为影像的时间基线，

为多普勒基线，B_c，T_c，和f_c为临界基线，γ^m max的SAR影像为主影像，则由式(1)可知，当空间基线和时间基线值最小时，γ^m值最大，因此主影像的时间空间基线应是多景SAR影像中最小的；

对于地形相位和平地相位，PS-InSAR流程中是采取加入DEM来减少影响，在干涉相位有5部分组成，在DEM校正中上述两种相位将被去除，其中高程上由于DEM选用的精度会存在误差，因此对于干涉图N的像元，其相位为：

式(6)中，

为参考DEM校正中由于DEM精度引起的残余地形相位，

为大气相位，

为形变相位，

为其余噪声，由式(3)可知，DEM精度越高则残余地形相位会偏小，则有：

将(4)，(5)中的

带入(3)可知：

对于(6)中的大气相位和噪声等非线性相位，在覆盖大量混凝土的城市、道路或高速铁路稳定地表的区域，二者的影像可视为线性形变；

经过上述步骤，会产生永久散射点即PS点，除了有形变信息的PS点外，还有许多对形变观测无用的点，所以需对PS点进行筛选；在本步骤中采取相干系数阈值法和振幅离差阈值法对PS点进行筛选，选取相干性高的PS点作为PS点，其中相干系数阈值法以一个局部矩阵进行像元的相干性计算，具体为：

式(7)中，γ是相干系数，M，S是两景SAR影像，m，n是局部矩阵的行列数，γ计算结果值为(0,1)，相干系数阈值是在进行运算之前给定阈值，大于阈值选为PS点，且相干系数值越接近于1表示其相干性越强；

相反的，振幅离差阈值是针对影像的单个像元，给定一个阈值范围，以像元振幅的稳定性来选取PS点，具体表现：

式(8)中，σ_φ和D_A表示相位标准差和振幅离差指数，μ_A为均值振幅，σ_A是振幅标准差，σ_nR为噪声实部标准差，σ_nl为噪声虚部标准差；

进行选点之后就可以进行解算，其中包括去除大气延迟、轨道误差等，然后进行线性形变估算，最终得到PS点的形变速率，对其进行地理编码即可得到形变结果。

进一步的，进行SBAS-InSAR流程处理的具体方法如下：

不同于PS-InSAR流程处理，SBAS-InSAR流程处理无需对N景SAR影像选取主影像，在用于生成干涉图的SAR影像时选取时间空间基线较小的SAR影像形成干涉组对，以减少失相干对结果的影响，N景SAR影像生成M幅干涉图，干涉图的相位表示：

式(9)中，x,r是干涉图坐标，t_A,t_B为生成干涉图的两景影像的获取时间，若t_A<t_B，式(9)中第一项即为形变相位，式中的d(t_B,x,r),d(t_A,x,r)为t_A,t_B时刻(x,r)相对于t0时刻视线方向的累积形变量；第二项即为残余高程相位；第三项为大气延迟相位；余式为噪声相位；若假设M景主影像和从影像的获取时间以序列IE＝[IE₁...IE_M]与IS＝[IS₁...IS_M],若对M幅干涉图先不考虑大气延迟，只求解高程误差和形变量，则：

且将形变相位转化为平均相位速度矢量[v]^T：

若速度矢量是线性的，即v＝Mp，则有矩阵方程：

DMp+CΔz＝Δφ (12)

式(12)中D是M×(N-1)的矩阵，p为速度模型的参数矢量，则参数个数i应小于等于(N-1)，M为(N-1)×i的矩阵；Δz为高程相位误差；对于式(12)采用最小二乘法求解，可得到形变相位和高程误差相位，并可以去除相位图中高程相位误差，则矩阵方程可修改为：

Dv+CΔz＝Δφ (13)

在SBAS处理中，由于考虑到空间失相干的问题，将影像分成若干基线对，导致式(13)中的D为秩亏矩阵，若基线对为L(L>1)，则D的秩为N-L，使用SVD法对式(13)求解，对得到的速度矢量Dv积分即可获得形变相位序列，需要注意的是此时结果中还有大气相位，一般采取高通滤波去除大气相位，得到研究区的形变序列,即研究区的地面形变信息。

进一步的，所述步骤S4的具体方法如下：

S401、对步骤S3中采集的地面形变数据结果进行校准，针对高速铁路采用PS-InSAR流程处理，将结果与二等水准测量数据进行校正；

S402、对形变区域进行高频低频滤波的小波变换分析，用于反映出沉降区域的细节部分与趋势部分的差异；

其中，通过高频滤波组的用于表达高铁线路及周围沉降区域细节性的信息；通过低频滤波的用于反映整区域趋势性相关信息。

进一步的，所述步骤S402中，小波变换分析的具体方法如下：

利用小波变换完成一次变换，将影像数据分成高频图组与低频图组，其中高频组会更好的反映锐化的部分；

对于一维的信号：

定义

其逆变换为：

上式中Ψ为傅里叶变换，C_Ψ取值为：

对于栅格影像的二维的小波变换，有一个前提条件是尺度可分离的，即信号函数Φ(x,y)＝Φ(x)Φ(y),其中Φ(x)Φ(y)分别对应的小波函数为Ψ(x)Ψ(y)其创立二维小波变换的三个基本小波可表达为：

ψ¹(x,y)＝φ(x)ψ(y),ψ²(x,y)＝φ(y)ψ(x),ψ³(x,y)＝ψ(x)ψ(y) (16)

式中构建的二维平方可积函数空间L²(R²)的正交归一，既有：

l＝1,2,3,...且j,l,m,n均为整数；

对于一幅N×N的影像f(x,y)开始，其中上标尺度N，N＝2n，对于j＝0，尺度2j＝20＝1，为源影像尺度；j每增大一次使尺度加倍，分辨率减半，影像分为三个高频组，一个低频组，低频组中包括大部分能量，还需进行分离，高频组则更好的反映锐化等凸显线条的影像；对于第一个层次(j＝1)，则有：

通过小波变换分解后重构，将InSAR影响处理结果中有关高速铁路的部分提取出来，高频组中将高铁及沿线研究区的信息更加突出的表现。

进一步的，所述步骤S5的具体方法如下：

S501、针对步骤S1中采集的数据，选择影响高铁行车地面形变的危害因子，包括地下水位、土壤压缩模量、交通流量、压缩层厚度；

S502、利用GIS软件工具，将数据进行统一的栅格化处理，对于数据点，首选克里金插值法进行插值，对于遥感地面沉降数据、光学遥感数据栅格数据，直接针对栅格灰度值进行栅格计算；

S503、将危害因子统一转为栅格数据后，根据MIC算法，对栅格灰度值之间进行计算，从而完成危害因子赋权重。

进一步的，所述步骤S6的具体方法如下：

S601、基于步骤S5中得到的赋权重的危害因子与地面沉降图并结合随机森林算法建立决策树，并且将处理后的光学遥感影像数据进行分类，用以对高速铁路沿线不同地段的危害性进行评估；

S602、结合Boruta算法辅助修改各级决策树，Boruta算法通过对决策树建立副本并对比新建随机副本与原始变量的重要性对决策树进行修改，用于对高铁沿线行车安全地段进行提取标记，此方法可对赋值的栅格影像进行计算，根据权重提取特征地段，实现了对高铁沿线沉降危害性的分区与评估。

进一步的，在执行步骤S601前，在建立决策树前，需对影像建立训练样本，通过构建训练样本以确定特征属性和目标变量的关系；采用CART算法构建决策树训练样本，CART算法通过基尼系数确定属性节点分割属性，其中选取属性基尼系数最小的进行二分化，最后构建一个二叉树的决策树：对于一个栅格数据集D，它的基尼系数为：

式中p_i是数据集D中任意数据被抽取成为训练样本的概率，对生成的二分数据集进行二分则有：

式中D₁，D₂为二分后的两个数据集，D₁为可进行二分的子数据集，并对D₁进行式22的计算再进行二分，并在节点处形成训练规则，直到所有二分子集都符合节点的规则时，结束算法，若有不符合的，可对其再次进行式22的计算，直到结束算法。

相对于现有技术，本发明所述的基于多源数据集成的高铁沉降危害性评估方法具有以下优势：

(1)本发明所述的基于多源数据集成的高铁沉降危害性评估方法将统计学的最大信息系数法(MIC)算法引入危害因子的权重赋值，以取代以往的基于专家先验知识的主观赋值，并将此算法引入影像计算，在影像计算中完成对危害性因子的评估，而不是对表格数据进行数字计算，相较于表格数字计算，直接对栅格影像进行影像计算，计算结果以更直观的，大图幅的栅格遥感影像展现，更便于直观上的判断；加之影像的空间属性，每一个栅格像素都可以进行定位分析，以便于在城市规划等领域的应用。

附图说明

构成本发明的一部分的附图用来提供对本发明的进一步理解，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。在附图中：

图1为本发明实施例所述的基于多源数据集成的高铁沉降危害性评估方法流程图。

具体实施方式

需要说明的是，在不冲突的情况下，本发明中的实施例及实施例中的特征可以相互组合。

在本发明的描述中，需要理解的是，术语“中心”、“纵向”、“横向”、“上”、“下”、“前”、“后”、“左”、“右”、“竖直”、“水平”、“顶”、“底”、“内”、“外”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。此外，术语“第一”、“第二”等仅用于描述目的，而不能理解为指示或暗示相对重要性或者隐含指明所指示的技术特征的数量。由此，限定有“第一”、“第二”等的特征可以明示或者隐含地包括一个或者更多个该特征。在本发明的描述中，除非另有说明，“多个”的含义是两个或两个以上。

在本发明的描述中，需要说明的是，除非另有明确的规定和限定，术语“安装”、“相连”、“连接”应做广义理解，例如，可以是固定连接，也可以是可拆卸连接，或一体地连接；可以是机械连接，也可以是电连接；可以是直接相连，也可以通过中间媒介间接相连，可以是两个元件内部的连通。对于本领域的普通技术人员而言，可以通过具体情况理解上述术语在本发明中的具体含义。

下面将参考附图并结合实施例来详细说明本发明。

如图1所示，一种基于多源数据集成的高铁沉降危害性评估方法，包括以下步骤：

步骤一：InSAR技术地面形变信息采集

采用星载InSAR对研究区进行观测，得到大幅图的研究区SAR影像资料；同时收集研究区多光谱遥感影像，研究区二级水准控制网监测数据。

步骤二：针对数据的预处理。

包括光学遥感影像的辐射定标，大气校正、几何校正；

对于水文地质资料的插值化，矢量化栅格化处理；

InSAR影像的去噪，多视，选取主影像。

步骤三：地面形变数据的获取

本发明中，针对从SAR影像中提取地面形变信息主要基于采用PS-InSAR+SBAS-InSAR两种时序InSAR技术进行针对性的地面形变信息提取。

1、其中注意针对有大量为稳定散射体的高铁沿线，优先进行PS-InSAR流程处理，在多时序SAR影像中选取时间空间基线小的影像为主影像，通过相干系数阈值和振幅离差阈值两种方法对影像进行选点处理，选取稳定的PS点；

2、针对有大量荒地或逐步建成的区域，建议进行SBAS-InSAR流程处理，除去必要的去噪之外，不必进行主影像选取，在影像基线优化处理时，引入最小二乘法优化干涉对组，在保障影像处理质量的同时，减少冗余干涉对，减轻计算负担。

将SAR影像进行PS-InSAR处理，一般需要在N+1景影像中选取1景主影像m，其余N景SAR影像将与主影像匹配，生成N幅干涉图。基于PS处理流程的特殊性，一般要求N≧20，且主影像的时间空间基线在N+1景影像中最小。PS流程处理的最大优势是不用进行LOS方向的滤波。在选取主影像时，要建立相关函数模型：

上述公式(1),(2)中，

是SAR影像N和M的空间基线，T^n,m为影像的时间基线，

为多普勒基线，B_c，T_c，和f_c为临界基线，γ^m max的SAR影像为主影像，则由式(1)可知，当空间基线和时间基线值最小时，γ^m值最大，因此主影像的时间空间基线应是多景SAR影像中最小的。

对于地形相位和平地相位，PS-InSAR流程中是采取加入DEM来减少影响，在干涉相位有5部分组成，在DEM校正中上述两种相位将被去除，其中高程上由于DEM选用的精度会存在误差，因此对于干涉图N的像元，其相位为：

式(6)中，

为参考DEM校正中由于DEM精度引起的残余地形相位，

为大气相位，

为形变相位，

为其余噪声，由式(3)可知，DEM精度越高则残余地形相位会偏小，则有：

将(4)，(5)中的

带入(3)可知：

对于(6)中的大气相位和噪声等非线性相位，在覆盖大量混凝土的城市、道路或高速铁路等稳定地表的区域，二者的影像可视为线性形变。

经过上述步骤，会产生大量的永久散射点(PS)点，除了有形变信息的PS点外，还有许多对形变观测无用的点，所以需对PS点进行筛选。在本步骤中采取相干系数阈值法和振幅离差阈值法对PS点进行筛选，选取相干性高的PS点作为PS点，其中相干系数阈值法以一个局部矩阵进行像元的相干性计算，具体为：

式(7)中，γ是相干系数，M，S是两景SAR影像，m，n是局部矩阵的行列数，γ计算结果值为(0,1)，相干系数阈值是在进行运算之前给定阈值，大于阈值选为PS点，且相干系数值越接近于1表示其相干性越强。

相反的，振幅离差阈值是针对影像的单个像元，给定一个阈值范围，以像元振幅的稳定性来选取PS点，具体表现：

式(8)中，σ_φ和D_A表示相位标准差和振幅离差指数，μ_A为均值振幅，σ_A是振幅标准差，σ_nR为噪声实部标准差，σ_nl为噪声虚部标准差。

进行选点之后就可以进行解算，其中包括去除大气延迟、轨道误差等，然后进行线性形变估算，最终得到PS点的形变速率，对其进行地理编码即可得到形变结果。

不同于PS-InSAR流程处理，SBAS-InSAR流程处理无需对N景SAR影像选取主影像，在用于生成干涉图的SAR影像时选取时间空间基线较小的SAR影像形成干涉组对，以减少失相干对结果的影响，N景SAR影像生成M幅干涉图，干涉图的相位表示：

式(9)中，x,r是干涉图坐标，t_A,t_B为生成干涉图的两景影像的获取时间，若t_A<t_B，式(9)中第一项即为形变相位，式中的d(t_B,x,r),d(t_A,x,r)为t_A,t_B时刻(x,r)相对于t0时刻视线方向的累积形变量；第二项即为残余高程相位；第三项为大气延迟相位；余式为噪声相位；若假设M景主影像和从影像的获取时间以序列IE＝[IE₁...IE_M]与IS＝[IS₁...IS_M],若对M幅干涉图先不考虑大气延迟，只求解高程误差和形变量，则：

且将形变相位转化为平均相位速度矢量[v]^T：

若速度矢量是线性的，即v＝Mp，则有矩阵方程：

DMp+CΔz＝Δφ (12)

式(12)中D是M×(N-1)的矩阵，p为速度模型的参数矢量，则参数个数i应小于等于(N-1)，M为(N-1)×i的矩阵；Δz为高程相位误差；对于式(12)采用最小二乘法求解，可得到形变相位和高程误差相位，并可以去除相位图中高程相位误差，则矩阵方程可修改为：

Dv+CΔz＝Δφ (13)

在SBAS处理中，由于考虑到空间失相干的问题，将影像分成若干基线对，导致式(13)中的D为秩亏矩阵，若基线对为L(L>1)，则D的秩为N-L，使用SVD法对式(13)求解，对得到的速度矢量积分即可获得形变相位序列，需要注意的是此时结果中还有大气相位，一般采取高通滤波去除大气相位，得到研究区的形变序列,即研究区的地面形变信息。

步骤四：高铁线路背景区域沉降分离。

首先需对地面形变数据结果校准，高速铁路本专利建议采用PS-InSAR流程处理，所得结果可直接与二等水准控制点结果进行校正。校正后需要对地面形变区域进行分离，本专利在获取地面形变数据的基础上采用小波变换分析，针对形变区域进行高频低频滤波的小波分析，除了可以将部分噪声去除，还可以将沉降结果进行重构，从而反映出沉降区域的细节部分与趋势部分的差异，其中，通过高频滤波组的可很好的表达高铁线路及周围沉降区域细节性的信息；通过低频滤波的可很好的反映整区域趋势性相关信息。

步骤五：最大信息系数赋值权重因子。

对InSAR处理最终得到的影像进行小波变换处理可以将线状的高速铁路沿线与背景区进行分离。其在与小波变换分析可对栅格影像进行分离重构。

小波变换分析是在一维离散小波变换演变而来，小波变换最大的特点是在完成一次变换，影像都会分成高频图组与低频图组，其中高频组会更好的反映锐化的部分。

对于一维的信号：

定义

其逆变换为：

上式中Ψ为傅里叶变换，C_Ψ取值为：

对于栅格影像的二维的小波变换，有一个前提条件是尺度可分离的，即信号函数Φ(x,y)＝Φ(x)Φ(y),其中Φ(x)Φ(y)分别对应的小波函数为Ψ(x)Ψ(y)其创立二维小波变换的三个基本小波可表达为：

ψ¹(x,y)＝φ(x)ψ(y),ψ²(x,y)＝φ(y)ψ(x),ψ³(x,y)＝ψ(x)ψ(y) (16)

式中构建的二维平方可积函数空间L²(R²)的正交归一，既有：

l＝1,2,3,...且j,l,m,n均为整数。

对于一幅N×N的影像f(x,y)开始，其中上标尺度N，N＝2n，对于j＝0，尺度2j＝20＝1，为源影像尺度。j每增大一次使尺度加倍，分辨率减半，影像分为三个高频组，一个低频组，低频组中包括大部分能量，还需进行分离，高频组则更好的反映锐化等凸显线条的影像。对于第一个层次(j＝1)，则有：

通过小波变换分解后重构，将InSAR处理结果中有关高速铁路的部分可以提取出来，高频组中将高铁及沿线研究区的信息更加突出的表现，有利于针对高铁沿线展开分析。

在步骤一中，已对影响高铁行车安全的危害性因素诸如地下水位、土壤压缩模量、交通流量、压缩层厚度等数据，利用GIS软件工具，将数据进行统一的栅格化处理，对于数据点，首选克里金插值法进行插值，对于遥感地面沉降数据、光学遥感数据等栅格数据，可直接针对栅格灰度值进行栅格计算。将危害因子统一转为栅格数据后，将根据最大信息系数(MIC)算法，对栅格灰度值之间进行计算，从而完成危害因子赋值。

MIC算法是由互信息最大系数算法(MIM,Mutual Information Maximization)算法优化而来，通过对熵值的最大化，优化多特征值中互信息法存在的偏差，针对两种随机变量X,Y，I(X；Y)为其熵值，H(X),H(Y)为两种变量的熵值，I(X；Y)的上界为H(X),H(Y)之间的最小值，Imax(X；Y)为两种变量的最大信息系数

I(X；Y)＝H(X)-H(X|Y)＝H(Y)-H(Y|X) (18)

0≦I(X；Y)≦MIN{H(X),H(Y)} (19)

I_max(X；Y)取值范围为[0,1]，当MIC值趋近于0时表示两因子之间相互独立；趋近于1时表示两因子相关性较强。将沉降值作为变量X，分别将危害性因素作为变量Y分别与沉降值进行MIC计算，每个危害因子所得出的MIC值将作为危害性评估的权重进行赋值，从而减少对于主观经验的依赖，并且可根据地形因子条件设置不同的危害因子进行权重赋值。其中注意的是该步骤所进行计算除了数字计算还进行影像灰度值计算，输出结果包括进行灰度值计算后的栅格影像图，并在此基础上进行综合评估。

S6地面沉降危害性综合评估，随机森林算法在遥感影像的应用中已经取得了很好的成果。

1、本专利中将根据随机森林算法，基于地面沉降信息建立CRAT决策树，CRAT决策树的优势在于可监督的建立训练样本，但建立训练样本的工作量小，特征提取的精度高。基于上述步骤将赋权重的危害性评估因子与地面沉降图建立决策树，并且将处理后的光学遥感影像辅助进行分类，用以对高速铁路沿线不同地段的危害性进行评估。

作为本专利最重要的内容，对于地面沉降的危害性评估采用的是基于机器学习的随机森林算法，思路是通过上述步骤的提取、分区、赋权重等处理，将对所得结果建立决策树，进行特征提取，从而对地面沉降区域进行提取，做到对于危害区域的划分，达到危害性评估的目的。

在建立决策树前，需对影像建立训练样本，选取决策树算法的优点是可以在缺少先验性专家知识的情况下建立有效的训练样本，通过构建训练样本以确定特征属性和目标变量的关系。构建决策树训练样本的方法有：CART算法，C4.5和ID3，本专利中采用的是CART算法。CART是通过基尼系数来确定属性节点分割属性，其中选取属性基尼系数最小的进行二分化，最后构建一个二叉树的决策树。对于一个栅格数据集D，它的基尼系数为：

式中p_i是数据集D中任意数据被抽取成为训练样本的概率，对生成的二分数据集进行二分则有：

式中D₁，D₂为二分后的两个数据集，D₁为可进行二分的子数据集，并对D₁进行式22的计算再进行二分，并在节点处形成训练规则，直到所有二分子集都符合节点的规则时，结束算法，若有不符合的，可对其再次进行式22的计算，直到结束算法。

此步骤中利用CART构建决策树是为随机森林特征提取构建决策树。

2、其中，结合Boruta算法辅助修改各级决策树，Boruta通过对决策树建立副本并对比新建随机副本与原始变量的重要性对决策树进行修改，从而将对高铁沿线行车安全最重要的地段进行提取标记。由于进行计算的影像都可已进行栅格计算，因此此方法可对赋值的栅格影像进行计算，根据权重提取特征地段，实现了对高铁沿线沉降危害性的分区与评估。

随机森林利用组合的思想，对决策树进行改进，由单棵决策树组合成多棵决策树，形成森林。其中每棵决策树构建所需的数据集，是通过随机有放回的对原始数据集抽样获取的，这样保证了训练数据集的不相关性。决策树构建过程中分裂节点所需的属性，也是利用随机抽取的方式，对所有属性抽取部分属性，形成属性子空间，由此可知，每棵决策树选取的属性子空间都不会相同。

将步骤五中构建的带有危害因子、地面沉降信息的决策树与光学遥感影像，土地年鉴等资料再构建决策树，构建方法仍为CART。

在构建随机森林时本专利引入Boruta算法辅助进行危害因子评估。Boruta算法可在随机森林算法中对特征变量的重要性进行评价，本专利基于此思路对输入的因子要素进行评价。如上文中所述，Boruta算法是对原变量增加随机变量副本，通过对比原变量与副本变量的重要性进行重要性评分提取特征，其中Z值为重要性评分，其中Boruta算法在进行重要性评价时主要基于RF算法的袋外误差定义：

式23中，MSE_OOB为随机森林的袋外误差；y_i为样本值；

为样本y的袋外样本的预测值，基于此，则有：

式24中，Z_socre即为Z值；

为袋外误差平均值；SDMSE_OOB为袋外误差的标准差。

重要性评分Z值取最大值，记为shadowMax，当特征变量的Z值大于shadowMax时，认为该特征变量重要，参与后续建模；反之，则认为不重要，不参与后续计算。

利用Boruta对输入构建随机森林的因子进行重要性计算，在基于MIC权重赋值的情况下，随机森林会将权重值较高的因子划为重要性变量进行特征提取，因此在特征提取时，在步骤3中被分离的高频组高铁沿线影像和步骤4中赋予较高权重值的因子会被视为重要变量被提取，分类，从而使得地面沉降严重的区域得以与其他地区分离，达到划分高危害性地区和影像因子评估的目的；其他权重值相比沉降因子低的因子，会在此步骤中与参与进来的光学遥感影像构建决策树，对其进行提取分类，由于权重值较低，在Boruta对随机森林辅助计算得到地类的分类图，结合危害因子，完成对非沉降严重区域的评估，最终完成研究区域内地面沉降的危害性评估。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (10)

1.一种基于多源数据集成的高铁沉降危害性评估方法，其特征在于，包括以下步骤：

S1、利用InSAR技术采集地面形变数据，包括采集InSAR影像数据、光学遥感影像数据、水文地质资料数据、二级水准测量数据；

S2、对步骤S1中采集的数据进行预处理；

S3、采用PS-InSAR+SBAS-InSAR两种时序InSAR技术分析提取地面形变结果；

S4、针对高铁线路背景区域沉降分离；

S5、对影响地面形变的危害因子采用MIC算法赋权重；

S6、地面沉降危害性评估。

2.根据权利要求1所述的基于多源数据集成的高铁沉降危害性评估方法，其特征在于：所述步骤S2中，对采集的数据进行预处理的方法包括：

对光学遥感影像数据的辐射定标、大气校正、几何校正；

对水文地质资料数据的插值化、矢量化栅格化处理；

对InSAR影像数据的去噪、多视、选取主影像。

3.根据权利要求1所述的基于多源数据集成的高铁沉降危害性评估方法，其特征在于，所述步骤S3中，针对有大量为稳定散射体的高铁沿线的处理方法为：进行PS-InSAR流程处理，在多时序SAR影像中选取时间空间基线小的影像为主影像，通过相干系数阈值和振幅离差阈值两种方法对影像进行选点处理，选取稳定的PS点；

针对有大量荒地或逐步建成的区域的处理方法为：进行SBAS-InSAR流程处理，除去必要的去噪之外，不必进行主影像选取，在影像基线优化处理时，引入最小二乘法优化干涉对组，在保障影像处理质量的同时，减少冗余干涉对，减轻计算负担。

4.根据权利要求3所述的基于多源数据集成的高铁沉降危害性评估方法，其特征在于，进行PS-InSAR流程处理的具体方法如下：

将InSAR影像数据进行PS-InSAR处理，需要在N+1景影像中选取1景主影像m，其余N景InSAR影像将与主影像匹配，生成N幅干涉图，基于PS处理流程的特殊性，一般要求N≧20，且主影像的时间空间基线在N+1景影像中最小，在选取主影像时，要建立相关函数模型：

上述公式(1)、(2)中，

是SAR影像N和M的空间基线，T^n,m为影像的时间基线，

为多普勒基线，B_c，T_c，和f_c为临界基线，γ^mmax的SAR影像为主影像，则由式(1)可知，当空间基线和时间基线值最小时，γ^m值最大，因此主影像的时间空间基线应是多景SAR影像中最小的；

对于地形相位和平地相位，PS-InSAR流程中是采取加入DEM来减少影响，在干涉相位有5部分组成，在DEM校正中上述两种相位将被去除，其中高程上由于DEM选用的精度会存在误差，因此对于干涉图N的像元，其相位为：

式(6)中，

为参考DEM校正中由于DEM精度引起的残余地形相位，

为大气相位，

为形变相位，

为其余噪声，由式(3)可知，DEM精度越高则残余地形相位会偏小，则有：

将(4)，(5)中的

带入(3)可知：

对于(6)中的大气相位和噪声等非线性相位，在覆盖大量混凝土的城市、道路或高速铁路稳定地表的区域，二者的影像可视为线性形变；

经过上述步骤，会产生永久散射点即PS点，除了有形变信息的PS点外，还有许多对形变观测无用的点，所以需对PS点进行筛选；在本步骤中采取相干系数阈值法和振幅离差阈值法对PS点进行筛选，选取相干性高的PS点作为PS点，其中相干系数阈值法以一个局部矩阵进行像元的相干性计算，具体为：

式(7)中，γ是相干系数，M，S是两景SAR影像，m，n是局部矩阵的行列数，γ计算结果值为(0,1)，相干系数阈值是在进行运算之前给定阈值，大于阈值选为PS点，且相干系数值越接近于1表示其相干性越强；

相反的，振幅离差阈值是针对影像的单个像元，给定一个阈值范围，以像元振幅的稳定性来选取PS点，具体表现：

式(8)中，σ_φ和D_A表示相位标准差和振幅离差指数，μ_A为均值振幅，σ_A是振幅标准差，σ_nR为噪声实部标准差，σ_nl为噪声虚部标准差；

进行选点之后就可以进行解算，其中包括去除大气延迟、轨道误差等，然后进行线性形变估算，最终得到PS点的形变速率，对其进行地理编码即可得到形变结果。

5.根据权利要求3所述的基于多源数据集成的高铁沉降危害性评估方法，其特征在于，进行SBAS-InSAR流程处理的具体方法如下：

不同于PS-InSAR流程处理，SBAS-InSAR流程处理无需对N景SAR影像选取主影像，在用于生成干涉图的SAR影像时选取时间空间基线较小的SAR影像形成干涉组对，以减少失相干对结果的影响，N景SAR影像生成M幅干涉图，干涉图的相位表示：

式(9)中，x,r是干涉图坐标，t_A,t_B为生成干涉图的两景影像的获取时间，若t_A<t_B，式(9)中第一项即为形变相位，式中的d(t_B,x,r),d(t_A,x,r)为t_A,t_B时刻(x,r)相对于t0时刻视线方向的累积形变量；第二项即为残余高程相位；第三项为大气延迟相位；余式为噪声相位；若假设M景主影像和从影像的获取时间以序列IE＝[IE₁...IE_M]与IS＝[IS₁...IS_M],若对M幅干涉图先不考虑大气延迟，只求解高程误差和形变量，则：

且将形变相位转化为平均相位速度矢量[v]^T：

若速度矢量是线性变换的，即v＝Mp，则有矩阵方程：

DMp+CΔz＝Δφ (12)

式(12)中D是M×(N-1)的矩阵，p为速度模型的参数矢量，则参数个数i应小于等于(N-1)，M为(N-1)×i的矩阵；Δz为高程相位误差；对于式(12)采用最小二乘法求解，可得到形变相位和高程误差相位，并可以去除相位图中高程相位误差，则矩阵方程可修改为：

Dv+CΔz＝Δφ (13)

在SBAS处理中，由于考虑到空间失相干的问题，将影像分成若干基线对，导致式(13)中的D为秩亏矩阵，若基线对为L(L>1)，则D的秩为N-L，使用SVD法对式(13)求解，对得到的速度矢量积分即可获得形变相位序列，需要注意的是此时结果中还有大气相位，一般采取高通滤波去除大气相位，得到研究区的形变序列,即研究区的地面形变信息。

6.根据权利要求1所述的基于多源数据集成的高铁沉降危害性评估方法，其特征在于，所述步骤S4的具体方法如下：

S401、对步骤S3中采集的地面形变数据结果进行校准，针对高速铁路采用PS-InSAR流程处理，将结果与二等水准测量数据进行校正；

S402、对形变区域进行高频低频滤波的小波变换分析，用于反映出沉降区域的细节部分与趋势部分的差异；

其中，通过高频滤波组的用于表达高铁线路及周围沉降区域细节性的信息；通过低频滤波的用于反映整区域趋势性相关信息。

7.根据权利要求6所述的基于多源数据集成的高铁沉降危害性评估方法，其特征在于，所述步骤S402中，小波变换分析的具体方法如下：

利用小波变换完成一次变换，将影像数据分成高频图组与低频图组，其中高频组会更好的反映锐化的部分；

对于一维的信号：

定义

其逆变换为：

上式中Ψ为傅里叶变换，C_Ψ取值为：

对于栅格影像的二维的小波变换，有一个前提条件是尺度可分离的，即信号函数Φ(x,y)＝Φ(x)Φ(y),其中Φ(x)Φ(y)分别对应的小波函数为Ψ(x)Ψ(y)其创立二维小波变换的三个基本小波可表达为：

ψ¹(x,y)＝φ(x)ψ(y),ψ²(x,y)＝φ(y)ψ(x),ψ³(x,y)＝ψ(x)ψ(y) (16)

式中构建的二维平方可积函数空间L²(R²)的正交归一，既有：

l＝1,2,3,...且j,l,m,n均为整数；

对于一幅N×N的二维栅格影像f(x,y)开始，其中上标尺度N，N＝2n，对于j＝0，尺度2^j＝2⁰＝1，为源影像尺度；j每增大一次使尺度加倍，分辨率减半，影像分为三个高频组，一个低频组，低频组中包括大部分能量，还需进行分离，高频组则更好的反映锐化等凸显线条的影像；对于第一个层次(j＝1)，则有：

通过小波变换分解后重构，将InSAR影响处理结果中有关高速铁路的部分提取出来，高频组中将高铁及沿线研究区的信息更加突出的表现。

8.根据权利要求1所述的基于多源数据集成的高铁沉降危害性评估方法，其特征在于，所述步骤S5的具体方法如下：

S501、针对步骤S1中采集的数据，选择影响高铁行车地面形变的危害因子，包括地下水位、土壤压缩模量、交通流量、压缩层厚度；

S502、利用GIS软件工具，将数据进行统一的栅格化处理，对于数据点，首选克里金插值法进行插值，对于遥感地面沉降数据、光学遥感数据栅格数据，直接针对栅格灰度值进行栅格计算；

S503、将危害因子统一转为栅格数据后，根据MIC算法，对栅格灰度值之间进行计算，从而完成危害因子赋权重。

9.根据权利要求1所述的基于多源数据集成的高铁沉降危害性评估方法，其特征在于，所述步骤S6的具体方法如下：

S601、基于步骤S5中得到的赋权重的危害因子与地面沉降图并结合随机森林算法建立决策树，并且将处理后的光学遥感影像数据进行分类，用以对高速铁路沿线不同地段的危害性进行评估；

S602、结合Boruta算法辅助修改各级决策树，Boruta算法通过对决策树建立副本并对比新建随机副本与原始变量的重要性对决策树进行修改，用于对高铁沿线行车安全地段进行提取标记，此方法可对赋值的栅格影像进行计算，根据权重提取特征地段，实现了对高铁沿线沉降危害性的分区与评估。

10.根据权利要求9所述的基于多源数据集成的高铁沉降危害性评估方法，其特征在于：在执行步骤S601前，在建立决策树前，需对影像建立训练样本，通过构建训练样本以确定特征属性和目标变量的关系；采用CART算法构建决策树训练样本，CART算法通过基尼系数确定属性节点分割属性，其中选取属性基尼系数最小的进行二分化，最后构建一个二叉树的决策树：对于一个栅格数据集D，它的基尼系数为：

式中p_i是数据集D中任意数据被抽取成为训练样本的概率，对生成的二分数据集进行二分则有：

式中D₁，D₂为二分后的两个数据集，D₁为可进行二分的子数据集，并对D₁进行式22的计算再进行二分，并在节点处形成训练规则，直到所有二分子集都符合节点的规则时，结束算法，若有不符合的，可对其再次进行式22的计算，直到结束算法。

Images