CN109828270B - 一种表征地面沉降时序演变的方法 - Google Patents

Abstract

本发明公开了一种表征地面沉降时序演变的方法，属于地质灾害领域。基于雷达影像数据，利用PS‑InSAR方法，获取区域地面沉降信息。根据沉降信息，将区域地面上的沉降时序演变特征进行分类，获取长时间序列下历史地面沉降量的发展趋势与类型，更加准确的掌握区域地面单元上的沉降的长时序演展规律，在应对地面沉降灾害时能够提出更加高效的防治措施，为城市灾害防治节约时间和经济成本，同时也可以作为城市规划的依据。

Description

一种表征地面沉降时序演变的方法

技术领域

本发明涉及地质灾害领域，尤其涉及一种表征地面沉降时序演变的方法。

背景技术

地面沉降是在自然因素和人为因素的影响下引发地表缓慢下降的一种环境地质现象。地面沉降会导致房屋、道路等基础设施的破坏，同时也可能引发水管、油气管的破裂、洪水频发和海水倒灌等一系列灾害。地面沉降造成的灾害会严重影响人类的生产和生活，阻碍经济的发展。对地面沉降进行长时序监测，掌握地面沉降在时序上的变化特征，由此分析地面沉降的演化规律，能够对地面沉降提出高效的防治措施。

目前，PS-InSAR技术已广泛应用于区域地面沉降的监测。为描述地面沉降在一定时期内的演化规律，现有研究大多仅对PS点多年平均沉降速率进行分析。实际上，某些区域地面沉降量在季节上的具有集中分布的特征，而且长时序上的沉降速率有时因为自然或者人为等因素发生变化，但是这些信息都无法通过PS-InSAR的方法获取。

针对上述问题，本发明提出一种新的区域地面沉降时序演变描述的方法。该方法用季节性强度来表示月平均沉降量的季节性，同时将平均沉降速率进行分类，归为匀速、加速和减速。最终根据重分类结果，可以准确掌握区域的地面沉降的时序演变特征。分析具有季节性的PS点，获取一年内地面沉降严重的时间，可以在相应的时间进行防控。根据变速的时间断点可以准确的掌握造成变速的原因，能够为地面沉降灾害的防治提供依据，对合理规划城市发展、预防和控制城市地面沉降具有重要的应用价值。

发明内容

本发明的目的是为了解决目前，现有的对地面沉降进行长时序监测，无法准确地掌握地面沉降在时序上的变化特征的问题，而提出的一种表征地面沉降时序演变的方法。

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

一种表征地面沉降时序演变的方法，包括以下步骤：

S1、获取覆盖研究区范围的长时序合成孔径雷达观测影像数据，假设影像获取时间为t₀，t₁，…t_n，利用永久散射体合成孔径雷达干涉测量PS-InSAR方法获取地表每一个PS点的形变信息，通过现场实测GPS和水准点数据对形变数据进行校准和精度验证，获取稳定PS点的高精度垂向形变数据库；

S2、将PS-InSAR结果进行预处理，将累积沉降量按照时间间隔进行插值，获取每月的累计沉降量X_m,将每相邻月份的沉降量进行相减，就获得了m个月的沉降量x₁，x₂...x_m，即x₁＝X₁-X₀，x₂＝X₂-X₁，…x_m＝X_m-X_m-1，再根据数据时间p，计算每月平均沉降量B_q；当数据跨度为p年时，p＝m/12,如果m/12不为整数，则，每月平均沉降量：B_q(X₁＝X₂₊₁₂+...X_1+12(p-1))/p；

S3、计算沉降的季节性，根据传统集中度的方法，用季节性强度M来表示地面沉降量集中于某几个月的程度，找到每一个PS点中B_q的最小值B_x，每个月的平均沉降量B_q加上B_x的绝对值∣B_x∣，获得变为正值的每月相对变化量C_q(_q＝1，2，…12)；其中，季节性强度M值在0-1之间，当M值在为1极大值时，说明沉降集中发生在某个月内，具有极强的季节性；当M为0极小值时，表示沉降量在1年内均匀的分布在12个月内；将M≥0.3表示沉降量的季节性显著,类型为季节性；反之，则表示季节性不显著，类型为非季节性；

S4、将沉降速率变化进行分类，将PS点将每月的累计沉降量X₀，X₁，…X_m分成两个部分，假设存在断点，断点为X_b(b≥5)，根据双线性模型公式计算时间断点X_b；

