CN115982968B - 一种揭示地面时空间演化差异的算法 - Google Patents

Abstract

本发明公开了一种揭示地面沉降时空间演化差异的算法，属于地质灾害预防技术领域，包括以下步骤：基于雷达影像数据，利用SBAS‑InSAR技术，获取地面沉降信息；在构建地面沉降生存分析模型时，首先结合沉降信息和生存分析原理，确定沉降生存指标；然后以不同分类依据确定多种空间分区，带入生存分析算法，得到地面沉降生存模型，量化揭示地面沉降时序演化过程，并同时实现沉降空间区域差异对比。本发明揭示并量化了地面沉降的时序发展过程和空间演化差异，进一步掌握研究区地面沉降的时空演化差异，为城市灾害防治节约时间和经济成本，也可为城市规划决策提供科学依据。

Description

一种揭示地面时空间演化差异的算法

技术领域

本发明属于地质灾害预防技术领域，具体涉及种揭示地面沉降时空间演化差异的算法。

背景技术

地面沉降是人类工程活动(如地下水开采)或自然因素引起的一种地质现象，在这种现象下，一定区域内的地面高程降低。地面沉降的发生几乎是不可逆转的。当发生严重的局部地表沉降时，可诱发一系列地质环境灾害，如基础下沉、房屋开裂、地下管线破坏等；沉降具有形成时间长、影响范围广、防治难度大等特点。目前，地面沉降问题主要关注其空间分布和形成因素，而在时间序列背景下进行地面沉降区域化差异的研究还很少。

发明内容

针对现有技术存在的技术问题，本发明的目的在于提供了一种揭示地面沉降时空间演化差异的算法。

为实现上述目的，本发明采用了如下技术方案：

本发明提供了一种揭示地面沉降时空间演化差异的算法，包括以下步骤：

S1、获取覆盖研究区范围的长时序合成孔径雷达影像数据，并利用时序合成孔径雷达干涉测量方法，获取地表视线向的形变信息，然后通过实测GPS和水准点数据对InSAR监测结果进行校准和精度验证，将视线向的值转化为垂向，得到垂向形变数据，进而获取地面沉降速率；

S2、获取沉降生存指标：根据回归合成孔径雷达回归周期和研究区地面沉降发育程度，确定生存分析的时间单位，进而确定生存分析中的关注事件，即沉降生存指标；其中，所述关注事件以累计沉降值或沉降速率为阈值定义；

S3、根据区域沉降发育情况，确定起始事件I、终点事件E，计算得到每个沉降监测点的生存时间t，统计样本结局状态；

S4、根据研究区实际情况，通过Kaplan-Meier方法绘制生存函数；

S5、对研究区结合Log-Rank检验和等比例风险模型(Cox Proportional HazardsModel)分析，确立分区方案以及地面沉降影响因子，构建区域地面沉降生存分析模型，揭示地面沉降时空间演化差异，量化展示地面沉降时序演化过程和空间分布差异。

优选地，步骤S4中，所述根据研究区实际情况，通过Kaplan-Meier方法绘制生存函数，具体包括以下步骤：

S41、对于具有n个样本的沉降数据集，根据合成孔径雷达回归周期T和研究区地面沉降发育程度，确定生存分析的时间单位t：首先将卫星回归周期T确定为分析时间单位t，计算T_n-T_n-1内地面形变速率，若在该T_n-T_n-1内地面形变速率内无法观察到毫米级形变，则将分析时间单位扩大至2T₀，以此类推；

S42、在步骤A1确定生存分析时间单位t的基础下，分别计算：t₁、t₂、t₃…t_n时刻的累积地面沉降值和当n＝1、2、3…n时的t_n-t_n-1内的沉降速率；

S43、根据研究区在研究时段的地面沉降发育情况，确定沉降速率υ和累积沉降值s阈值，视为生存分析中的死亡事件，t_υ和t_s视为死亡时刻，并分别统计每个样本的生存时间t_υ-t₀和t_s-t₀；

S44、根据可能的影响因素对沉降样本点进行分组，并结合样本的生存时间t_υ-t₀绘制每组的Kaplan-Meier生存曲线S_(t)，绘制方法如下：

其中，d为t_n时刻发生死亡事件的个体数，m为t_n时刻仍存活的样本数。

优选地，步骤S5中，所述Log-Rank检验，用于比较步骤S4中得到的两组生存曲线之间有无差别，具体步骤如下：

