CN117761716A - 融合CR和GNSS的InSAR滑坡形变时序监测方法和存储介质 - Google Patents
融合CR和GNSS的InSAR滑坡形变时序监测方法和存储介质 Download PDFInfo
- Publication number
- CN117761716A CN117761716A CN202311775084.1A CN202311775084A CN117761716A CN 117761716 A CN117761716 A CN 117761716A CN 202311775084 A CN202311775084 A CN 202311775084A CN 117761716 A CN117761716 A CN 117761716A
- Authority
- CN
- China
- Prior art keywords
- deformation
- radar
- gnss
- sight
- line
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
- 238000000034 method Methods 0.000 title claims abstract description 50
- 238000012544 monitoring process Methods 0.000 title claims abstract description 26
- 230000001427 coherent effect Effects 0.000 claims abstract description 30
- 238000004458 analytical method Methods 0.000 claims abstract description 17
- 238000004364 calculation method Methods 0.000 claims description 17
- 238000001914 filtration Methods 0.000 claims description 14
- 238000005259 measurement Methods 0.000 claims description 10
- 230000000007 visual effect Effects 0.000 claims description 10
- 238000012300 Sequence Analysis Methods 0.000 claims description 7
- 238000012545 processing Methods 0.000 claims description 7
- 238000001228 spectrum Methods 0.000 claims description 4
- 238000005305 interferometry Methods 0.000 claims description 3
- 238000012795 verification Methods 0.000 claims description 3
- 238000000605 extraction Methods 0.000 abstract description 5
- 238000010586 diagram Methods 0.000 description 7
- 230000002123 temporal effect Effects 0.000 description 4
- 230000002452 interceptive effect Effects 0.000 description 3
- 241001085205 Prenanthella exigua Species 0.000 description 2
- 238000006243 chemical reaction Methods 0.000 description 2
- 230000004927 fusion Effects 0.000 description 2
- 230000003287 optical effect Effects 0.000 description 2
- 238000012731 temporal analysis Methods 0.000 description 2
- 238000000700 time series analysis Methods 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
- 238000004804 winding Methods 0.000 description 1
Landscapes
- Radar Systems Or Details Thereof (AREA)
Abstract
本发明公开了一种融合CR和GNSS的InSAR滑坡形变时序监测方法和存储介质。根据本发明的实施例，通过将至少一套CR和GNSS设备布设在滑坡的同一位置处，并将GNSS设备记录的地理坐标转换为SAR距离多普勒坐标，从而识别CR在SAR主影像上的位置，然后将GNSS设备观测值转换为雷达视线向观测值并确定时间基线阈值，进而根据时间基线阈值筛选确定既能够满足相干性、又能恢复出正确形变信息的差分干涉图，基于多幅差分干涉图进行干涉点目标分析得到相干点目标的形变速率和形变时间序列。因此，采用本发明能够将CR和GNSS定位、GNSS点上形变时序观测以及高分辨率InSAR提取面上形变信息进行无缝衔接，实现大变形滑坡形变演化过程分析。
Description
技术领域
本发明涉及地质灾害监测技术领域，更为具体而言，涉及一种融合CR和GNSS的InSAR滑坡形变时序监测方法和存储介质。
背景技术
InSAR形变监测的原理为利用变化前后的两幅星载SAR影像差分干涉处理得到缠绕相位，经相位解缠后，得到绝对相位来表征地表高程变化信息。如果小范围内出现较大变形时，往往会造成干涉图中相位急剧变化，使得干涉条纹变得模糊，从而无法准确的恢复形变相位。换言之，现有技术中缺少对大变形滑坡形变演化过程进行监测分析的方法。
发明内容
本发明实施方式提供了一种融合CR和GNSS的InSAR滑坡形变时序监测方法和存储介质，将CR和GNSS定位、GNSS点上形变时序观测以及高分辨率InSAR提取面上形变信息进行无缝衔接，提升了所能探测的最大形变梯度能力，测量的滑坡最大形变量级由原来的几十cm提高至几m，可更为真实的反映滑坡形变特征，实现了大变形滑坡形变演化过程分析。
根据本发明的第一方面，本发明实施方式提供了一种融合CR和GNSS的InSAR滑坡形变时序监测方法，该方法包括：获取覆盖滑坡区的卫星雷达数据，并在所述卫星雷达数据中选取SAR主影像；获取从滑坡形变显著区域布设的至少一套CR和GNSS设备发送的数据，其中，至少一套CR和GNSS设备布设在滑坡的同一位置用于精度验证；将所述GNSS设备记录的地理坐标转换为所述SAR主影像的距离多普勒坐标，并识别CR在所述SAR主影像上的位置；根据所述SAR主影像的参数文件计算覆盖所述滑坡区的卫星轨道信息，所述卫星轨道信息包括卫星轨道方位角和雷达入射角；根据所述卫星轨道信息将GNSS设备观测值转换为雷达视线向观测值；基于所述雷达视线向观测值确定雷达视线向形变时间序列，根据所述雷达视线向形变时间序列计算理论雷达卫星拍摄时间间隔，并设置时间基线阈值；将所述卫星雷达数据中的所有影像进行亚像元配准到所述主影像的距离多普勒坐标系下，得到配准后的SAR影像；在所述配准后的SAR影像中选取时间间隔不超过所述时间基线阈值的影像作为干涉像对，并对所述干涉像对进行差分干涉处理得到差分干涉图集；根据所述差分干涉图集计算初始视线向形变速率；对于所述配准后的SAR影像中时间间隔超过所述时间基线阈值的相邻影像，根据所述初始视线向形变速率和所述相邻影像的时间间隔确定所述相邻影像的替代差分干涉图；将所述差分干涉图集和替代差分干涉图作为干涉点目标分析方法IPTA的标准数据输入，并选取包括CR点的相干点目标，对所述差分干涉图集和替代差分干涉图进行形变时序分析，得到所述相干点目标的形变速率和形变时间序列。
本发明上述实施方式通过将至少一套CR和GNSS设备布设在滑坡的同一位置处，并将GNSS设备记录的地理坐标转换为SAR距离多普勒坐标，可以识别CR在SAR主影像上的位置，然后将GNSS设备观测值转换为雷达视线向观测值并确定时间基线阈值，进而根据时间基线阈值确定既能够满足相干性，又能恢复出正确形变信息的差分干涉图，基于多幅差分干涉图进行干涉点目标分析得到相干点目标的形变速率和形变时间序列。由此，将CR和GNSS定位、GNSS点上形变时序观测以及高分辨率InSAR提取面上形变信息进行无缝衔接，可真实的反映滑坡形变特征，实现了大变形滑坡形变演化过程分析。
在本发明的一些实施方式中，所述卫星雷达数据是覆盖同一滑坡区的不同模式数据，并且，所述不同模式数据的入射角基本一致。
在本发明的一些实施方式中，根据距离方程R＝|Rs-Rt|，以及多普勒方程将所述GNSS设备记录的地理坐标转换为所述SAR主影像的距离多普勒坐标；其中，R为传感器到点目标的距离，Rs为传感器的位置，vs为传感器的速度矢量，Rt为点目标的位置，vt为点目标的速度矢量，λ为雷达波长，fd为多普勒频率。
在本发明的一些实施方式中，根据下述计算公式将GNSS设备观测值转换为雷达视线向观测值：
[Unsinα-Uecosα]sinθ+Uucosθ+δlos＝glos
其中，Ue、Un、Uu分别为GNSS设备东西、南北、垂向形变分量，α为卫星轨道方位角，θ为雷达入射角，glos为GNSS设备观测值转换得到的雷达视线向观测值。
在本发明的一些实施方式中，根据所述雷达视线向形变时间序列计算理论雷达卫星拍摄时间间隔包括：
根据下述公式计算雷达视线向形变速率：
其中，t0为起始观测时间，ti为截止观测时间，Δvglos为起始观测时间t0到截止观测时间ti之间的雷达视线向形变速率；
将Δvglos代入下述计算公式中的Δv计算理论雷达卫星拍摄时间间隔：
其中，Δv为同一像元内雷达视线向形变速率变化值，λ为雷达波长，Δt为理论雷达卫星拍摄时间间隔。
在本发明的一些实施方式中，根据下述形变速率计算公式计算初始视线向形变速率：
其中，Δtk为所述差分干涉图集中的第k组干涉像对的时间间隔，φk为所述差分干涉图集中的第k组干涉像对的解缠相位，N为所述干涉像对的总数值，p为所述初始视线向形变速率。
在本发明的一些实施方式中，根据所述初始视线向形变速率和所述相邻影像的时间间隔确定所述相邻影像的替代差分干涉图包括：所述初始视线向形变速率和所述相邻影像的时间间隔的乘积为相邻影像的替代视线向形变量；根据所述替代视线向形变量确定所述相邻影像的替代差分干涉图。
在本发明的一些实施方式中，所述方法还包括：采用Goldstein-Werner滤波方法对所述差分干涉图集和替代差分干涉图进行滤波；利用最小费用流算法对滤波后的差分干涉图进行空间相位解缠；采用幅度离散指数和相干系数阈值综合选取SAR影像频谱特性稳定的所述相干点目标；根据下述公式计算得到所有相干点目标的形变速率和形变时间序列：
其中，为第k个干涉像对中相干点目标i和j之间的相位差，θ为入射角，/>为垂直基线，tk为时间间隔，/>为邻点目标i和j间的相对高程差，/>为邻点目标i和j间的相对形变速率，/>为其它相位，λ为雷达波长，R为传感器到点目标的距离；针对所述相干点目标的形变速率和形变时间序列，采用干涉点目标分析方法IPTA进行形变时序分析。
在本发明的一些实施方式中，所述方法还包括：根据下述公式计算同一时间段内CR点上InSAR测量值与GNSS转换为视线向形变观测值的互差的均方根误差σ：
其中，xi为CR点上InSAR测量值，yi为CR点上GNSS转换为视线向形变观测值，n为统计样本个数；利用所述CR点上InSAR测量值与GNSS转换为视线向形变观测值的互差的均方根误差评判InSAR结果的测量精度。
根据本发明的第二方面，本发明实施方式提供一种计算机可读存储介质，其上存储有计算机可读指令，所述计算机可读指令被处理器执行时，使得计算机执行如下操作：所述操作包括如上任意一种实施方式所述方法所包含的步骤。
综上所述，本发明实施方式提供的一种融合CR和GNSS的InSAR滑坡形变时序监测方法和存储介质，将CR和GNSS定位、GNSS点上形变时序观测以及高分辨率InSAR提取面上形变信息进行无缝衔接，能够实现大变形滑坡形变演化过程分析。
附图说明
图1是根据本发明一种实施方式的融合CR和GNSS的InSAR滑坡形变时序监测方法的流程示意图；
图2是采用图1所示融合CR和GNSS的InSAR滑坡形变时序监测方法的具体示例的流程示意图；
图3是在雷达影像上识别的CR点图；
图4是24景SAR影像形成的35个干涉像对的空间和时间基线分布图，横轴为时间(年/月/日)，纵轴为垂直基线，Bperp min和max为垂直基线最小和最大值，Delta_T min和max为时间基线的最小和最大值；
图5a是2023/04/04-2023/04/08(4天间隔)滤波后差分干涉图，-π到π是指一个2π周期的相位变化；
图5b是2023/04/04-2023/04/12(8天间隔)滤波后差分干涉图，-π到π是指一个2π周期的相位变化；
图5c是2023/04/04-2023/04/16(12天间隔)滤波后差分干涉图，-π到π是指一个2π周期的相位变化；；
图6是LT-1卫星多模式InSAR交互差分干涉图(2023/06/19_2023/06/27)，即3m(2023/06/19)和6m(2023/06/27)分辨率SAR影像进行差分干涉处理；
图7是泄流坡滑坡InSAR形变速率图；
图8a是选取的泄流坡滑坡后缘典型相干点目标InSAR形变时间序列图，InSAR_ori和InSAR_cor分别为超过8天的相邻影像的原始差分干涉图和替代差分干涉图参与计算得到的形变时间序列；
图8b是选取的泄流坡滑坡中部典型相干点目标InSAR形变时间序列图，InSAR_ori和InSAR_cor分别为超过8天的相邻影像的原始差分干涉图和替代差分干涉图参与计算得到的形变时间序列；
图8c是选取的泄流坡滑坡前缘典型相干点目标InSAR形变时间序列图，InSAR_ori和InSAR_cor分别为超过8天的相邻影像的原始差分干涉图和替代差分干涉图参与计算得到的形变时间序列；
图9是CR-2点上InSAR测量值与GNSS转换为视线向形变观测值对比形变时序图。
具体实施方式
以下结合附图和具体实施方式对本发明的各个方面进行详细阐述。本领域技术人员应当理解，下述的各种实施方式只用于举例说明，而非用于限制本发明的保护范围，在不脱离本发明实质的情况下，还可以对本发明的各种实施方式进行各种不同的组合。
下面对本文中使用的术语进行简要说明。
CR：Corner Reflector，角反射器。
GNSS：Global Navigation Satellite System，全球导航卫星系统。
InSAR：Interferometric Synthetic Aperture Radar，合成孔径雷达干涉测量，指采用合成孔径雷达的干涉测量技术。
SAR：Synthetic Aperture Radar，合成孔径雷达。
IPTA：Interferometric Point Target Analysis，干涉点目标分析方法。
Perpendicular Baseline：垂直基线。
SLC：Single Look Complex，单视复数影像，属于Level-1影像。
图1是根据本发明一种实施方式的融合CR和GNSS的InSAR滑坡形变时序监测方法的流程示意图。
如图1所示，在本发明的一种实施方式中，所述滑坡形变时序监测方法可包括：步骤S101、步骤S102、步骤S103、步骤S104、步骤S105、步骤S106、步骤S107、步骤S108、步骤S109、步骤S110和步骤S111，下面对上述步骤进行具体的描述。
在步骤S101中，获取覆盖滑坡区的卫星雷达数据，并在所述卫星雷达数据中选取SAR主影像。在一些实施方式中，所述卫星雷达数据是覆盖同一滑坡区的不同模式数据，并且，所述不同模式数据的入射角基本一致。
在本实施例中，所述卫星雷达数据为LT-1卫星雷达数据，进一步的，选取LT-1卫星雷达数据中模式1数据作为SAR主影像，首次实现了利用LT-1卫星星座不同模式SAR影像的交互干涉进行大变形滑坡形变监测。其中，陆地探测-1号(LT-1)01组卫星由2颗完全相同的L波段全极化雷达卫星组成。卫星设计有6种观测模式，以6m/100km(模式2)和3m/50km(模式1)作为形变监测主要模式，理论严格回归周期单星8天、双星4天，实际的重访周期取决于拍摄规划以及选用的模式，具备单/双星重轨干涉以及双星单航过干涉能力。并且，LT-1卫星是迄今为止重访周期最短的L波度雷达卫星星座，适用于实现对大变形滑坡的监测。
在步骤S102中，获取从滑坡形变显著区域布设的至少一套CR和GNSS设备发送的数据，其中，至少一套CR和GNSS设备布设在滑坡的同一位置用于精度验证。在一些实施方式中，GNSS设备发送的数据包括GNSS设备记录的地理坐标以及东西Ue、南北Un、垂向Uu形变分量；根据先验知识在滑坡形变显著(或最大)区域布设CR和GNSS设备。
在步骤S103中，将所述GNSS设备记录的地理坐标转换为所述SAR主影像的距离多普勒坐标，并识别CR在所述SAR主影像上的位置。
在一些实施方式中，根据下述公式(1)将所述GNSS设备记录的地理坐标转换为所述SAR主影像的距离多普勒坐标：
距离方程：R＝|Rs-Rt|
多普勒方程：
其中，R为传感器到点目标的距离，Rs为传感器的位置，vs为传感器的速度矢量，Rt为点目标的位置，vt为点目标的速度矢量，λ为雷达波长，fd为多普勒频率。
至少一套CR和GNSS设备布设在滑坡的同一位置处，通过将GNSS设备记录的地理坐标转换为SAR距离多普勒坐标，初步识别CR在SAR主影像上的位置，然后在SAR主影像上初步识别出CR行列号的位置附近，通过目视解译确认SAR主影像上亮白方形点为CR点。
在步骤S104中，根据所述SAR主影像的参数文件计算覆盖所述滑坡区的卫星轨道信息，所述卫星轨道信息包括卫星轨道方位角和雷达入射角。
在步骤S105中，根据所述卫星轨道信息将GNSS设备观测值转换为雷达视线向观测值。在一些实施方式中，根据下述计算公式(2)将GNSS设备观测值转换为雷达视线向观测值：
[Unsinα-Uecosα]sinθ+Uucosθ+δlos＝glos (2)
其中，Ue、Un、Uu分别为GNSS设备东西、南北、垂向形变分量，α为卫星轨道方位角，θ为雷达入射角，glos为GNSS设备观测值转换得到的雷达视线向观测值。
在步骤S106中，基于GNSS雷达视线向观测值确定雷达视线向形变时间序列，根据所述雷达视线向形变时间序列计算理论雷达卫星拍摄时间间隔，并设置时间基线阈值。
在一些实施方式中，根据所述雷达视线向形变时间序列计算理论雷达卫星拍摄时间间隔包括如下步骤：
首先，根据下述计算公式(3)计算雷达视线向形变速率：
其中，t0为起始观测时间，ti为截止观测时间，Δvglos为起始观测时间t0到截止观测时间ti之间的雷达视线向形变速率；
然后，将根据计算公式(3)得到的Δvglos代入下述计算公式(4)中的Δv计算理论雷达卫星拍摄时间间隔：
其中，Δv为同一像元内雷达视线向形变速率变化值，λ为雷达波长，Δt为理论雷达卫星拍摄时间间隔。根据上述公式(4)可知，雷达卫星的重访周期越短，越有利于监测大变形现象。
进一步的，根据所述理论雷达卫星拍摄时间间隔结合SAR数据获取情况设置时间基线阈值为m天。
在步骤S107中，将所述卫星雷达数据中的所有影像进行亚像元配准到所述主影像的距离多普勒坐标系下，得到配准后的SAR影像。
在步骤S108中，在所述配准后的SAR影像中选取时间间隔不超过所述时间基线阈值的影像作为干涉像对，并对所述干涉像对进行差分干涉处理得到差分干涉图集。结合公式(4)得到的理论雷达卫星拍摄时间间隔确定的时间基线阈值确定既能够满足相干性，又能恢复出正确形变信息的差分干涉图，进而基于差分干涉图集中的数据进行初始形变速率估计。
在步骤S109中，根据所述差分干涉图集计算初始视线向形变速率。在一些实施方式中，根据下述形变速率计算公式(5)计算初始视线向形变速率：
其中，Δtk为所述干涉像对集中的第k组干涉像对的时间间隔，φk为所述干涉像对集中的第k组干涉像对的解缠相位，N为所述干涉像对的总数值，p为所述初始视线向形变速率。
在步骤S110中，对于所述配准后的SAR影像中时间间隔超过所述时间基线阈值的相邻影像，根据所述初始视线向形变速率和所述相邻影像的时间间隔确定所述相邻影像的替代差分干涉图。在一些实施方式中，根据所述初始视线向形变速率和所述相邻影像的时间间隔确定所述相邻影像的替代差分干涉图包括：所述初始视线向形变速率和所述相邻影像的时间间隔的乘积为相邻影像的替代视线向形变量；根据所述替代视线向形变量确定所述相邻影像的替代差分干涉图。
在步骤S111中，将所述差分干涉图集和替代差分干涉图作为干涉点目标分析方法IPTA的标准数据输入，并选取包括CR点的相干点目标，对所述差分干涉图集和替代差分干涉图进行形变时序分析，得到所述相干点目标的形变速率和形变时间序列。
在一些实施方式中，根据下述步骤获取相干点目标的形变速率和形变时间序列：
步骤(1)，采用Goldstein-Werner滤波方法对所述差分干涉图集和替代差分干涉图进行滤波；
步骤(2)，利用最小费用流算法对滤波后的差分干涉图进行空间相位解缠；
步骤(3)，采用幅度离散指数和相干系数阈值综合选取SAR影像频谱特性稳定的所述相干点目标；
步骤(4)，根据下述公式(6)计算得到所有相干点目标的形变速率和形变时间序列：
其中，为第k个干涉像对中相干点目标i和j之间的相位差，θ为雷达入射角，为垂直基线，tk为时间间隔，/>为邻点目标i和j间的相对高程差，/>为邻点目标i和j间的相对形变速率，λ为雷达波长，R为传感器到点目标的距离，/>为其它相位(残余相位)，其他相位例如包括轨道误差相位、大气延迟相位、噪声相位等；
步骤(5)，针对所述相干点目标的形变速率和形变时间序列，采用干涉点目标分析方法IPTA进行形变时序分析。
在本实施方式中，步骤S101和S102之间没有固定顺序。
采用本发明实施方式的上述方法，通过将至少一套CR和GNSS设备布设在滑坡的同一位置处，并将GNSS设备记录的地理坐标转换为SAR距离多普勒坐标，可以识别CR在SAR主影像上的位置，然后将GNSS设备观测值转换为雷达视线向观测值并确定时间基线阈值，进而根据时间基线阈值确定既能够满足相干性，又能避免相位模糊以恢复出正确形变信息的差分干涉图，基于多幅差分干涉图进行干涉点目标分析得到相干点目标的形变速率和形变时间序列。由此，将CR和GNSS定位、GNSS点上形变时序观测以及高分辨率InSAR提取面上形变信息进行无缝衔接，提升了所能探测的最大形变梯度能力，测量的滑坡最大形变量级由原来的几十cm提高至几m，可更为真实的反映滑坡形变特征，实现大变形滑坡形变演化过程分析。
在进一步的实施方式中，所述滑坡形变时序监测方法还包括：根据下述公式(7)计算同一时间段内CR点上InSAR测量值与GNSS转换为视线向形变观测值的互差的均方根误差σ：
其中，xi为CR点上InSAR测量值，yi为CR点上GNSS转换为视线向形变观测值，n为统计样本个数；利用所述CR点上InSAR测量值与GNSS转换为视线向形变观测值的互差的均方根误差评判InSAR结果的测量精度。
本发明给出一种采用图1所示融合CR和GNSS的InSAR滑坡形变时序监测方法实现大变形滑坡形变时序监测的具体示例，该示例以甘肃舟曲泄流坡滑坡监测为例，采用LT-1卫星模式1和模式2数据以及布设在滑坡体上的CR和GNSS设备联合观测，下面结合图2至图9进行具体说明：
图2是采用图1所示融合CR和GNSS的InSAR滑坡形变时序监测方法的具体示例的流程示意图；图3是在雷达影像上识别的CR点图；图4是24景SAR影像形成的35个干涉像对的空间和时间基线分布图，横轴为时间(年/月/日)，纵轴为垂直基线，Bperp min和max为垂直基线最小和最大值，Delta_T min和max为时间基线的最小和最大值；图5a是2023/04/04-2023/04/08(4天间隔)滤波后差分干涉图；图5b是2023/04/04-2023/04/12(8天间隔)滤波后差分干涉图；图5c是2023/04/04-2023/04/16(12天间隔)滤波后差分干涉图；图6是LT-1卫星多模式InSAR交互差分干涉图(2023/06/19_2023/06/27)，即3m(2023/06/19)和6m(2023/06/27)分辨率SAR影像进行差分干涉处理；图7是泄流坡滑坡InSAR形变速率图；图8a是选取的泄流坡滑坡后缘典型相干点目标InSAR形变时间序列图，InSAR_ori和InSAR_cor分别为超过8天的相邻影像的原始差分干涉图和替代差分干涉图参与计算得到的形变时间序列；图8b是选取的泄流坡滑坡中部典型相干点目标InSAR形变时间序列图，InSAR_ori和InSAR_cor分别为超过8天的相邻影像的原始差分干涉图和替代差分干涉图参与计算得到的形变时间序列；图8c是选取的泄流坡滑坡前缘典型相干点目标InSAR形变时间序列图，InSAR_ori和InSAR_cor分别为超过8天的相邻影像的原始差分干涉图和替代差分干涉图参与计算得到的形变时间序列；图9是InSAR结果与GNSS观测值对比形变时序图。
如图2所示，所述方法包括：
步骤一：LT-1卫星单视复数影像选取
选取覆盖泄流坡滑坡的LT-1卫星雷达数据，其中模式1数据6景(2023年6月27日-2023年8月18日)，模式2数据18景(2023年3月15日-2023年6月19日)，雷达数据主要参数如下表1所示。
表1
| 雷达参数 | 模式1 | 模式2 |
| 载波频率 | 1.26e+009Hz | 1.26e+009Hz |
| 中心入射角 | 31.2851° | 31.2911° |
| 轨道方位角 | 192.477621° | 192.4766964° |
| 距离向像元大小 | 1.665514m | 1.665514m |
| 方位向像元大小 | 1.281204m | 3.046507m |
步骤二：雷达影像CR识别
根据现场及以往掌握的滑坡滑动情况，在滑坡体中下部同一位置布设2套CR(图7中CR-1和CR-2)和GNSS设备。选取LT-1卫星雷达数据中的2023/06/27影像(模式1)作为SAR主影像，将GNSS设备记录的地理坐标按照公式(1)转换为2023/06/27影像(SAR主影像)的距离多普勒坐标框架下，初步识别出CR在2023/06/27影像上的行列号，最终通过目视解译确认图3中的亮白方形点为CR点。
步骤三：大变形滑坡多分辨率InSAR形变时序分析
(1)利用2023/06/27影像参数文件，计算获取覆盖泄流坡滑坡区的卫星轨道方位角为192.4786°，入射角为31.2672°。根据公式(2)将CR-1处的GNSS观测值转换为雷达视线向形变时间序列，并得到GNSS雷达视线向形变速率约为1.1cm/d，根据公式(3)和(4)计算得到泄流坡滑坡中、下部监测约需10天LT-1雷达卫星重访周期。将所有的雷达影像配准到2023/06/27影像的距离多普勒坐标下(配准精度优于0.1个像元)，并且，由于滑坡后缘以及中部的形变速率超过CR-1位置的量级，设置时间间隔阈值(即步骤S106中的时间基线阈值)为8天，在配准后的SAR影像中选取时间间隔不超过8天的影像作为干涉像对，得到33个干涉像对，然后根据公式(5)计算出初始视线向形变速率(即前述步骤S109中的初始视线向形变速率)。
其中，根据CR-1处的GNSS观测数据和公式(3)计算得到测量滑坡中、下部约需10天重访周期，但滑坡中后部形变量更大，所以需要更快的重访周期。结合LT-1卫星4天或8天的重访周期，同时考虑到让更多的数据参与计算，设置时间间隔为8天。
(2)LT-1卫星自2023年6月22日起，由在轨测试模式(模式2)转入模式1开展中国陆域常规拍摄，同一轨道的实际重访周期会出现超过8天的情况，如2023/07/09与2023/08/10相邻的2期时间间隔为32天，导致解缠相位不连续，因此，InSAR形变时间序列分析时，根据公式(5)计算得到的初始视线向形变速率乘以32天代替2023/07/09_2023/08/10的原始差分干涉图。由此，根据时间间隔不超过8天的干涉像对集以及时间间隔超过8天的相邻的影像对构建如图4所示的35个干涉像对，以形成连续的形变时间序列。其中，用于做时间序列分析的差分干涉图集包括2部分：一部分是时间间隔不超过8天的差分干涉图(有33个)；另一部分是卫星拍摄过程中，会出现重访周期超过8天的情况，为了能让形变分析时间连续，时间间隔超过8天的相邻影像对也参与计算(有2个)，相邻影像的差分干涉图是利用初始形变速率乘以相邻影像的时间间隔计算得到。
然后设置相干系数阈值对35个干涉像对进行差分干涉处理，距离向和方位向的多视比为2:4，利用机载LiDAR生产的2m分辨率DEM模拟干涉相位补偿地形相位后得到差分干涉图；采用Goldstein-Werner滤波方法对差分干涉图进行滤波，得到滤波后的差分干涉图(例如图5a至5c所示)；利用最小费用流算法(minimum cost flow，MCF)对滤波后的差分干涉图进行空间相位解缠；采用幅度离散指数和相干系数阈值综合选取SAR影像频谱特性稳定的相干点目标，并保证相干点目标中有CR点；然后利用公式(6)计算得到所有相干点目标的形变速率和形变时间序列，得到如图7和图8a至8c所示的泄流坡滑坡InSAR形变速率图和形变时间序列图。
步骤四：CR-InSAR与GNSS数据互检校
自2023年3月15日至8月18日共收集24期数据，因此CR点的形变时间序列共有24个时间节点，从CR-2处的GNSS测量值转换为视线向形变观测序列中找到与雷达拍摄日期相对应的观测数据，求取CR-InSAR(即CR点上InSAR测量值)与GNSS转换为视线向形变观测值之间的互差，并利用公式(7)分别计算其均方根误差。如图9所示，计算的均方根误差为7mm，CR-InSAR与GNSS观测序列趋势非常一致，达到了预期监测效果。
通过以上实施方式的描述，本领域的技术人员可以清楚地了解到本发明可借助软件结合硬件平台的方式来实现。基于这样的理解，本发明的技术方案对背景技术做出贡献的全部或者部分可以以软件产品的形式体现出来，该计算机软件产品可以存储在存储介质中，如ROM/RAM、磁碟、光盘等，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行本发明各个实施方式或者实施方式的某些部分所述的方法。
对应的，本发明实施方式还提供一种计算机可读存储介质，其上存储有计算机可读指令或程序，所述计算机可读指令或程序被处理器执行时，使得处理器执行如下操作：所述操作包括如上任意一种实施方式所述方法所包含的步骤，在此不再赘述。其中，所述存储介质可以包括：例如，光盘、硬盘、软盘、闪存、磁带等。本发明上述实施方式的部分或全部可以通过计算机设备实现，所述计算机设备包括上述存储介质和处理器，所述处理器执行所述存储介质中的计算机可读指令以执行上述任一实施方式所记载的步骤、操作、处理。
最后应说明的是：以上实施方式仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施方式对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施方式所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施方式技术方案的精神和范围。因此本发明的保护范围应以权利要求为准。
Claims (10)
1.一种融合角反射器CR和全球导航卫星系统GNSS的合成孔径雷达干涉测量InSAR滑坡形变时序监测方法，所述方法包括：
获取覆盖滑坡区的卫星雷达数据，并在所述卫星雷达数据中选取SAR主影像；
获取从滑坡形变显著区域布设的至少一套CR和GNSS设备发送的数据，其中，至少一套CR和GNSS设备布设在滑坡的同一位置用于精度验证；
将所述GNSS设备记录的地理坐标转换为所述SAR主影像的距离多普勒坐标，并识别CR在所述SAR主影像上的位置；
根据所述SAR主影像的参数文件计算覆盖所述滑坡区的卫星轨道信息，所述卫星轨道信息包括卫星轨道方位角和雷达入射角；
根据所述卫星轨道信息将GNSS设备观测值转换为雷达视线向观测值；
基于所述雷达视线向观测值确定雷达视线向形变时间序列，根据所述雷达视线向形变时间序列计算理论雷达卫星拍摄时间间隔，并设置时间基线阈值；
将所述卫星雷达数据中的所有影像进行亚像元配准到所述主影像的距离多普勒坐标系下，得到配准后的SAR影像；
在所述配准后的SAR影像中选取时间间隔不超过所述时间基线阈值的影像作为干涉像对，并对所述干涉像对进行差分干涉处理得到差分干涉图集；
根据所述差分干涉图集计算初始视线向形变速率；
对于所述配准后的SAR影像中时间间隔超过所述时间基线阈值的相邻影像，根据所述初始视线向形变速率和所述相邻影像的时间间隔确定所述相邻影像的替代差分干涉图；
将所述差分干涉图集和替代差分干涉图作为干涉点目标分析方法IPTA的标准数据输入，并选取包括CR点的相干点目标，对所述差分干涉图集和替代差分干涉图进行形变时序分析，得到所述相干点目标的形变速率和形变时间序列。
2.如权利要求1所述的方法，其特征在于，所述卫星雷达数据是覆盖同一滑坡区的不同模式数据，并且，所述不同模式数据的入射角基本一致。
3.如权利要求1所述的方法，其特征在于，根据距离方程R＝|Rs-Rt|，以及多普勒方程将所述GNSS设备记录的地理坐标转换为所述SAR主影像的距离多普勒坐标；
其中，R为传感器到点目标的距离，Rs为传感器的位置，vs为传感器的速度矢量，Rt为点目标的位置，vt为点目标的速度矢量，λ为雷达波长，fd为多普勒频率。
4.如权利要求1所述的方法，其特征在于，根据下述计算公式将GNSS设备观测值转换为雷达视线向观测值：
[Unsinα-Uecosα]sinθ+Uucosθ+δlos＝glos
其中，Ue、Un、Uu分别为GNSS设备东西、南北、垂向形变分量，α为卫星轨道方位角，θ为雷达入射角，glos为GNSS设备观测值转换得到的雷达视线向观测值。
5.如权利要求4所述的方法，其特征在于，根据所述雷达视线向形变时间序列计算理论雷达卫星拍摄时间间隔包括：
根据下述公式计算雷达视线向形变速率：
其中，t0为起始观测时间，ti为截止观测时间，Δvglos为起始观测时间t0到截止观测时间ti之间的雷达视线向形变速率；
将Δvglos代入下述计算公式中的Δv计算理论雷达卫星拍摄时间间隔：
其中，Δv为同一像元内雷达视线向形变速率变化值，λ为雷达波长，Δt为理论雷达卫星拍摄时间间隔。
6.如权利要求1所述的方法，其特征在于，根据下述形变速率计算公式计算初始视线向形变速率：
其中，Δtk为所述差分干涉图集中的第k组干涉像对的时间间隔，φk为所述差分干涉图集中的第k组干涉像对的解缠相位，N为所述干涉像对的总数值，p为所述初始视线向形变速率。
7.如权利要求6所述的方法，其特征在于，根据所述初始视线向形变速率和所述相邻影像的时间间隔确定所述相邻影像的替代差分干涉图包括：
所述初始视线向形变速率和所述相邻影像的时间间隔的乘积为相邻影像的替代视线向形变量；
根据所述替代视线向形变量确定所述相邻影像的替代差分干涉图。
8.如权利要求1所述的方法，其特征在于，所述方法还包括：
采用Goldstein-Werner滤波方法对所述差分干涉图集和替代差分干涉图进行滤波；
利用最小费用流算法对滤波后的差分干涉图进行空间相位解缠；
采用幅度离散指数和相干系数阈值综合选取SAR影像频谱特性稳定的所述相干点目标；
根据下述公式计算得到所有相干点目标的形变速率和形变时间序列：
其中，为第k个干涉像对中相干点目标i和j之间的相位差，θ为入射角，/>为垂直基线，tk为时间间隔，/>为邻点目标i和j间的相对高程差，/>为邻点目标i和j间的相对形变速率，/>为其它相位，λ为雷达波长，R为传感器到点目标的距离；
针对所述相干点目标的形变速率和形变时间序列，采用干涉点目标分析方法IPTA进行形变时序分析。
9.如权利要求1所述的方法，其特征在于，所述方法还包括：
根据下述公式计算同一时间段内CR点上InSAR测量值与GNSS转换为视线向形变观测值的互差的均方根误差σ：
其中，xi为CR点上InSAR测量值，yi为CR点上GNSS转换为视线向形变观测值，n为统计样本个数；
利用所述CR点上InSAR测量值与GNSS转换为视线向形变观测值的互差的均方根误差评判InSAR结果的测量精度。
10.一种计算机可读存储介质，所述计算机可读存储介质储存计算机可读指令，其特征在于，所述计算机可读指令由处理器执行以实现如权利要求1-9中任一项所述的方法。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311775084.1A CN117761716A (zh) | 2023-12-21 | 2023-12-21 | 融合CR和GNSS的InSAR滑坡形变时序监测方法和存储介质 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311775084.1A CN117761716A (zh) | 2023-12-21 | 2023-12-21 | 融合CR和GNSS的InSAR滑坡形变时序监测方法和存储介质 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN117761716A true CN117761716A (zh) | 2024-03-26 |
Family
ID=90311824
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202311775084.1A Pending CN117761716A (zh) | 2023-12-21 | 2023-12-21 | 融合CR和GNSS的InSAR滑坡形变时序监测方法和存储介质 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN117761716A (zh) |
-
2023
- 2023-12-21 CN CN202311775084.1A patent/CN117761716A/zh active Pending
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN113624122B (zh) | 融合GNSS数据与InSAR技术的桥梁变形监测方法 | |
| Catalão et al. | Merging GPS and atmospherically corrected InSAR data to map 3-D terrain displacement velocity | |
| CN109917356B (zh) | 一种机载激光扫描系统误差标定方法 | |
| CN105550639B (zh) | 对地观测激光测高卫星高程控制点自动提取方法和数据处理方法 | |
| CN101770027B (zh) | 基于InSAR与GPS数据融合的地表三维形变监测方法 | |
| CN102654576B (zh) | 基于sar图像和dem数据的图像配准方法 | |
| EP2413158B1 (en) | A method for monitoring terrain and man-made feature displacements using ground-based synthetic aperture radar (GBSAR) data | |
| CN110487241B (zh) | 卫星激光测高提取建筑区高程控制点方法 | |
| CN112711021B (zh) | 一种多分辨率InSAR交互干涉时序分析方法 | |
| CN105487065A (zh) | 一种时序星载雷达数据处理方法和装置 | |
| Zhang et al. | Satellite SAR geocoding with refined RPC model | |
| CN105824022A (zh) | 一种电网不良地质体三维形变监测方法 | |
| CN110109112A (zh) | 一种基于InSAR的填海地区机场形变监测方法 | |
| CN112986949B (zh) | 针对角反射器的sar高精度时序形变监测方法和装置 | |
| CN116047519A (zh) | 一种基于合成孔径雷达干涉测量技术的选点方法 | |
| CN108983231B (zh) | 一种基于视频合成孔径雷达的干涉视频测量方法 | |
| CN116774264B (zh) | 基于低轨卫星机会信号多普勒的运动目标定位方法 | |
| CN109946682B (zh) | 基于ICESat/GLAS的GF3数据基线估计方法 | |
| Feng et al. | An improved geometric calibration model for spaceborne SAR systems with a case study of large-scale Gaofen-3 images | |
| CN113238228B (zh) | 基于水准约束的三维地表形变获取方法、系统及装置 | |
| CN112179314B (zh) | 一种基于三维网格投影的多角度sar高程测量方法及系统 | |
| CN112946647A (zh) | 大气误差改正InSAR干涉图堆叠地质灾害普查方法和装置 | |
| WO1998002761A1 (en) | Terrain elevation measurement by interferometric synthetic aperture radar (ifsar) | |
| CN117761716A (zh) | 融合CR和GNSS的InSAR滑坡形变时序监测方法和存储介质 | |
| CN114527465A (zh) | 一种基于永久散射体的地面控制点自动选取方法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination |