地表变形监测方法、终端设备及计算机可读存储介质

1.本发明公开了一种地表变形监测方法、终端设备及计算机可读存储介质，具体涉及基于sbas-insar相位恢复与改正技术的高精度地表变形监测方法、终端设备及计算机可读存储介质，属于地表形变监测技术领域。


背景技术：

2.基于sar影像的遥感技术能够有效应用于针对诸如地震、滑坡、火山、地裂缝等地质灾害的调查、监测及预警工作中。目前，sar技术主要有两大类，包括基于强度信息的偏移量跟踪和基于相位的sar干涉测量(insar)。
3.基于强度信息的sar偏移量跟踪技术利用两个sar影像之间强度信息，基于互相关技术估计距离向和方位向的偏移量，可以监测大梯度形变，但是测量精度在米级以上，很难精确监测小梯度形变。基于相位的sar干涉测量方法利用两部sar天线(或一部天线重复观测)来获取同一地区具有一定视角差的两幅单数复数影像，对两幅影像干涉得到的相位进行dem相位去除、滤波、相位解缠、大气相位去除等步骤后，实现对地表微小形变的观测。
4.sbas-insar技术是在基于相位的sar干涉测量方法的基础上，根据多景sar影像的成像时间和轨道矢量，生成符合一定时间和空间基线阈值的小基线集(sbas)，利用小基线集提供的多余观测进行平差来获取地表毫米级精度的形变时间序列。sbas-insar技术中的关键步骤为相位解缠，其目的是将sar影像干涉后缠绕在[-ππ]之间的相位根据itoh假设恢复为相对于参考点的真实相位。但是，在相位解缠过程中，受到地表植被以及其它失相干因素的影响，解缠相位中可能存在误差，当失相干现象严重时，影响insar监测结果精度。


技术实现要素：