假设步骤S4中得到的两组生存曲线一样即零假设，基于此假设，得到两组内的事件发生个体数之比等于两组样本数之比，由此计算出事件发生的期望数，分别将两组所有时间点的期望数相加，与所有观察数进行比较；如果零假设成立，那么总的事件发生数d应该按样本量等比例分布在两组之间，那么每组期望数(即期望的事件发生数)可以通过下式计算：

计算Log-Rank统计量：

根据得到的Log-Rank统计量数值判断是否拒绝零假设，进而判断两组生存曲线之间有统计学差异。

优选地，步骤S5中，所述对研究区进行等比例风险模型(Cox ProportionalHazards Model)分析，确立分区方案以及地面沉降影响因子，构建区域地面沉降生存分析模型，具体包括以下步骤：

S51、确定Cox回归分析要求比例满足风险假定的前提条件，即风险比不随时间的变化而变化：

HR＝EXP{β}

其中，HR代表即为假定风险比，β为自变量的偏回归系数；

S52、以λ(t，x)为因变量，等比例风险模型的基本形式为：

λ(t，x)＝λ₀(t)exp(β₁x₁+β₂x₂+...+β_mx_m)

其中，β为自变量的偏回归系数；λ(0)代表x为零向量时，λ(t，x)的基准危险率；

S53、因素i和j影响下的生存风险度相比，求得相对危险度RR：

与现有技术相比，本发明具有的有益效果如下：

(1)本发明首先基于雷达影像数据，利用SBAS-InSAR技术，获取区域地面沉降信息，然后根据沉降信息，确定研究区分区原则，带入生存分析算法，构建得到地面沉降生存模型，量化揭示了地面沉降时序演化过程，并同时实现了沉降空间区域差异对比。

(2)本发明将生存分析的方法应用于地面沉降研究中，首次提出将生存分析方法运用于地面沉降长时序多区域对比研究，将地面沉降时序发展演化与空间分布在不同影响下的差异进行量化揭示，为城市灾害防治节约时间和经济成本，也可为城市规划提供科学依据。

附图说明

图1为本发明提出的一种揭示地面沉降时空间演化差异的算法的流程示意图；

图2为京津冀地区的垂向形变数据结果图；

图3为京津冀地区地下水储量不同的区域间地面沉降生存情况结果图。

具体实施方式

下面将结合本发明实施例及附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。

实施例1

下面结合图1-图3，以京津冀地区作为研究区域，描述一种揭示地面沉降时序演变区域差异的算法，包括以下步骤：

S1、获取京津冀地区长时序合成孔径雷达影像数据，利用时序合成孔径雷达干涉测量方法，获取地表视线向的形变信息，通过实测GPS和水准点数据对InSAR监测结果进行校准和精度验证，将视线向的值转化为垂向，得到垂向形变数据，进而获取地面沉降速率，如图2所示；

S2、根据回归合成孔径雷达回归周期T₀和研究区地面沉降发育程度，确定生存分析的时间单位t，进而确定生存分析中的关注事件即沉降生存指标，本发明中的关注事件是以累计沉降值或沉降速率为阈值定义；

S3、根据研究区域沉降发育情况，确定起始事件I、终点事件E，计算得到每个沉降监测点的生存时间t，并统计样本结局状态；以累积沉降值为例，京津冀地区起始事件为雷达影像开始监测，终点事件为累计沉降值监测达到400mm；

S4、根据研究区域实际情况，通过Kaplan-Meier方法绘制生存函数，在生存分析中，持续时间T被假设为是随机的，则得到累积分布函数为：

F(t)＝Pr(T≤t) (1)

上述式(1)的含义为在时间t或t以内的失败(死亡)的概率，由此可得关注事件发生概率对应的密度函数为：

则生存函数S(t)定义为：

S(t)＝1-F(t) (3)

上述式(3)的含义为生存时间超过某个时间t的概率，即关注事件在时间t或t以内都未发生的概率。

生存函数S(t)为常用的生存分析方法，将其应用于地面沉降研究，能够在整体尺度上量化时间序列下地面沉降发育过程，并结合多条不同区域生存曲线探究地面沉降空间分布差异。图3为京津冀地区地下水储量不同的区域间地面沉降生存情况结果图；

其中，通过Kaplan-Meier方法绘制生存曲线，具体包括以下步骤：

S41、对于具有n个样本的沉降数据集，根据合成孔径雷达回归周期T和研究区地面沉降发育程度，确定生存分析的时间单位t：首先将卫星回归周期T确定为分析时间单位t，计算T_n-T_n-1内地面形变速率，若在该T_n-T_n-1内地面形变速率内无法观察到毫米级形变，则将分析时间单位扩大至2T₀，以此类推；

