联合DInSAR的3种下沉时序模型关键点缺失问题研究

石晓宇，魏祥平，杨可明，王 剑，姚树一

(1.中国矿业大学(北京) 地球科学与测绘工程学院，北京 100083；2.中国四维测绘技术有限公司，北京 100094；3. 淮北矿业股份有限公司通防地测部，安徽 淮北 235000)

0 引 言

地下煤炭资源的大规模开采，会危及采空区上方及周边基础建(构)筑物的安全，同时也会带来一系列新的衍生问题[1]。因此，在开展地下采煤中，为控制采动损害等级和保护采区周围各类建(构)筑物免受或少受采动损害影响，需实时监测并提前预估煤炭开采导致的地表变形情况，以便分析和掌握矿区地表动态沉降规律。

目前，传统的各类开采沉陷预计方法只是针对地下采矿引起的地表沉陷现象进行研究(不考虑全采工作面上覆岩层移动机理和过程)[2]，所以此类方法仅能较好地描述地表受采动影响稳定后的沉陷状态，并不能描述地表沉陷的动态过程。因此，众多学者以局部观测点作为切入点，建立沉陷区单点下沉量时间序列模型，从而描述矿区地表单点下沉随时间变化的趋势。已有的下沉时序模型包括Knothe模型、双曲线模型、Gompertz模型、幂指数Knothe模型、Logistic模型和Weibull模型等，其中幂指数Knothe模型、Logistic模型和Weibull模型是3种常用的非线性时间函数模型，已有较多文献研究表明3种时间函数模型均能较好地描述矿山开采沉陷过程[3-5]。此外，传统监测方式是以水准测量为代表，而近年来合成孔径雷达差分干涉测量(Differential Interferometric Synthetic Aperture Radar，DInSAR)技术成为了新兴的矿区地表形变监测手段[6]，2种监测方法瑕瑜互见。前者监测精度较高，但易受观测环境以及观测人员的影响,获取的监测数据周期不稳定且时间跨度较大，从而致使矿区地表点动态沉降过程中的某些“关键点”(地表单点下沉曲线的拐点)未被测量，后者由于其固定的重访周期以及特别时间段(夏季)失相干严重等问题也会导致地表点下沉时间序列中某些关键点无监测数据。当前众多研究[7-10]主要基于单一的水准测量数据建立时间函数，分析矿区地表点动态形变趋势，然而某些时段“关键点”缺失，会导致预估的下沉曲线偏离真实下沉状态，导致不能准确地建立下沉时间序列模型来分析地表点动态沉降规律。因此，若联合水准测量与DInSAR技术，基于多源监测手段获取多源监测数据实现时间序列的加密，建立相应的预计时间函数，可提高拟合曲线与真实下沉曲线的吻合度，进一步提高沉陷预计精度。

为克服以上不足，选取幂指数Knothe时间模型、Logistic模型和Weibull模型，首先以仿真数据分析模型参数大小以及关键点缺失对3种模型拟合精度影响；之后，以淮北袁二煤矿某一采区为试验对象，采用DInSAR技术监测采区阶段的地表形变，获得各水准监测点2017-11-04—2018-02-08的InSAR时序监测值；最后，以InSAR监测值作为水准测量加密值，整合InSAR与水准数据作为原始观测序列，建立幂指数Knothe时间模型、Logistic模型和Weibull模型，结合矿区地表单点动态沉降规律，对比并分析基于多源和单一数据分别求解的3种模型的拟合结果差异。

1 沉陷区地表单点下沉特征与时间序列模型

1.1 沉陷区地表观测点下沉变形特征

