[0001] 本发明属于地质灾害监测技术领域，具体为一种基于多因子随机森林回归模型的地表连续形变图的重建方法，可应用于大范围地质灾害连续形变场回复，用于地质灾害调查、监测与防治。

[0002] 传统的地表形变监测方法，如全球导航卫星系统(Global Navigation Satellite System，GNSS)和水准测量，受限于其空间覆盖率和分辨率较低，不适合大规模连续形变监测。合成孔径雷达干涉测量(Interferometric Synthetic Aperture Radar，InSAR)在高空间分辨率和大规模形变监测方面具有显著优势。目前，InSAR技术在城市地表沉降监测中的应用已十分广泛。然而，InSAR技术在进行连续地表形变监测时面临两大技术难题，一是较长的时空基线、植被以及地表扰动导致去相干，二是较大的形变梯度导致相位解缠错误，导致InSAR地表形变图中存在监测数据缺失，极大地限制了InSAR变形图的时空特征分析和机理反演研究。

[0003] 为解决该问题，现有技术通常采用空间插值方法或机器学习方法来预测和重建地表形变图。首先常规插值方法，如反距离(IDW)插值法、克里金(Kriging)插值法等是基于已知和未知数据具有相似的统计或几何结构进行预测的，预测结果仅具有数学意义。但是，插值方法的预测结果很容易受到已知测量数据空间分布的影响。在测量值均匀分布时效果最佳，如果观测数据稀疏且分布不均，内插结果的精度就会降低。而常规的机器学习地表变形预测模型一般忽略由地下水位变化和地裂缝分布引起的地表沉降的空间异质性，导致地表变形预测精度较低。

[0004] 然而，InSAR大范围地表形变测量中，相干点的分布并不均匀，而且不同区域的地表形变影响因素和机制也不相同，导致现有方法难以获得高精度的地表形变重建结果。

[0052] 以下描述中，为了说明而不是为了限定，提出了诸如特定系统结构、技术之类的具体细节，以便透彻理解本发明实施例。然而，本领域的技术人员应当清楚，在没有这些具体细节的其它实施例中也可以实现本发明。在其它情况中，省略对众所周知的系统、装置、电路以及方法的详细说明，以免不必要的细节妨碍本发明的描述。

[0053] 本发明实施例的第一方面提供了一种地表连续形变图的重建方法，如图1所示，包括：

[0054] 步骤1、根据研究区域的地表形变图，确定待重建区域，具体包括：

[0055] 步骤1.1、将研究区域的地表形变图划分为多个区域。

[0056] 本发明实施例中，首先获取研究区域的SAR影像数据集，然后利用差分干涉测量短基线集时序分析技术(SmallBaselineSubsetInSAR，简称SBAS‑InSAR)获取地表沉降速率值，得到地表形变图。

[0057] 本发明实施例中，将研究区域的地表形变图划分为多个区域，具体可以利用网片对地表形变图进行划分，划分为连续的多个网格，网格的大小主要取决于后续使用的影响因子数据集的分辨率。本发明实施例中，使用网格单元内所有地表沉降速率值的平均值作为网格测量值。

[0058] 示例性地，将地表形变图通过Fishnet划分为100m*100m的连续的多个网格。

[0059] 步骤1.2、将不包含地表形变数据的区域确定为待重建区域。

[0060] 示例性地，不包含地表形变数据的区域为不包含地表形变数据的网格，将没有地表沉降速率值的网格，确定为待重建区域。

[0061] 在步骤1之后，还包括：

[0062] 步骤S1：获取历史多源异构数据集。其中历史多源异构数据集中包括研究区域的多源异构的影响因子对应的历史数据。

[0063] 示例性地，多源异构的影响因子包括地下水、降雨、地裂缝、地层岩性、地貌、水文地质、工程地质、用地类型、土壤类型、GDP、DEM等。

[0064] 为保证后续模型训练的精度，本发明实施例将上述多源异构影响因子的历史数据统一投影到WGS84坐标系中。然后将如地下水、地裂缝的矢量因子历史数据转换为栅格数据，并将所有栅格数据集重采样至分辨率为100*100m。进一步地，对定类因子进行分级量化处理，再与栅格数据同时进行归一化处理。

[0065] 步骤S2：从历史多源异构数据集中，确定与地表形变数据关联度大于预设第一阈值的多个不同源的历史异构数据，具体包括：

[0066] 步骤S21：获取历史多源异构数据集中每种历史异构数据在地表形变图中的频率比，具体包括：

[0067] 获取历史多源异构数据集中每种历史异构数据影响的沉降区域数与沉降区域总数之间的第一比值；

[0068] 获取历史多源异构数据集中每种历史异构数据在地表形变图中所占区域与区域总数之间的第二比值；

[0069] 根据第一比值和第二比值，确定每种历史异构数据在地表形变图中的频率比。

[0070] 示例性地，当使用网格划分研究区域的地表形变图时，步骤S21中的区域即为网格，例如：沉降区域数为沉降网格数。

[0071] 本发明实施例中，忽略沉降率值小于5mm/年的区域，叠加多源异构因子和沉降区域，计算沉降位置的数量和与每个类别关联的像素。