S42、在步骤A 41确定生存分析时间单位t的基础下，分别计算：t₁、t₂、t₃…t_n时刻的累积地面沉降值和当n＝1，2，3…n时的t_n-t_n-1内的沉降速率；

S43、根据研究区域在研究时段的地面沉降发育情况，确定沉降速率υ和累积沉降值s阈值，视为生存分析中的死亡事件，t_υ和t_s视为死亡时刻，并分别统计每个样本的生存时间t_υ-t₀和t_s-t₀；

S44、根据可能的影响因素对沉降样本点进行分组(此处的影响因素必须为分类变量)并结合样本的生存时间t_υ-t₀绘制每组的Kaplan-Meier生存曲线S_(t)，绘制方法如下：

其中，d为t_n时刻发生死亡事件的个体数，m为t_n时刻仍存活的样本数；

S5、对研究区结合Log-Rank检验以及等比例风险模型(Cox ProportionalHazards Model)回归分析方法，确立分区方案以及地面沉降影响因子，建立区域地面沉降生存分析模型，揭示地面沉降时空间演化差异，量化展示地面沉降时序演化过程和空间分布差异；同时对不同因素影响下的地面沉降差异进行同步展示。上述分析方法均可通过SPSS软件、Python、R语言等方式实现。

其中，上述Log-Rank检验是用于比较两组生存曲线之间有无差别，具体原理如下：

将Log-Rank检验进行零假设，即假设组A和组B的生存曲线一样，基于此假设，组A和组B内的事件发生个体数之比等于该两组样本数之比，由此计算出事件发生的期望数，分别将组A和组B两组所有时间点的期望数相加，与所有观察数进行比较(见表1)；

表1

	事件发生个体数	生存个体数	风险个体数
				组A	dA	nA-dA	nA
组B	dB	nB-dB	nB
				合计	d	n-d	n

如果零假设成立，那么总的事件发生个体数d应该按样本量等比例分布在组A和组B两组之间，那么每组期望数，即期望的事件发生数E_At通过下式(4)计算：

进而计算Log-Rank统计量为：

根据计算的Log-Rank统计量结果判断是否拒绝零假设，进而判断两组生存曲线之间有统计学差异；本实施例1中计算图2所示生存曲线的Log-Rank统计量，结果如表2所示，由表2结果可知，多组地下水储量的不同水平检验生存分析分布等同性。

表2.整体比较

	卡方	df	Sig.
				Log Rank(Mantel-Cox)	1101.223	1	.000
Breslow(Generalized Wilcoxon)	1093.332	1	.000
				Tarone-Ware	1101.002	1	.000

等比例风险模型(Cox Proportional Hazards Model)，可同时分析众多因素对地面沉降生存期的影响，对于地面沉降生存模型的确立具有重要意义，其具体步骤如下：

S51、确定Cox回归分析要求比例满足风险假定的前提条件，即假定风险比HR不随时间的变化而变化：

HR＝EXP{β}

其中，β为自变量的偏回归系数。

Cox模型由λ(t)表示的风险函数表示，风险函数λ(t)表示关注事件在t时刻发生的概率P(t，s)，即t时刻关注事件发生的风险，具体计算如下：

P(t，s)＝Pr(t＜T＜s|T≥t)

上述风险函数可大致分为指数模型(Exponential model)冈珀茨模型(Gompertzmodel)和威布尔模型(Weibull model)，统称为比例风险模型，其具体原理如下：

指数模型：

λ(t)＝λ

其中，λ为参数；

冈珀茨模型：

λ(t)＝λ₀e^αt whereλ₀＝e^μ

其中，λ和α为参数

威布尔模型：

λ(t)＝λ₀t^αwh ereα＞1(λ和α为参数)

在使用比例风险模型的同时，需要满足以下比例风险假定，具体风险函数形式如下所示：

λ(t，x)＝λ₀(t)h(x)＝λ₀(t)e^x′β

由于比例风险模型和加速失效时间模型均为参数模型，需要对风险函数的具体形式作出假设，再用最大似然法估计。但是截堵数据可能会导致风险函数设定错误，此时可应用等比例风险模型。它是一种半参数模型，利用个体i和个体j的风险函数相除，只剩下与解释变量x相关的项。

S52、以、λ(t，x)为因变量，得到等比例风险回归模型的基本形式为：

λ(t，x)＝λ₀(t)exp(β₁x₁+β₂x₂+…+β_mx_m)

其中，β为自变量的偏回归系数，λ(0)代表为零向量时，λ(t，x)的基准危险率；

A53、将i因素和j影响下的生存风险度相比，求得相对危险度RR：