地下采煤引起的采空区地表动态沉陷变形是一个随时间和空间变化的复杂四维问题[11]，随着回采工作面的推进，以下沉速度为划分大致可分为3个时段：下沉开始阶段(下沉速度v1≤1.67 mm/d)、下沉活跃阶段(下沉速度v2>1.67 mm/d)和下沉衰退阶段(下沉速度v3≤1.67 mm/d)[12]。相应地，3个时段也反映了沉陷区地表单点下沉量Wt与时间t函数曲线变化特征，下沉量随时间变化的曲线形态如图1所示。图1中，矿区地表单点下沉曲线最终近似呈“S型”，每一阶段的跃迁都存在相应的拐点，这3个拐点由左往右分别为正最大加速度点、最大速度点和负最小加速度点[13]，与下沉曲线的方向、斜率以及变化率密切相关，不仅对沉降曲线形态起约束作用，而且也是常用时间序列模型是否能较好拟合观测点下沉动态过程的3个关键点。

Wmax—累计最大下沉量

1.2 3种常用下沉量时间函数模型

现有文献[14-15]中已有较多代表性的时间函数模型被应用于研究矿区地表沉陷动态过程，常用的3种时间函数模型各有千秋，分别为幂指数Knothe时间模型、Logistic模型和Weibull模型。

1.2.1 幂指数Knothe时间模型

Knothe模型是波兰学者Knothe于1952年提出的可用于描述地表动态下沉的时间函数模型[16]，幂指数Knothe模型[17]为修正的Knothe时间函数模型，其基本函数表达式为

W(t)=Wmax(1-e-bt)m

(1)

式中:W(t)为t时刻某地表沉陷点的瞬时下沉量；Wmax为地表沉陷稳定后达到的最大下沉量；b和m为Knothe模型待估参数，定义为影响地表下沉规律且与采场上覆岩层特性有关的参数。

对式(1)时间函数进行一阶和二阶求导运算，可得到下沉速度v以及下沉加速度a关于时间t的函数，函数表达式分别为

v(t)=Wmaxbme-bt(1-e-bt)m-1

(2)

a(t)=Wmax[-mb2e-bt(1-e-bt)m-1+m(m-1)b2e-2bt(1-e-bt)m-2]

(3)

式中：v(t)为t时刻某地表沉陷点的瞬时下沉速度；a(t)为t时刻某地表沉陷点的瞬时下沉加速度。

1.2.2 Logistic模型

Logistic模型是荷兰数学家Verhulst于1938年最早提出的生物增长模型[18]，其“S型”增长曲线所描述的某指标时间序列变化过程与矿区地表沉陷单点下沉的物理过程较相似。下沉量函数表达式W(t)、下沉速度时间函数v(t)以及下沉加速度时间函数a(t)分别为

(4)

(5)

(6)

式中：c和d为Logistic模型待估参数，定义为影响地表下沉规律且与采场上覆岩层特性有关的参数。

1.2.3 Weibull模型

Weibull模型是韦布尔于1951年提出的统计分布模型[19]，其时间函数可较准确地拟合和描述沉陷区地表观测点的下沉曲线形态。下沉量时间函数W(t)、下沉速度时间函数v(t)以及下沉加速度时间函数a(t)的表达形式分别为

W(t)=Wmax(1-e-ftk)

(7)

v(t)=Wmaxfktk-1e-ftk

(8)

a(t)=Wmaxfke-ftk[(k-1)tk-2-t2(k-1)fk]

(9)

式中：f和k为Weibull模型待估参数，定义为影响地表下沉规律且与采场上覆岩层特性有关的参数。

2 影响各时间函数模型曲线拟合的因素分析

2.1 时间序列模型地表下沉待估参数分析