[0072] 示例性地，利用如下公式(1)确定每种历史异构数据在地表形变图中的频率比：

[0073]

[0074] 式中，FR表示每种历史异构数据在研究区域的地表形变图中的频率比，A是每种历史异构数据影响的地表沉降网格数；A'是地表沉降网格的总数；B是每种历史异构数据在地表形变图中所占网格数；B'是地表形变图中网格的总数。

[0075] 步骤S22：根据频率比和地表沉降速率值的关联度值，确定与地表沉降速率值关联度大于预设第一阈值的多个不同源的历史异构数据。

[0076] 本发明实施例中，计算关联度值的方法可使用现有技术中的常规方法，例如可利用灰色关联分析法计算频率比和地表沉降速率值的关联度值，筛选出与地表沉降关联度大的因子，剔除关联度小的因子。

[0077] 步骤S3、对大于预设第一阈值的多个不同源的异构数据进行聚类，得到聚类多源异构数据。

[0078] 本发明对大于预设第一阈值的多个不同源的异构数据进行聚类，所采用的聚类方法可为K‑Means聚类方法、层次聚类方法、SOM聚类方法或者FCM聚类方法等。本发明实施例优选使用K‑Means聚类方法。

[0079] 采用K‑Means聚类方法对不同源的异构数据进行聚类的过程为：样本以迭代方式分组为K类，聚类质心逐个更新，以最小化每个样本与其所属类的中心或均值之间的距离之和d，直到获得最佳聚类效果，其中d根据公式(2)确定：

[0080]

[0081] 式中，xij为第j类的第i个影响因素，cj为第j类的中心，k为聚类簇个数。

[0082] 待聚类趋于稳定后，根据肘部法则公式确定最优分类，将研究区域划分为若干个同质区域。

[0083] 步骤2、基于研究区域的历史地表形变数据和历史聚类多源异构数据训练地表沉降预测模型。

[0084] 本发明实施例中的沉降预测模型为随机森林回归模型。

[0085] 本发明实施例在训练过程中，将聚类后的历史聚类多源异构数据和研究区域的历史地表形变数据划分为70％的训练集和30％的测试集，迭代测试确定最优参数构建地表沉降预测模型。其中研究区域的历史地表形变数据为历史地表沉降速率值。

[0086] 步骤3、将待重建区域的当前聚类多源异构数据输入至训练后的沉降预测模型中，得到待重建区域的预测形变图，具体包括：

[0087] 步骤31、将待重建区域的当前聚类多源异构数据输入至训练后的随机森林回归模型中，得到多个回归决策树的预测形变图；

[0088] 步骤32、基于综合学习的思想，计算多个回归决策树的预测形变值的均值，得到待重建区域的预测形变图，如公式(3)：

[0089]

[0090] 式中： 为待重建区域的预测形变图；h(x，θt)为基于x和θt的输出，x为自变量，θt为独立同分布的随机向量；T为回归决策树的个数。

[0091] 本发明实施例利用R2、平均绝对误差(MAE)、均方根误差(RMSE)评估模型的精度。2
R、平均绝对误差(MAE)、均方根误差(RMSE)的常规公式分别为公式(5)至公式(7)。

[0092]

[0093]

[0094]

[0095] 式中，yi为形变真实值， 为形变值均值， 为形变预测值。

[0096] 步骤4、将地表形变图和预测形变图进行融合，得到大范围研究区域内的高精度的地表连续形变图。

[0097] 进一步地，为避免共线性因子数据影响建模精度，本发明实施例在步骤S3之前还包括：

[0098] 获取与地表沉降速率值关联度大于预设第一阈值的多个不同源的历史异构数据中，任意两种历史异构数据之间的相关性；

[0099] 当相关性大于预设第二阈值时，从对应的两种历史异构数据中选取其中一种，得到相关性处理后的多个不同源的历史异构数据；

[0100] 相应的，步骤S3具体为：

[0101] 对相关性处理后的多个不同源的历史异构数据进行聚类，得到聚类多源异构数据。

[0102] 本发明计算任意两种历史异构数据之间的相关性可采用皮尔逊相关系数法或者斯皮尔曼等级相关系数法，本发明实施例优选为斯皮尔曼等级相关系数法。使用Spearman相关系数法来描述两向量的相关性，此处的两变量为两种历史异构数据。计算2个n维向量x，y的Spearman相关系数方法如式：

[0103]

[0104] 式中Ri和Si各代表向量x，y进行排序后观测值i的次序， 和 则表示向量x和y的平均次序，N是i的总数量，di＝Ri‑Si表示两个变量中i的次序之差。Spearman相关系数中ρs是[‑1,1]中的实数。当相关系数大于零时，两变量成正相关，反之，则负相关。

[0105] 本发明实施例中，取第二阈值为0.8。

[0106] 本发明提出了一种多因子随机森林回归模型重建地表连续形变图的方法，该方法综合考虑影响地表沉降的多种影响因素(多源异构数据)及其各因素在空间的异质性，高精度重建大范围地表连续形变场。

[0107] 本发明的第二方面提供了一种地表连续形变图的重建系统，包括区域确定模块、模型训练模块、预测模块和重建模块。