[0005]
本技术的目的在于，提供一种地表变形监测方法、终端设备及计算机可读存储介质，以解决现有技术中存在的失相干现象严重，导致变形监测精度低的技术问题。
[0006]
本发明的第一方面提供了一种地表变形监测方法，包括：
[0007]
获取待监测区域不同监测时间的sar影像，得到多幅sar影像；
[0008]
根据所述多幅sar影像构建小基线集，并对所述小基线集中的多幅sar影像差分干涉后进行相位解缠，得到多幅相位解缠图；
[0009]
判断所述多幅相位解缠图的解缠相位是否满足预设条件，如否，则更新对应的相位解缠图；
[0010]
根据更新后的相位解缠图，确定所述待监测区域的变形序列。
[0011]
优选地，判断所述多幅相位解缠图的解缠相位是否满足预设条件，具体包括：
[0012]
从所述小基线集中的多幅sar影像中，确定基线满足预设要求的多组图像组，每组所述图像组中包括两幅sar影像；
[0013]
判断每组中所述两幅sar影像对应的相位解缠图的解缠相位是否为非空值。
[0014]
优选地，更新对应的相位解缠图，具体包括：
[0015]
以所述两幅sar影像的中间监测时间为中心，获取距所述中心设定时长内的多幅所述sar影像对应的相位解缠图；
[0016]
根据距所述中心设定时长内的多幅所述sar影像对应的相位解缠图，确定平均相位变化速率；
[0017]
利用所述平均相位变化速率，补全空值对应的所述解缠相位，更新对应的相位解缠图。
[0018]
优选地，所述平均相位变化速率根据第一公式确定，所述第一公式为：
[0019][0020]
式中，v为所述平均相位变化速率，ph
l
为第l幅相位解缠图中的解缠相位，δt
l
为第l幅相位解缠图对应的两幅sar影像的监测时间间隔，n为距所述中心设定时长内的多幅所述sar影像对应的相位解缠图的总数。
[0021]
优选地，利用所述平均相位变化速率，补全空值对应的所述解缠相位，具体包括：
[0022]
将所述平均相位变化速率和所述两幅sar影像的监测时间间隔相乘；
[0023]
利用所述相乘后的结果补全空值对应的所述解缠相位。
[0024]
优选地，更新对应的相位解缠图之后，还包括：
[0025]
根据所述小基线集中的多幅sar影像构建多个相位闭合环，确定存在解缠误差的相位；
[0026]
将所述解缠误差从对应的相位中去除；
[0027]
相应的，根据更新后的相位解缠图，确定所述待监测区域的变形序列，具体为：
[0028]
根据去除解缠误差后的多幅相位解缠图，确定所述待监测区域的变形序列。
[0029]
优选地，根据所述小基线集中的多幅sar影像构建多个相位闭合环，确定存在解缠误差的相位，具体包括：
[0030]
以所述小基线集中的多幅sar影像中的三幅构建一个相位闭合环，得到多个相位闭合环；
[0031]
获取每个相位闭合环的残差，结合预设的残差阈值，确定对应相位解缠图中存在解缠误差的相位。
[0032]
优选地，获取每个相位闭合环的残差，具体包括：
[0033]
根据第二公式确定每个相位闭合环的残差，所述第二公式为：
[0034][0035]
式中，为由第i、j、k监测时间获取的sar影像构建的相位闭合环的残差，为第i监测时间获取的sar影像与第j监测时间获取的sar影像之间的解缠相位，为第j监测时间获取的sar影像与第k监测时间获取的sar影像之间的解缠相位，为第k监测时间获取的sar影像与第i监测时间获取的sar影像之间的解缠相位。
[0036]
优选地，将所述解缠误差从对应的相位中去除，具体包括：
[0037]
将所述解缠误差视为粗差，确定所述粗差的估值；
[0038]
将所述粗差的估值约束到预设整数倍；
[0039]
将约束后的粗差的估值从对应的相位中去除。
[0040]
本发明的第二方面提供了一种终端设备，包括存储器、处理器以及存储在所述存储器中并可在所述处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现上述方法的步骤。
[0041]
本发明的第三方面提供了一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，所述计算机程序被处理器执行时实现上述方法的步骤。
[0042]
本发明的地表变形监测方法、终端设备及计算机可读存储介质，相较于现有技术，具有如下有益效果：
[0043]
本发明的地表变形监测方法及终端设备，能够充分利用多景sar影像通过sbas(small baseline subset)策略生成的冗余相位解缠图，自动识别并改正存在解缠误差与空值的解缠相位，提供精确的地表变形时间序列与完整的地表形变场，克服时间与空间失相干现象对sbas-insar监测结果的影响，特别是长时间失相干现象导致平差网络断开的影响。
附图说明
[0044]
图1为本发明地表变形监测方法的流程示意图；
[0045]
图2为本发明方法与现有方法的模拟数据对比图；
[0046]
图3为本发明实例提供的滑坡地表形变高精度反演小基线集时空基线分布图；
[0047]
图4为本发明实例提供的解缠相位恢复与改正sentinel-1绵沙湾滑坡2019年1月4日至2022年3月25日部分地表形变图；
[0048]
图5为本发明实例提供的解缠相位恢复与改正sentinel-1绵沙湾滑坡2019年1月4日至2022年3月25日地表形变时间序列图。
具体实施方式
[0049]
以下描述中，为了说明而不是为了限定，提出了诸如特定系统结构、技术之类的具体细节，以便透彻理解本发明实施例。然而，本领域的技术人员应当清楚，在没有这些具体细节的其它实施例中也可以实现本发明。在其它情况中，省略对众所周知的系统、装置、电路以及方法的详细说明，以免不必要的细节妨碍本发明的描述。
[0050]
如图1所示，本发明实施例的地表变形监测方法包括：
[0051]
步骤1、获取待监测区域不同监测时间的sar影像，得到多幅sar影像。
[0052]
本发明实施例中，为便于对多幅sar影像进行配准，在获取待监测区域的多幅sar影像的同时，还获取了待监测区域的数字高程模型(dig ital elevationmodel)数据，简称为dem数据。
[0053]
然后将所有的sar影像处理为单视复数图像。
[0054]
本发明实施例中，对于来自同一轨道的sar影像，考虑时间基线、空间基线和多普勒中心频率的变化，确定该轨道数据集中的主影像，并将数据集中所有sar影像共配准至主影像，得到配准后的多幅sar影像。
[0055]
后续步骤均基于配准后的sar影像进行操作。
[0056]
步骤2、根据多幅sar影像构建小基线集，并对小基线集中的多幅sar影像差分干涉后进行相位解缠，得到多幅相位解缠图，具体包括：
[0057]
获取满足预设时间基线阈值和预设空间基线阈值的sar影像并构建小基线集；例如可以获取小于预设时间基线阈值和预设空间基线阈值的sar影像并构建小基线集。
[0058]
对所述小基线集中的多幅sar影像进行差分干涉、多视、dem相位去除、自适应滤波、相位解缠等，得到多幅相位解缠图，多幅相位解缠图中包括冗余的解缠相位。
[0059]
步骤3、判断多幅相位解缠图的解缠相位是否满足预设条件，如否，则更新对应的相位解缠图，具体为：
[0060]
步骤31、从小基线集中的多幅sar影像中，确定基线满足预设要求的多组图像组，每组图像组中包括两幅sar影像。
[0061]
本发明实施例中的预设要求可以为最短时间基线，或者小于预设基线阈值的时间基线，优选为最短时间基线。最短时间基线通常为12天。
[0062]
步骤32、判断每组中两幅sar影像对应的相位解缠图的解缠相位是否满足预设条件，如否，则更新对应的相位解缠图。
[0063]
本发明实施例中，判断每组中两幅sar影像对应的相位解缠图的解缠相位是否满足预设条件，具体为：
[0064]
判断每组中两幅sar影像对应的相位解缠图的解缠相位是否为非空值。本发明实施例中的空值表明该点相干性较差，解缠相位不可靠。故在相位解缠步骤前对此像素掩膜，不进行相位解缠，因此相位为空。
[0065]
进一步地，本发明实施例中更新对应的相位解缠图的方法为：
[0066]
(a)以空值对应的两幅sar影像的中间监测时间为中心，获取距中心设定时长内的多幅sar影像对应的相位解缠图。具体为：以空值对应的两幅sar影像的中间监测时间为中心，在时间上开一个前后各设定时长的窗口，例如可为30天、50天、60天等，然后获取前后各设定时长窗口内的多幅sar影像对应的相位解缠图。
[0067]
(b)根据距中心设定时长内的多幅sar影像对应的相位解缠图，确定平均相位变化速率。
[0068]
在一个具体的实施例中，确定平均相位变化速率，具体为使用干涉图堆叠技术获取平均相位变化速率：
[0069][0070]
式(1)中，v为平均相位变化速率，ph
l
为第l幅相位解缠图中的解缠相位，δt
l
为第l幅相位解缠图对应的两幅sar影像的监测时间间隔，n为距中心设定时长内的多幅sar影像对应的相位解缠图的总数。
[0071]
(c)根据平均相位变化速率，补全空值对应的解缠相位，更新对应的相位解缠图，具体为：
[0072]
将平均相位变化速率与空值对应的两幅sar影像的监测时间间隔相乘，利用相乘后的结果补全空值对应的解缠相位。
[0073]
遍历所有的相位解缠图，直至补全所有的最短时间解缠相位，得到更新后的相位解缠图。
[0074]
步骤4、根据更新后的相位解缠图，确定待监测区域的变形序列。
[0075]
拟准检定(quad)是一种传统的大地测量数据处理方法，其通过建立误差估计模型来有效识别和修正粗差。quad可以通过选取拟准干涉相位，计算解缠相位的真误差，之后根据真误差分层聚束的现象确定出包含解缠粗差的相位。最后将估计出的解缠误差约束到2π整数倍，进而对含有解缠误差的相位进行修正并对空相位进行恢复。但是quad拟准相位的选取从m+2(m为sar影像数量)开始迭代，直到真误差分层的位置与拟准相位的数量相等，时间复杂度很大。并且当某个解缠误差的数值很大时，真误差分层的位置会出现确定错误的现象，导致非稳相位误选为拟稳相位，降低解缠误差改正效果甚至产生额外的解缠误差。
[0076]
为进一步提升变形监测精度，在另一个实施例中，在步骤3之后，还包括：
[0077]
(a)根据小基线集中的多幅sar影像构建多个相位闭合环，确定存在解缠误差的相位，具体为：
[0078]
以小基线集中的多幅sar影像中的三幅构建一个相位闭合环，得到多个相位闭合环。
[0079]
示例性地，假设在i,j,k三个时间点各自获取了一幅sar影像，三幅sar影像构成了一个时间上的相位闭合环。为第i监测时间获取的sar影像与第j监测时间获取的sar影像之间的解缠相位，为第j监测时间获取的sar影像与第k监测时间获取的sar影像之间的解缠相位，为第k监测时间获取的sar影像与第i监测时间获取的sar影像之间的解缠相位。每一个解缠相位均由同样的多个分量构成。以为例，说明其构成分量。由等分量构成。其中表示地表在视线向的形变量；表示因传感器位置不精确引入的残余轨道误差；为用于地形改正的外部dem误差；为相位解缠图上的解缠误差；为相位解缠图上的大气误差；表示噪声。
[0080]
获取每个相位闭合环的残差，结合预设的残差阈值，从补全后的相位解缠图中确定存在解缠误差的相位。
[0081]
本发明实施例中，相位闭合环的残差如公式(2)：
[0082][0083]
式(2)中，为由第i、j、k监测时间获取的sar影像构建的相位闭合环的残差，为闭合环残差中的形变相位分量，为闭合环残差中的形变相位分量，为大气相位残差，为大气相位残差，为残余轨道误差，为残余轨道误差，为dem相位分量，为dem相位分量，为解缠误差分量，为解缠误差分量，为随机噪声，
实际上，在闭合环残差部分，大气相位残差分量与形变分量恒为零，dem相位分量和噪声分量表现出随机性，残余轨道误差分量表现为趋势项分布，且对闭合环残差的贡献通常小于1弧度。而相位解缠误差分量在空间上表现为2π整周数的跳变。故设定π为相位闭合环的残差阈值，大于该阈值残差对应的相位闭合环视为存在解缠误差。然后通过多个相位闭合环的交集确定存在解缠误差的相位，其余相位认为不存在解缠误差，作为初始拟稳相位。
[0084]
然后判断初始拟稳相位是否为空值，若为空，则视为非拟稳相位，其余作为最终拟稳相位。
[0085]
(b)将解缠误差从对应的相位中去除，得到去除解缠误差后的相位解缠图。
[0086]
由于解缠误差通常数值较大，且非随机分布，因此将解缠误差视为粗差处理，然后确定粗差的估值；
[0087]
假设最终确定了n个解缠相位中含有b个粗差，可得到b个n维单位向量ej＝(0,
…
,0,1,0,
…
,0)
t
，对应第j个观测相位有粗差，即第j个分量为1，其余为0，cb＝(e1,
…
,eb)，则含有粗差的观测方程可用如下模型表示：
[0088][0089]
式(3)中，a是n
×
m维系数矩阵，秩为m；为m维未知参数的估值；为粗差估值；l是n维观测相位向量；v是观测残差。利用最小二乘法，可从式(3)中求出粗差估值为：
[0090][0091]
式(4)中，为粗差的估值，cb为存在解缠误差的多个相位解缠图组成的矩阵，为cb的转置矩阵，p为观测权阵，此处为单位矩阵，r为平差因子，l是n维观测相位向量。
[0092]
解缠误差只会以2π的整周倍出现，将式(4)求出的粗差估计值约束到2π的整周倍，将约束后的粗差的估值从对应的相位中去除。
[0093]
相应地，步骤4具体为：
[0094]
根据去除解缠误差后的多幅相位解缠图，确定所述待监测区域的变形序列，具体为使用秩亏自由网平差确定待监测区域的累积形变场，并从中提取变形序列：
[0095]
判断经过上述步骤的恢复与改正后的解缠相位是否为空，若为空，则忽略该观测值，即忽略为空的解缠相位，其余解缠相位经式(5)运算后，得到地表时间变形序列。
[0096][0097]
式中，为形变序列，a
+
为系数矩阵的伪逆，l为不为空的处理后的解缠相位，为sar影像获取日期地表的形变量，unw1为相位解缠图的相位值，m为sar影像的数量，s为不为
空的解缠相位的数量。
[0098]
本发明的第二方面提供了一种终端设备，包括存储器、处理器以及存储在存储器中并可在处理器上运行的计算机程序，处理器执行计算机程序时实现上述方法的步骤。
[0099]
本发明的第三方面提供了一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，所述计算机程序被处理器执行时实现上述方法的步骤。
[0100]
为验证本发明方法及终端设备的有效性，下面将以具体的模拟实验和实施例详述本发明的方法。
[0101]
模拟数据对比实验
[0102]
首先模拟单点的地表形变时间序列，由趋势项形变和线性项形变组成(ts＝(-2t+cos(10t)*t)/4；t∈(0,100])，共100个时间节点，其中tsi表示第i个时间点的累积形变。以四个时间点作为时间基线阈值生成小基线集，对累积形变做差得到390个解缠相位，ph
i,i+k
代表第i+k(k≤4)时间点与第i时间点的解缠相位。对解缠相位添加随机噪声、粗差(2π整周倍)以及空值，之后进行直接平差、拟准检定恢复与改正解缠相位后平差。
[0103]
图2为本发明方法与现有方法的模拟数据对比图。其中ts1是模拟的单点时间形变序列，由线性项和趋势项两部分组成。对模拟的时间序列根据sbas策略生成冗余解缠相位，并随机添加相位解缠误差(包括空值)。ts2为不进行相位恢复与改正的冗余解缠相位平差得到的时间形变序列；ts3为通过拟准检定改正的冗余解缠相位平差得到的时间形变序列；ts4为通过本发明的方法恢复与改正的冗余解缠相位平差得到的时间形变序列。从图2中可以发现本发明提出方法得到的时间形变序列与模拟时间形变序列基本吻合，得到的时序曲线抖动性更小，精确度更高。
[0104]
实施例
[0105]
本实施例选择位于我国金沙江流域的绵沙湾滑坡。金沙江是长江的上游，流经中国西部的青海、西藏、四川和云南四省。白鹤滩水电站是金沙江下游开发的第二座梯级水电站。2021年4月22日，白鹤滩水电站开始蓄水，水位从开始的660.35米上升到了916.51米。库区水位的快速上升破坏了绵沙湾坡体的稳定性，诱发了滑坡的发生。
[0106]
实验选用的sar数据为2019年1月4日至2022年3月25日覆盖金沙江流域绵沙湾滑坡的94景sentinel-1数据。由于滑坡表面覆盖有茂密的植被，严重的信号失相干导致sar影像的解缠相位中出现空值与误差。利用本发明改进的sbas-insar相位恢复与改正技术对该滑坡进行高精度时序解算。
[0107]
首先，选取2020年10月1日获取的sar影像作为主影像，其余图像作为从影像，并将所有从影像精确配准至主影像。其次，基于精确配准的sar影像通过设置时空基线阈值方式进行干涉对的组合，随后对获取的干涉对经过相位干涉、多视、dem相位去除、自适应滤波、相位解缠等流程后得到相位解缠图。对冗余的解缠相位恢复与改正后，使用秩亏自由网平差反演滑坡表面的高精度时间序列。
[0108]
图3所示为本发明所选取的用于滑坡表面反演的干涉对时空基线分布图。设定的时间基线为40天，空间基线为150米，94景sar影像构成了262个干涉对。
[0109]
图4所示为本发明所提出方法计算获得的sentinel-1sar影像绵沙湾滑坡部分日期相较于起始日期(2019年1月4日)的地表累积形变图。传统及现有的sbas-insar方法由于解缠误差以及空值的影响，无法得到完整的形变场。从图4中可以看到，本发明所提方法准
确、完整的获得了受失相干影响严重的绵沙湾滑坡2019年1月4日至2022年3月25日的累积形变场，其最大累积形变达到25cm。
[0110]
图5所示为本发明所提出方法计算获得的sentinel-1sar影像绵沙湾滑坡2019年1月4日至2022年3月25日地表形变时间序列。本发明所提方法采用sentinel-1sar影像精确反演了绵沙湾滑坡地表形变时间演化特征，并成功捕获到该滑坡受库区蓄水影响地表形变加速过程，相较传统方法提升了地表形变反演精度。
[0111]
本发明的地表变形监测方法及终端设备，能够充分利用多景sar影像通过sbas(small baseline subset)策略生成的冗余相位解缠图，自动识别并改正存在解缠误差与空值的解缠相位，提供精确的地表变形时间序列与完整的地表形变场，克服时间与空间失相干现象对sbas-insar监测结果的影响，特别是长时间失相干现象导致平差网络断开的影响。在恢复和改正解缠相位后，根据卫星的航向角、入射角以及地形因子可以获得滑坡、地面沉降、火山、路基形变等的坡向形变，能够准确地得到形变特征。为滑坡、地面沉降、火山、路基形变等的高精度监测提供了一种改进的方法，此工作对于滑坡监测具有重要意义。
[0112]
以上所述，仅是本技术的几个实施例，并非对本技术做任何形式的限制，虽然本技术以较佳实施例揭示如上，然而并非用以限制本技术，任何熟悉本专业的技术人员，在不脱离本技术技术方案的范围内，利用上述揭示的技术内容做出些许的变动或修饰均等同于等效实施案例，均属于技术方案范围内。