1)幂指数Knothe模型待估参数b和m。以模拟值研究幂指数Knothe模型中待估参数b和m对模型曲线的影响，首先假设式(1)中Wmax=-2 000 mm，t=0～100 d，m=10，则当b取值为0.08，0.1，0.15，0.2时，幂指数Knothe模型随待估参数b变化的下沉曲线、下沉速度曲线以及下沉加速度曲线如图2a—图2c所示;同理假设式(1)中Wmax=-2 000 mm，t=0～100 d，b=0.1，当m取值为1、5、10、100时，幂指数Knothe模型随待估参数m变化的下沉曲线、下沉速度曲线以及下沉加速度曲线如图2d—图2f所示。由图2可知，随着参数b的增大，下沉曲线越来越陡，地表单点开始下沉至衰减阶段所需时间越来越短，最大下沉速度与加速度推前且大小变化明显；相反，待估参数m的增大只会致使地表下沉达到最大速度的时间迟延，而最大下沉速度与加速度的大小相差无几。

图2 幂指数Knothe模型待估参数对下沉量、速度和加速度的影响

2)Logistic模型待估参数c和d。以模拟值研究Logistic模型中待估参数c和d对模型曲线的影响，首先假设式(4)中Wmax=-2 000 mm，t=0～200 d，c=1 000，则当d取值为0.08，0.1，0.15，0.2时，Logistic模型随待估参数d变化的下沉曲线、下沉速度曲线以及下沉加速度曲线如图3a—图3c所示;同理假设式(4)中Wmax=-2 000 mm，t=0～200 d，d=0.1，当c取值为1 000，5 000，8 000，10 000时，Logistic模型随待估参数c变化的下沉曲线、下沉速度曲线以及下沉加速度曲线如图3d—图3f所示。由图3中可以看出，参数d对各曲线(图3a—图3c)的形态以及数值上影响较显著，最大下沉速度和加速度均存在时间滞后；而参数c的变化只会致使各曲线(图3d—图3f)发生平移变换而无伸缩变换。

图3 Logistic模型待估参数对下沉量、速度和加速度的影响

3)Weibull模型待估参数f和k。以模拟值研究Weibull模型中待估参数f和k对模型曲线的影响，首先假设式(7)中Wmax=-2 000 mm，t=0～200 d，k=2.1，则当f取值为0.95×10-4，3×10-4，5×10-4，1×10-3时，Logistic模型随待估参数f变化的下沉曲线、下沉速度曲线以及下沉加速度曲线如图4a—图4c所示;同理假设式(7)中Wmax=-2 000 mm，t=0～200 d，f=1×10-4，当k取值为2，2.1，2.2，2.5时，Weibull模型随待估参数k变化的下沉曲线、下沉速度曲线以及下沉加速度曲线如图4d—图4f所示。由图4中可以看出，随着参数f的减小，地表单点沉陷的三阶段时期延长，沉陷最大速度和加速度也随之减小；随着k的增大，地表下沉达到最大下沉速度和加速度的时间会推前，沉陷稳定所需时间缩短。

图4 Weibull模型待估参数对下沉量、速度和加速度的影响

综上所述，各时间序列模型地表下沉待估参数的大小对曲线拟合形态有较大影响，需要选择合适的下沉参数才能更准确地拟合沉陷区观测点下沉时间序列。因此，综合考虑到3种时间函数模型的非线性特征以及算法求解难易程度，选用LM(Levenberg-Marquard)[20]算法求解非线性最小二乘问题，求得待估参数最优值。

2.2 关键点缺失对各时间函数拟合的影响

依据上述模拟数据分析关键点缺失对3种时间函数曲线拟合的影响，首先假定不同监测周期情况下，将各组模拟数据中加入正态随机误差，使其模拟环境更接近真实情况；然后采用LM算法求解各待估参数，以建立时间序列模型，最终将得到的最佳下沉参数拟合的曲线结果与模拟曲线进行对比，分析关键点缺失对曲线函数拟合的影响。

幂指数Knothe模型不同监测周期下(20、50 d)的拟合结果(图5)。

1—模拟曲线；2—拟合曲线；—观测点

Logistic模型不同监测周期(20、50、60、75 d)下的拟合结果如图6所示；Weibull模型不同监测周期(20、60、75 d)下的拟合结果如图7所示。