[0108] 其中，区域确定模块用于根据研究区域的地表形变图，确定待重建区域；

[0109] 模型训练模块用于基于研究区域的历史地表形变数据和历史聚类多源异构数据训练地表沉降预测模型；

[0110] 预测模块用于将待重建区域的当前聚类多源异构数据输入至训练后的沉降预测模型中，得到待重建区域的预测形变图；

[0111] 重建模块用于将地表形变图和预测形变图进行融合，得到研究区域的地表连续形变图。

[0112] 本发明考虑多因子异质性的随机森林回归模型方法重建(即在空间域预测)大范围连续地表形变信息，能够充分利用不同源的异构环境因子数据，减少对于已有历史数据的依赖性；同时，利用机器学习建模预测增加了预测结果的可解释性，预测结果更加可靠、精度更高；此外，采用K‑Means聚类方法对多源异构数据进行聚类，可以大大降低不同源因子空间异质性对于地表沉降预测模型构建的影响，提高建模精度。该方法适用于大范围空间形变预测，该方法得出的连续地表形变结果有利于时空地表变形特征描述、机理解释和预防决策。

[0113] 下面将以更为具体的实施例详述本发明的方法。

[0114] 本实施例选择位于我国秦岭北部的西安市。西安地貌复杂，构造作用下形成的黄土洼地多达10个，活动性裂缝数个，其中以渭河断裂和临潼—长安断裂最为典型。由于人口快速增长和工业用水大幅增加，地表水难以满足需求，地下水成为西安市重要的供水来源。地下水的长期过度开采是地表沉降的原因之一。随着西安经济的快速发展和城市建设的扩大，地下工程和高层建筑的频繁建设成为地表沉降的另一个原因。地质裂缝、断裂带等自然因素与地下水开采等人为因素共同作用，导致西安部分地区呈现区域性地表沉降。

[0115] 在本实施例中，为获取研究区域的形变特征，我们以2012年1月‑2015年5月的40景降轨的德国TerraSAR‑X数据为基准数据集，利用SBAS‑InSAR技术获得了研究区地表沉降速率值。同时获取了承压水水位、承压水变化、潜水水位、潜水变化、降雨量五种动态影响因素，以及地裂缝、地层岩性、地貌、水文地质、工程地质、用地类型、土壤类型、GDP、DEM共9种静态影响因素。

[0116] 采用SBAS‑InSAR技术计算了西安市2012、2013、2014和2015年的年平均变形率图，如图2所示。然而，测量结果中存在许多测量信息缺失区域，这使得分析时空特征变得困难。因此，利用一种多因子随机森林回归模型重建地表连续形变图的方法(又称K‑RFR方法)进行重建西安市连续地表形变图，流程如图3。将平均变形率图离散化为100m*100m网格单元，测量区域和不含测量信息的区域分别含有47810和42496个单元，即大约47％的面积需要重建才能获得连续变形图。网格离散化示意图如图4。

[0117] 图5为利用本发明所选取方法重建后的连续地表形变图。重建区域可以达到100m的空间分辨率，预测结果与测量结果具有良好的空间连续性。为了观察重建形变结果的细节，图6展示了三个典型区域。图6(a)显示了InSAR技术获得的形变速率，图6(b)显示了使用本发明所选方法在不含测量信息区域预测的2012‑2015年的形变速率，图6(c)显示了测量结果与预测结果融合后的结果，图6(d‑m)显示了三个典型区域的形变细节图。

[0118] 为了进一步证明K‑Means聚类分析提高了地基变形预测模型的精度，分别计算了聚类前后随机森林回归模型构建的模型预测结果与InSAR测量值之间的绝对误差和绝对误差分布，如图7所示。两种聚类方法都能很好地重建地表变形。大多数预测结果的绝对误差值小于6mm，误差小于沉降率的10％，能够满足预测的要求。与聚类前构建的模型的预测结果相比，本发明所选取的方法预测结果的标准差从2.352mm降低到1.782mm，验证了本发明所选方法具有更高的预测精度。

[0119] 图8为在I‑Ⅲ区域去除已有沉降点后，本发明所选方法重建的形变结果与IDW和普通Kriging重建的累积形变结果对比。其中，第一列为IDW方法，第二列为普通Kriging方法，第三列为本发明所选的方法。与IDW和普通Kriging等传统插值方法相比，本发明所选方法能很好地反映沉降的空间变化。图9为本发明所选方法和Kriging、IDW在I‑III区域重建结2 2
果精度对比图。以去除的高相干点作为验证数据，与IDW(R ＝0.63)和普通Kriging(R ＝
2
0.69)相比，本发明所选方法预测(R＝0.94)结果与InSAR监测结果的一致性高于两种空间插值方法。

[0120] 以上所述，仅是本申请的几个实施例，并非对本申请做任何形式的限制，虽然本申请以较佳实施例揭示如上，然而并非用以限制本申请，任何熟悉本专业的技术人员，在不脱离本申请技术方案的范围内，利用上述揭示的技术内容做出些许的变动或修饰均等同于等效实施案例，均属于技术方案范围内。