本发明通过RADARSAT-2获取2012-2019年京津冀地区地表变形信息；利用空间统计查询方法整理得到沉降样本的生存时间和结局状态；随后利用Python对所有沉降样本进行生存分曲线绘制，得到时间序列下区域沉降发育差异；最后地下水储量、区域动载荷等影响因素，对沉降发育差异的风险影响因素进行揭示。本发明揭示并量化了地面沉降的时序发展过程和空间演化差异，进一步掌握了研究区地面沉降的时空变化差异，使得在应对地面沉降灾害时能够提出更加高效的防治措施，为城市灾害防治节约时间和经济成本，也为城市规划提供了科学依据。

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，根据本发明的技术方案及其改进构思加以等同替换或改变，都应涵盖在本发明的保护范围内。

Claims (1)

1.一种揭示地面沉降时空间演化差异的方法，其特征在于，包括以下步骤：

S1、获取覆盖研究区范围的长时序合成孔径雷达影像数据，并利用时序合成孔径雷达干涉测量方法，获取地表视线向的形变信息，然后通过实测GPS和水准点数据对InSAR监测结果进行校准和精度验证，将视线向的值转化为垂向，得到垂向形变数据，进而获取地面沉降速率；

S2、获取沉降生存指标：根据回归合成孔径雷达回归周期和研究区地面沉降发育程度，确定生存分析的时间单位，进而确定生存分析中的关注事件，即沉降生存指标；其中，所述关注事件以累计沉降值或沉降速率为阈值定义；

S3、根据区域沉降发育情况，确定起始事件I、终点事件E，计算得到每个沉降监测点的生存时间t，统计样本结局状态；

S4、根据研究区实际情况，通过Kaplan-Meier方法绘制生存函数；

S5、对研究区结合Log-Rank检验和等比例风险模型分析，确立分区方案以及地面沉降影响因子，构建区域地面沉降生存分析模型，揭示地面沉降时空间演化差异，量化展示地面沉降时序演化过程和空间分布差异；

步骤S4中，所述根据研究区实际情况，通过Kaplan-Meier方法绘制生存函数，具体包括以下步骤：

S41、对于具有n个样本的沉降数据集，根据合成孔径雷达回归周期T和研究区地面沉降发育程度，确定生存分析的时间单位t：首先将卫星回归周期T确定为分析时间单位t，计算T_n-T_n-1内地面形变速率，若在该T_n-T_n-1内地面形变速率内无法观察到毫米级形变，则将分析时间单位扩大至2T₀，以此类推；

S42、在步骤A1确定生存分析时间单位t的基础下，分别计算：t₁、t₂、t₃...t_n时刻的累积地面沉降值和当n＝1、2、3...n时的t_n-t_n-1内的沉降速率；

S43、根据研究区在研究时段的地面沉降发育情况，确定沉降速率v和累积沉降值s阈值，视为生存分析中的死亡事件，t_v和t_s视为死亡时刻，并分别统计每个样本的生存时间t_υ-t₀和t_s-t₀；

S44、根据可能的影响因素对沉降样本点进行分组，并结合样本的生存时间t_υ-t₀绘制每组的Kaplan-Meier生存曲线S_(t)，绘制方法如下：

其中，d为t_n时刻发生死亡事件的个体数，m为t_n时刻仍存活的样本数；

步骤S5中，所述Log-Rank检验，用于比较步骤S4中得到的两组生存曲线之间有无差别，具体步骤如下：

假设步骤S4中得到的两组生存曲线一样即零假设，基于此假设，得到两组内的事件发生个体数之比等于两组样本数之比，由此计算出事件发生的期望数，分别将两组所有时间点的期望数相加，与所有观察数进行比较；如果零假设成立，那么总的事件发生数d应该按样本量等比例分布在两组之间，那么每组期望数通过下式计算：

计算Log-Rank统计量：

根据得到的Log-Rank统计量数值判断是否拒绝零假设，进而判断两组生存曲线之间有统计学差异；

所述对研究区进行等比例风险模型分析，确立分区方案以及地面沉降影响因子，构建区域地面沉降生存分析模型，具体包括以下步骤：

S51、确定Cox回归分析要求比例满足风险假定的前提条件，即风险比不随时间的变化而变化：

HR＝EXP{β}

其中，HR代表即为假定风险比，β为自变量的偏回归系数；

S52、以λ(t，x)为因变量，等比例风险模型的基本形式为：

λ(t，x)＝λ₀(t)exp(β₁x₁+β₂x₂+...+β_mx_m)

其中，β为自变量的偏回归系数；λ(0)代表x为零向量时，λ(t，x)的基准危险率；

S53、因素i和j影响下的生存风险度相比，求得相对危险度RR：

Images