1—模拟曲线；2—拟合曲线；—观测点

1—模拟曲线；2—拟合曲线；—观测点

结合图5—图7可发现，当观测点几乎涵盖曲线上任意关键点且较均匀分布在动态下沉曲线上时(图5a、图6a以及图7a)，监测周期为20 d的各模型拟合曲线与模拟曲线吻合度较高(最大均方根差约为13 mm)；然而，随着监测周期增大(图5b、图6b和6c)，观测点未覆盖曲线特征点处，拟合曲线开始偏离模拟曲线，吻合度和精度相比监测周期20 d的拟合结果均较低，最大均方根差达293 mm；此外，对比图6b和6d发现，监测周期为75 d的观测点个数虽小于观测周期为50 d的监测次数，但前者相对于后者涵盖的特征点位置明显优于后者，所以其拟合结果也明显优于后者，就均方根误差而言，减少了7 mm(75 d的均方根误差约为53 mm；50 d的均方根误差约为60 mm)。综上所述，矿区地表沉降监测周期是否涵盖下沉过程的关键点对时间函数模型的拟合结果有较大影响，因此，联合多源监测数据加密观测点监测周期，可提高各下沉时间函数曲线的拟合精度。

3 联合多源监测数据的矿区下沉时序模型试验

3.1 袁二矿区DInSAR监测

袁二煤矿位于安徽省涡阳县内，采区7221工作面上覆较多的人造水渠和农田，工作面走向长544 m，倾向长110 m，采高达4.5 m，于2017-12-06开始回采。在开采工作面的同时，淮北矿业集团建立岩移观测站并布置走向倾向观测线(水准首次联测时间为2017-10-08，监测周期不等)，研究7221工作面开采对采区上方地表沉降变形以及周边建构筑物的影响。然而，这些水准点受地理位置以及自然环境的影响，监测数量、监测范围和监测周期都相对有限，因此，需要在该煤矿工作面回采阶段进行简便有效的DInSAR监测。选取袁二煤矿7221工作面开采期间9景C波段的哨兵1 A影像图(第1幅雷达影像为2017-11-04，重访周期为12 d)，共形成8个干涉对，影像下载详细参数见表1，干涉对详细参数见表2。

表1 影像下载参数

表2 干涉对详细参数

结合研究区DEM数据(SRTM90)，将8对干涉组合分别进行二轨法DInSAR差分处理，获取了7221工作面2017-11-04—2018-02-08地表时序沉降形变，2017-11-04—2018-02-08的地表累计沉降形变如图8所示，并从InSAR时序形变图中提取工作面走向的剖面下沉序列曲线，如图9所示。为验证DInSAR技术测量的精度，将水准各观测线路与DInSAR监测的结果进行对比，由于两者时间基准不一致，故通过对水准测量时序值进行非线性插值拟合，获得干涉对研究时段的沉降值，以尽可能减少因时间间隔不统一而带来的误差，最终得到RMSE为0.017 7 m，水准点H2处的相对误差为4.02%，水准点G7处的相对误差为6.80%。

图8 地表沉降形变

图9 工作面走向剖面下沉时序图

3.2 联合InSAR与水准数据建立下沉时序模型

为验证基于多源监测数据建立的下沉时序模型，可提高曲线拟合精度，以上述基于DInSAR技术获取相应水准点的InSAR时序监测数据(2017-11-04—2018-02-08，监测周期12 d)和岩移观测站联测的水准数据(2017-10-08—2018-07-22，前期监测周期为1个月，后期监测周期为7 d)为联合对象，分别建立幂指数Knothe时间模型、Logistic模型和Weibull模型，并分析各曲线拟合结果。在SAR影像成像期间内，选取的试验水准点主要位于下沉开始和活跃阶段(结合图9)，因沿工作面走向剖面的北部各监测点沉降量均很小，并不能准确描述地表单点下沉三阶段的变形特征，故从获取的7221工作面地表沉降时序形变图上，选取邻近下沉盆地中心H观测线上的水准点H2和下沉盆地边沿G观测上的水准点G7作为试验对象(具体位置如图8所示)，分别基于联合多源监测数据(水准和InSAR数据)和仅基于单一监测数据(水准数据)建立了3种下沉量、下沉速度时间函数模型，2种情况的对比结果如图10、图11所示。