S5、使用ANOVAF检验法计算二次项回归系数为0的概率为P2，如果P2≤0.01，则二次项对回归具有明显的贡献，应被保留，时间序列也因此被分类为二次型；相反，若P2＞0.01，二次项对回归没有明显的贡献，时间序列被归类为线性；

S6、将S4中存在时间断点的数据判断数据是否存在纵向阶跃，将断点X_b前后的数据进行线性拟合，利用断点前线性方程进行预测，预测值为Y’_b+1…Y’_m，断点后的方程预测值为Y’₁…Y’_b，判断断点前后分段线性模型95％置信度的预测区间是否重叠，如果两个区间有重叠，则不存在纵向阶跃，时间序列被归类为双线性，将返回断点前后的斜率，b₀，b₁；反之，数据存在纵向阶跃；

S7、将S6中存在纵向阶跃的数据，利用方差分析方法，如果Pv>0.05表明断点两端数据之间斜率未发生变化，则被归类为不连续匀速类型，反之表示斜率发生变化则被归类为不连续变速并返回此类型前后的斜率b₀，b₁；

S8、将变化类型进行重分类，将不相关、线性、不连续匀速的PS点归为匀速；将二次型以及双线性和不连续变速中b₀>b₁的PS点归为减速；将二次型以及双线性和不连续变速中b₀<b₁的PS点归为加速；

S9、将累计沉降量时序变化分类，将S3中季节性的分类结果与S8的长时序变化分类结果进行叠加重分类，将区域地面累计表沉降量分布时序划分为季节性匀速、非季节性匀速、季节性减速、非季节性减速、季节性加速和非季节性加速等6类。

优选地，所述S1中雷达影像的处理软件包括常用的StamPS、Sarproz和GAMMA等软件。

优选地，所述S2中PS-InSAR结果的预处理过程，包括以下步骤：

A1、将累积沉降量S₀，S₁，…S_n按照时间间隔进行插值，假设共有m+1个月，获取每月的累计沉降量X₀，X₁，…X_m；

A2、将每相邻月份的沉降量进行相减，就获得了m个月的沉降量x₁，x₂…x_m，即x₁＝X₁-X₀，x₂＝X₂-X₁，…x_m＝X_m-X_m-1；

A3、根据数据时间，计算每月平均沉降量B_q(_q＝1，2，…12)。

优选地，所述S3中季节性强度的计算方法如下：首先计算出每个月相对变化量在全年变化量的比例r_q：

计算季节性强度M：

优选地，所述S2和S3中的计算通过ArcGIS、Matlab和Excel等软件计算所得。

优选地，所述S4中双线性模型的计算公式为：

Y_i＝γ(X_i)+ε_i,i＝1,…m

r(X_i)＝a+b₀[X₁I(X_i＜X_b)+X_bI(X_i≥X_b)]+b₁(X₁-X_b)I(X_i≥X_b)

其中，表示独立正态分布误差，b₀，b₁分段拟合的线性斜率；I为指示函数，当I后面的公式成立时，I＝1，反之，I＝0；根据贝叶斯信息准则BIC评价拟合优度，X_i处的BIC的计算公式为：

其中，SR²为残差平方和，RSS由SR²计算所得；如果BIC_min值小于线性回归BIC值，即BIC_L和二次型回归模型的BIC值，即BIC_Q，并且证据比B_w大于预定的阈值B_th，B_th≥1时，则说明存在显著的时间断点，反之，则不存在时间断点，其中B_w的计算公式如下：

其中，i＝1,2,3分别表示分段、线性和二次型三种回归模型，分别为分段、线性和二次型模型的权重，分别表示BIC_L、BIC_Q和BIC_min之间的差值(Δ₁＝0；Δ₂＝BIC_min；Δ₃＝BIC_Q-BIC_min)。

优选地，所述S6中断点X_b前后数据的线性公式：

Y_i＝a₀+b₀X_i X_i≤X_b

Y_i＝a₁+b₁X_i X_i＞X_b。

优选地，所述S4和S6中的计算通过Matlab、Python和C++等计算获得。

与现有技术相比，本发明提供了一种表征地面沉降时序演变的方法，具备以下有益效果：