1—单一水准数据拟合曲线；2—多源数据拟合曲线；☆—InSAR加密点；+—水准点

1—单一水准数据拟合曲线；2—多源数据拟合曲线

3.3 结果对比与分析

结合图10和图11试验水准点H2和G7建立的3种时间函数模型拟合结果作以下对比并分析：

1) 无论基于多源或单一监测数据而言，幂指数Knothe时间模型、Logistic模型和Weibull模型均能较好地拟合地表单点下沉的时间序列，以及较准确地描述沉陷区地表单点动态下沉三阶段(开始阶段、活跃阶段和衰退阶段)，最大均方根误差约为5.941 mm。下沉盆地中心试验水准点H2下沉序列使用Logistic模型拟合效果最佳；下沉边沿水准点G7下沉序列使用幂指数Knothe时间模型拟合均方根误差最小。

2) 图10a、10b中，水准点H2下沉时间序列分别建立幂指数Knothe时间模型、Logistic模型，各时间函数模型拟合的曲线1和曲线2在时间段100～300 d上均无较大差异，吻合度较高，在时间段0～50 d上，由于加入的4期InSAR加密点，使得均方根误差相较于单一水准拟合分别减少了10.22%和8.11%。而图10c中，Weibull模型受基于多源和单一数据拟合结果影响较小，仅减少了3.35%。图11中，曲线1和曲线2相比，虽最大沉降速度极值出现无明显差异，但出现最大沉降速度极值增大等问题，不利于评估地表动态破坏。试验对象点H2的监测周期较好地涵盖了下沉三阶段中的各关键点，各水准点在下沉曲线上呈均匀分布状态，并未对各时间函数模型的曲线拟合造成较大影响。

3) 图10d、图10e中，水准点G7下沉时间序列基于Logistic和Weibull时间函数模型拟合的曲线1和曲线2吻合度较低，自时间轴上100 d后，二者出现一定偏离。其主要原因是：试验对象点G7的监测周期变化较大(最长监测周期达68 d)，仅基于单一水准数据拟合的曲线1均出现最大下沉值、最大速度时间迟延现象(图10、11e和11f)，同样不利于评估地表受采动损害影响；相反，对于联合多源数据拟合的曲线2而言，加入的后3期InSAR加密点(2018-01-15—2018-02-08)对各模型的下沉曲线活跃阶段起到了较好的约束作用，拟合结果相比曲线1更确实。图10d中，幂指数Knothe模型相对于其他模型受关键点缺失影响较小，联合InSAR与水准的均方根误差较单一水准仅减少了3.38%。

4 结 论

1) 针对难以全面研究沉陷区下沉变形的动态过程和单一监测手段可能会错失地表沉降阶段某些关键点的问题，通过模拟数据分析下沉待估参数以及关键点缺失对幂指数Knothe时间模型、Logistic模型和Weibull模型3种时间函数曲线拟合的影响，发现随着监测周期间隔越长，曲线拟合精度越低，保证一定的监测周期，均方根误差最大可提高95.56%。

2) 沉陷区观测点监测周期是否涵盖动态沉降中系列关键点对幂指数Knothe时间模型、Logistic模型和Weibull模型拟合精度有较大影响，将InSAR监测数据作为加密点与水准监测数据整合，增加多余观测数据，生成新下沉时间序列建立下沉时间序列模型，可较好地约束模型曲线形态，提高曲线拟合精度和分析矿区地表动态沉降规律的可靠性。