(1)本发明利用雷达影像数据，基于雷达影像数据，利用PS-InSAR方法，获取区域地面沉降信息。根据沉降信息，将区域地面上的沉降时序演变特征进行分类，获取长时间序列下历史地面沉降量的发展趋势与类型，更加准确的掌握区域地面单元上的沉降的长时序演展规律，在应对地面沉降灾害时能够提出更加高效的防治措施，为城市灾害防治节约时间和经济成本，同时也可以作为城市规划的依据。

(2)本发明中精确描述地表形变的时间序列变化特征，并且处理的PS点数量庞大，可以很好的扩展到区域尺度，能够掌握其空间上的分布特征。结合时序变化特征和空间分布状况，能够为依此建立了模型和数据库，可以有效的对未来地面沉降变化进行预测。对灾害集中区域加强治理，减少灾害造成的损失。

(3)本发明采用季节性强度来表示月平均沉降量的季节性，同时将平均沉降速率进行分类，归为匀速、加速和减速。最终根据重分类结果，可以准确掌握区域的地面沉降的时序演变特征。分析具有季节性的PS点，获取一年内地面沉降严重的时间，可以在相应的时间进行防控。根据变速的时间断点可以准确的掌握造成变速的原因，能够为地面沉降灾害的防治提供依据，对合理规划城市发展、预防和控制城市地面沉降具有重要的应用价值。

附图说明

图1为本发明提出的一种表征地面沉降时序演变的方法的流程示意图；

图2为本发明提出的一种表征地面沉降时序演变的方法的PS-InSAR流程图；

图3为本发明提出的一种表征地面沉降时序演变的方法的沉降速率分类流程图。

具体实施方式

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

在本发明的描述中，需要理解的是，术语上、下、前、后、左、右、顶、底、内、外等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。

实施例1：一种表征地面沉降时序演变的方法，其特征在于，包括以下步骤：

S1、获取覆盖研究区范围的长时序合成孔径雷达观测影像数据，假设影像获取时间为t₀，t₁，…t_n，利用永久散射体合成孔径雷达干涉测量PS-InSAR方法获取地表每一个PS点的形变信息，通过现场实测GPS和水准点数据对形变数据进行校准和精度验证，获取稳定PS点的高精度垂向形变数据库；

S2、将PS-InSAR结果进行预处理，将累积沉降量按照时间间隔进行插值，获取每月的累计沉降量X_m,将每相邻月份的沉降量进行相减，就获得了m个月的沉降量x₁，x₂...x_m，即x₁＝X₁-X₀，x₂＝X₂-X₁，…x_m＝X_m-X_m-1，再根据数据时间p，计算每月平均沉降量B_q；当数据跨度为p年时，p＝m/12,如果p＝m/12不为整数，则每月平均沉降量：B_q(X₁＝X₂₊₁₂+...X_1+12(p-1))/p；

S3、计算沉降的季节性，根据传统集中度的方法，用季节性强度M来表示地面沉降量集中于某几个月的程度，找到每一个PS点中B_q的最小值B_x，每个月的平均沉降量B_q加上B_x的绝对值∣B_x∣，获得变为正值的每月相对变化量C_q(_q＝1，2，…12)；其中，季节性强度M值在0-1之间，当M值在为1极大值时，说明沉降集中发生在某个月内，具有极强的季节性；当M为0极小值时，表示沉降量在1年内均匀的分布在12个月内；将M≥0.3表示沉降量的季节性显著,类型为季节性；反之，则表示季节性不显著，类型为非季节性；

S4、将沉降速率变化进行分类，将PS点将每月的累计沉降量X₀，X₁，…X_m分成两个部分，假设存在断点，断点为X_b(b≥5)，根据双线性模型公式计算时间断点X_b；

S5、使用ANOVAF检验法计算二次项回归系数为0的概率为P2，如果P2≤0.01，则二次项对回归具有明显的贡献，应被保留，时间序列也因此被分类为二次型；相反，若P2＞0.01，二次项对回归没有明显的贡献，时间序列被归类为线性；

S6、将S4中存在时间断点的数据判断数据是否存在纵向阶跃，将断点X_b前后的数据进行线性拟合，利用断点前线性方程进行预测，预测值为Y’_b+1…Y’_m，断点后的方程预测值为Y’₁…Y’_b，判断断点前后分段线性模型95％置信度的预测区间是否重叠，如果两个区间有重叠，则不存在纵向阶跃，时间序列被归类为双线性，将返回断点前后的斜率，b₀，b₁；反之，数据存在纵向阶跃；

S7、将S6中存在纵向阶跃的数据，利用方差分析方法，如果Pv>0.05表明断点两端数据之间斜率未发生变化，则被归类为不连续匀速类型，反之表示斜率发生变化则被归类为不连续变速并返回此类型前后的斜率b₀，b₁；

S8、将变化类型进行重分类，将不相关、线性、不连续匀速的PS点归为匀速；将二次型以及双线性和不连续变速中b₀>b₁的PS点归为减速；将二次型以及双线性和不连续变速中b₀<b₁的PS点归为加速；

S9、将累计沉降量时序变化分类，将S3中季节性的分类结果与S8的长时序变化分类结果进行叠加重分类，将区域地面累计表沉降量分布时序划分为季节性匀速、非季节性匀速、季节性减速、非季节性减速、季节性加速和非季节性加速等6类。

S1中雷达影像的处理软件包括常用的StamPS、Sarproz和GAMMA等软件。

实施例2：基于实施例1但有所不同的是，S2中PS-InSAR结果的预处理过程，包括以下步骤：

A1、将累积沉降量S₀，S₁，…S_n按照时间间隔进行插值，假设共有m+1个月，获取每月的累计沉降量X₀，X₁，…X_m；

A2、将每相邻月份的沉降量进行相减，就获得了m个月的沉降量x₁，x₂…x_m，即x₁＝X₁-X₀，x₂＝X₂-X₁，…x_m＝X_m-X_m-1；

A3、根据数据时间，计算每月平均沉降量B_q(_q＝1，2，…12)。

本发明利用雷达影像数据，基于雷达影像数据，利用PS-InSAR方法，获取区域地面沉降信息。根据沉降信息，将区域地面上的沉降时序演变特征进行分类，获取长时间序列下历史地面沉降量的发展趋势与类型，更加准确的掌握区域地面单元上的沉降的长时序演展规律，在应对地面沉降灾害时能够提出更加高效的防治措施，为城市灾害防治节约时间和经济成本，同时也可以作为城市规划的依据。

实施例3：基于实施例1和2但有所不同的是，本方案先通过Sarproz软件处理2010-2016年的Terrasar-X雷达遥感影像，利用水准点数据验证PS-InSAR数据结果的精度，获取高精度地面沉降信息。将北京区域作为研究区，利用监测方法获得地面沉降季节性强度和平均沉降速率的变化，分析2010-2016年北京地区地面沉降的时序演变规律。主要步骤如下：

步骤一：将2010-2016年覆盖北京平原区的53景Terrasar-X数据利用Sarproz数据处理软件进行处理，获取北京地区2010-2016年的地面沉降信息，利用研究区内六个水准点监测数据对PS-InSAR结果进行精度验证，获取高精度的地面沉降信息。

步骤二：将PS-InSAR中的累计沉降量结果进行插值，获得2010年5月1-2016年6月1日的累计沉降量，通过计算，获得每个月的沉降量以及每月平均沉降量。首先根据季节性强度的计算公式，得到每个PS点的季节性强度，将所有的PS点分为季节性和非季节性两类。

步骤三：将2010年5月1-2016年6月1日的累计沉降量利用MATLAB通过ANOVA F检验，BIC评价等方法将沉降速率的变化进行分类，主要分成不相关、线性、二次型、双线性、不连续匀速和不连续变速六种类型。变化类型进行重分类，将不相关、线性、不连续匀速的PS点归为匀速；将二次型以及双线性和不连续变速中b₀>b₁的PS点归为减速；将二次型以及双线性和不连续变速中b₀<b₁的PS点归为加速。

步骤四：利用ArcGIS根据步骤二和步骤三将点进行叠加分析，将PS点进行重分类，共分为6类：季节性匀速，非季节性匀速，季节性减速，非季节性减速季节性加速非季节性加速。

步骤五：根据上述分类结果。利用ArcGIS显示，可以在区域空间上找到每种类型的分布状况。尤其是季节性变速的类型，根据其季节性特征，分析沉降主要集中的季节，以及断点时间，变速的状况，可以分析其变速原因，为地面沉降灾害防治提供依据。

本发明中精确描述地表形变的时间序列变化特征，并且处理的PS点数量庞大，可以很好的扩展到区域尺度，能够掌握其空间上的分布特征。结合时序变化特征和空间分布状况，能够为依此建立了模型和数据库，可以有效的对未来地面沉降变化进行预测。对灾害集中区域加强治理，减少灾害造成的损失。

实施例4：基于实施例1、2和3但有所不同的是，S3中季节性强度的计算方法如下：首先计算出每个月相对变化量在全年变化量的比例r_q：

计算季节性强度M：

S2和S3中的计算通过ArcGIS、Matlab、Excel等软件计算所得。

S4中双线性模型的计算公式为：

Y_i＝γ(X_i)+ε_i,i＝1,…m

r(X_i)＝a+b₀[X₁I(X_i＜X_b)+X_bI(X_i≥X_b)]+b₁(X₁-X_b)I(X_i≥X_b)

其中，表示独立正态分布误差，b₀，b₁分段拟合的线性斜率；I为指示函数，当I后面的公式成立时，I＝1，反之，I＝0；根据贝叶斯信息准则BIC评价拟合优度，X_i处的BIC的计算公式为：

其中，SR²为残差平方和，RSS由SR²计算所得；如果BIC_min值小于线性回归BIC值，即BIC_L和二次型回归模型的BIC值，即BIC_Q，并且证据比B_w大于预定的阈值B_th，B_th≥1时，则说明存在显著的时间断点，反之，则不存在时间断点，其中B_w的计算公式如下：

其中，i＝1,2,3分别表示分段、线性和二次型三种回归模型，分别为分段、线性和二次型模型的权重，分别表示BIC_L、BIC_Q和BIC_min之间的差值(Δ₁＝0；Δ₂＝BIC_min；Δ₃＝BIC_Q-BIC_min)。

S6中断点X_b前后数据的线性公式：

Y_i＝a₀+b₀X_i X_i≤X_b

Y_i＝a₁+b₁X_i X_i＞X_b。

S4和S6中的计算通过Matlab、Python、C++等计算获得。

本发明采用季节性强度来表示月平均沉降量的季节性，同时将平均沉降速率进行分类，归为匀速、加速和减速。最终根据重分类结果，可以准确掌握区域的地面沉降的时序演变特征。分析具有季节性的PS点，获取一年内地面沉降严重的时间，可以在相应的时间进行防控。根据变速的时间断点可以准确的掌握造成变速的原因，能够为地面沉降灾害的防治提供依据，对合理规划城市发展、预防和控制城市地面沉降具有重要的应用价值。

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

Claims (8)

1.一种表征地面沉降时序演变的方法，其特征在于，包括以下步骤：

S1、获取覆盖研究区范围的长时序合成孔径雷达观测影像数据，假设影像获取时间为t₀，t₁，…t_n，利用永久散射体合成孔径雷达干涉测量PS-InSAR方法获取地表每一个PS点的形变信息，通过现场实测GPS和水准点数据对形变数据进行校准和精度验证，获取稳定PS点的高精度垂向形变数据库；

S2、将PS-InSAR结果进行预处理，将累计沉降量按照时间间隔进行插值，获取每月的累计沉降量X_m,将每相邻月份的沉降量进行相减，就获得了m个月的沉降量x₁，x₂...x_m，即x₁＝X₁-X₀，x₂＝X₂-X₁，…x_m＝X_m-X_m-1，再根据数据时间p，计算每月平均沉降量B_q；当数据跨度为p年时，p＝m/12，p＝m/12不为整数时，p＝(m/12)+1，则每月平均沉降量：B_q＝(X₁+X₁₊₁₂+...X_1+12(p-1))/p；

S3、计算沉降的季节性，根据传统集中度的方法，用季节性强度M来表示地面沉降量集中于某几个月的程度，找到每一个PS点中B_q的最小值B_x，每个月的平均沉降量B_q加上B_x的绝对值∣B_x∣，获得变为正值的每月相对变化量C_q(_q＝1，2，…12)；其中，季节性强度M值在0-1之间，当M值在为1极大值时，说明沉降集中发生在某个月内，具有极强的季节性；当M为0极小值时，表示沉降量在1年内均匀的分布在12个月内；将M≥0.3表示沉降量的季节性显著,类型为季节性；反之，则表示季节性不显著，类型为非季节性；

S4、将沉降速率变化进行分类，将PS点将每月的累计沉降量X₀，X₁，…X_m分成两个部分，假设存在断点，断点为X_b(b≥5)，根据双线性模型公式计算时间断点X_b；

S5、使用ANOVAF检验法计算二次项回归系数为0的概率为P2，如果P2≤0.01，则二次项对回归具有明显的贡献，应被保留，时间序列也因此被分类为二次型；相反，若P2＞0.01，二次项对回归没有明显的贡献，时间序列被归类为线性；

S6、将S4中存在时间断点的数据判断数据是否存在纵向阶跃，将断点X_b前后的数据进行线性拟合，利用断点前线性方程进行预测，预测值为Y’_b+1…Y’_m，断点后的方程预测值为Y’₁…Y’_b，判断断点前后分段线性模型95％置信度的预测区间是否重叠，如果两个区间有重叠，则不存在纵向阶跃，时间序列被归类为双线性，将返回断点前后的斜率，b₀，b₁；反之，数据存在纵向阶跃；

S7、将S6中存在纵向阶跃的数据，利用方差分析方法，如果Pv>0.05表明断点两端数据之间斜率未发生变化，则被归类为不连续匀速类型，反之表示斜率发生变化则被归类为不连续变速并返回此类型前后的斜率b₀，b₁；

S8、将变化类型进行重分类，将线性、不连续匀速的PS点归为匀速；将二次型以及双线性和不连续变速中b₀>b₁的PS点归为减速；将二次型以及双线性和不连续变速中b₀<b₁的PS点归为加速；

S9、将累计沉降量时序变化分类，将S3中季节性的分类结果与S8的长时序变化分类结果进行叠加重分类，将区域地面累计表沉降量分布时序划分为季节性匀速、非季节性匀速、季节性减速、非季节性减速、季节性加速和非季节性加速6类。

2.根据权利要求1所述的一种表征地面沉降时序演变的方法，其特征在于：所述S1中雷达影像的处理软件包括常用的StamPS、Sarproz和GAMMA这些软件。

3.根据权利要求1所述的一种表征地面沉降时序演变的方法，其特征在于：所述S2中PS-InSAR结果的预处理过程，包括以下步骤：

A1、将累计沉降量S₀，S₁，…S_n按照时间间隔进行插值，假设共有m+1个月，获取每月的累计沉降量X₀，X₁，…X_m；

A2、将每相邻月份的沉降量进行相减，就获得了m个月的沉降量x₁，x₂…x_m，即x_1＝X₁-X₀，x_2＝X₂-X₁，…x_m＝X_m-X_m-1；

A3、根据数据时间，计算每月平均沉降量B_q(_q＝1，2，…12)。

4.根据权利要求1所述的一种表征地面沉降时序演变的方法，其特征在于：所述S3中季节性强度的计算方法如下：首先计算出每个月相对变化量在全年变化量的比例r_q：

计算季节性强度M：

5.根据权利要求1所述的一种表征地面沉降时序演变的方法，其特征在于：所述S2和S3中的计算通过ArcGIS、Matlab、Excel这些软件计算所得。

6.根据权利要求1所述的一种表征地面沉降时序演变的方法，其特征在于：所述S4中双线性模型的计算公式为：

Y_i＝γ(X_i)+ε_i,i＝1,…m

r(X_i)＝a+b₀[X₁I(X_i＜X_b)+X_bI(X_i≥X_b)]+b₁(X₁-X_b)I(X_i≥X_b)

其中，表示独立正态分布误差，b₀，b₁分段拟合的线性斜率；I为指示函数，当I后面的公式成立时，I＝1，反之，I＝0；根据贝叶斯信息准则BIC评价拟合优度，X_i处的BIC的计算公式为：

其中，SR²为残差平方和，RSS由SR²计算所得；如果BIC_min值小于线性回归BIC值，即BIC_L和二次型回归模型的BIC值，即BIC_Q，并且证据比B_w大于预定的阈值B_th，B_th≥1时，则说明存在显著的时间断点，反之，则不存在时间断点，其中B_w的计算公式如下：

其中，i＝1,2,3分别表示分段、线性和二次型三种回归模型，分别为分段、线性和二次型模型的权重，分别表示BIC_L、BIC_Q和BIC_min之间的差值(Δ₁＝0；Δ₂＝BIC_min；Δ₃＝BIC_Q-BIC_min)。

7.根据权利要求1所述的一种表征地面沉降时序演变的方法，其特征在于：所述S6中断点X_b前后数据的线性公式：

Y_i＝a₀+b₀X_i X_i≤X_b

Y_i＝a₁+b₁X_i X_i＞X_b。

8.根据权利要求1所述的一种表征地面沉降时序演变的方法，其特征在于：所述S4和S6中的计算通过Matlab、Python、C++的计算获得。

Images