CN113866765B - 基于多成分时间相干模型的PS-InSAR测量方法 - Google Patents

基于多成分时间相干模型的PS-InSAR测量方法 Download PDF

Info

Publication number
CN113866765B
CN113866765B CN202111121640.4A CN202111121640A CN113866765B CN 113866765 B CN113866765 B CN 113866765B CN 202111121640 A CN202111121640 A CN 202111121640A CN 113866765 B CN113866765 B CN 113866765B
Authority
CN
China
Prior art keywords
phase
time
unwrapping
gradient
deformation
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.)
Active
Application number
CN202111121640.4A
Other languages
English (en)
Other versions
CN113866765A (zh
Inventor
杨波
江利明
徐华平
周志伟
汪汉胜
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Institute of Precision Measurement Science and Technology Innovation of CAS
Original Assignee
Institute of Precision Measurement Science and Technology Innovation of CAS
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Institute of Precision Measurement Science and Technology Innovation of CAS filed Critical Institute of Precision Measurement Science and Technology Innovation of CAS
Priority to CN202111121640.4A priority Critical patent/CN113866765B/zh
Publication of CN113866765A publication Critical patent/CN113866765A/zh
Application granted granted Critical
Publication of CN113866765B publication Critical patent/CN113866765B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
    • G01S13/88—Radar or analogous systems specially adapted for specific applications
    • G01S13/89—Radar or analogous systems specially adapted for specific applications for mapping or imaging
    • G01S13/90—Radar or analogous systems specially adapted for specific applications for mapping or imaging using synthetic aperture techniques, e.g. synthetic aperture radar [SAR] techniques
    • G01S13/9021—SAR image post-processing techniques

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Physics & Mathematics (AREA)
  • Electromagnetism (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • General Physics & Mathematics (AREA)
  • Radar Systems Or Details Thereof (AREA)

Abstract

本发明公开了基于多成分时间相干模型的PS‑InSAR测量方法，初步筛选差分干涉相位堆栈数据中的PS候选点；拟合周期运动初始偏差常数；对PS候选点的差分干涉相位进行三维相位解缠；滤波分离大气与轨道相位；获得SAR观测地形的高程与形变结果；估计PS候选点的相干系数；迭代输出地形的高程与形变结果。本发明获得的监测点位密度高且耗时少；本发明的高程不确定度的像素占比、平均形变速率不确定度更优。

Description

基于多成分时间相干模型的PS-InSAR测量方法
技术领域
本发明属于干涉合成孔径雷达数据处理技术领域，具体基于多成分时间相干模型的PS-InSAR测量方法。
背景技术
合成孔径雷达(Synthetic Aperture Radar，SAR)是一种全天时、全天候、主动式微波成像雷达，对地表具有一定的穿透能力。在SAR基础上发展起来的干涉合成孔径雷达(SAR Interferometry,InSAR)是一种高效的、能够获取高程信息的前沿大地测量技术，可以高精度获取连续覆盖的地表高程与形变信息。永久散射体干涉合成孔径雷达(Persistent scatterers InSAR,PS-InSAR)是一种先进的InSAR后续发展技术，其测高精度可达亚米级、形变监测精度可达毫米/年，在高精度数字高程模型构建，火山、地震、滑坡监测，冰川冻土变化监测，城市沉降等地表精密测量中发挥了重要作用。
PS-InSAR方法基于时间相干模型最大化获取相邻像素的高程与形变速率梯度，进而逐个像素估计相干系数以筛选高相干PS像素。传统的时间相干模型假设残余高程相位与形变相位线性，然而，我国北方部分局部区域受地下水位季节性变化影响，地表高程呈周期性变化，夏季出现抬升、冬季出现沉降趋势；此外，在冰川冻土区域，地表高程也常出现冬季冻胀抬升、夏季融化沉降的季节性变化。因此，传统的线性时间相干模型无法准确表征具有周期形变信号的PS像素的时间相干性，导致较大周期形变的PS像素被剔除，较小周期形变的PS像素的结果误差较大。
发明内容
本发明针对地表高程季节性周期变化信号的场景，提出了基于多成分时间相干模型的PS-InSAR测量方法，达到了提高监测点密度与精度的目的。
本发明的上述目的通过以下技术方案实现：
基于多成分时间相干模型的PS-InSAR测量方法，还包括以下步骤：
步骤1，基于差分干涉相位堆栈数据，利用幅度离差准则初步筛选差分干涉相位堆栈数据中的PS候选点；
步骤2，根据月平均温度拟合周期运动初始偏差常数；
步骤3，将步骤2估计的周期运动初始偏差常数代入多成分时间相干模型，采用基于多成分时间相干模型的三维相位解缠方法，对步骤1筛选的PS候选点的差分干涉相位进行三维相位解缠，并计算时间相干系数、残余相位梯度、平均形变速率梯度、周期形变幅度梯度；
步骤4，去除PS候选点的差分干涉相位中的残余地形相位、线性形变相位、周期形变相位等模拟成分后，利用PS候选点的差分干涉相位中大气与轨道相位、 SAR系统热噪声相位、以及非模拟形变相位的不同时空特性，采用高通时间、低通空间滤波分离大气与轨道相位；
步骤5，将解缠的PS候选点的差分干涉相位减去大气与轨道相位，线性拟合残余地形、平均形变信息以及总的时序形变，进而获得SAR观测地形的高程与形变结果；
步骤6，将解缠的差分干涉相位减去大气与轨道相位、残余地形相位、以及时序形变相位，获得噪声相位，进而估计PS候选点的相干系数，
步骤7，若步骤6计算的相干系数大于相干阈值或迭代次数超过所设最大更新次数，则输出步骤5中地形的高程与形变结果；反之，剔除低于相干阈值的 PS候选点，返回步骤3，重复第3-6步骤。
如上所述的步骤1中幅度离差准则为
Figure BDA0003277439860000021
其中DA为幅度离差指数，
Figure BDA0003277439860000022
为设置的指数阈值。
如上所述的步骤3包括以下步骤：
步骤3.1，根据各个PS候选点在SAR图像的行列位置，构建空间三角网；
步骤3.2，对空间三角网边连接的两个PS候选点的相位差分，最大化多成分时间相干模型：
步骤3.3，在基线-时间二维平面构建约束三角网，约束主要由时间与基线约束构成，剔除基线-时间三角网中基线与时间大于设定阈值的三角形；
步骤3.4，计算第x个PS候选点的缠绕差分干涉相位梯度
Figure BDA0003277439860000023
的真实相位梯度
Figure BDA0003277439860000024
的相位梯度时间解缠值
Figure BDA0003277439860000025
步骤3.5，设置低相干阈值，依据相干阈值约束，剔除步骤3.2计算的相干系数小于低相干阈值对应的三角网边，在方位-距离二维平面构建高相干约束的最大连通三角网；
步骤3.6，依据步骤3.4获得的相位梯度时间解缠值
Figure BDA0003277439860000031
采用稀疏最小费用流方法，在时间维度解缠步骤3.4中真实相位梯度和对应的相位梯度时间解缠值的差值
Figure BDA0003277439860000032
的倍数
Figure BDA0003277439860000033
进而获得解缠后的真实相位梯度
Figure BDA0003277439860000034
采用基于稀疏最小费用流方法的优化模型获得真实相位梯度在时间与空间维的解缠值；
步骤3.7，对步骤3.6获得的解缠后的真实相位梯度
Figure BDA0003277439860000035
沿步骤3.5获得的最大连通三角网作为最优路径进行积分，从而获得缠绕的第x个PS候选点在第i 幅干涉图中的差分干涉相位
Figure BDA0003277439860000036
的解缠相位值ψx,i，根据步骤3.6获得的解缠后的真实相位梯度
Figure BDA0003277439860000037
估计值，在方位-距离平面展开三角网每条边的表达式，联立三角网所有边的表达式组成边的方程组，将步骤3.2中的时间相干系数作为边的表达式的权重，重加权最小二乘求解任意第x个PS候选点在第i幅干涉图中的差分干涉相位
Figure BDA0003277439860000038
的解缠相位值ψx,i。
如上所述的步骤3.2中最大化多成分时间相干模型基于以下公式：
Figure BDA0003277439860000039
获得任意相邻PS候选点的时间相干系数
Figure BDA00032774398600000310
以及残余相位梯度
Figure BDA00032774398600000311
平均形变速率梯度
Figure BDA00032774398600000312
周期形变幅度梯度
Figure BDA00032774398600000313
x1∈x，x2∈x，
Figure BDA00032774398600000314
为任意相邻PS候选点的差分算子，M表示干涉图幅数，j表示虚数单位，
Figure BDA00032774398600000315
表示第x个PS候选点在第i幅干涉图中的差分干涉相位，
Figure BDA00032774398600000316
与
Figure BDA00032774398600000317
分别表示第x个PS候选点在第i幅干涉图中的模拟残余地形、线性形变干涉相位以及周期形变干涉相位，
Figure BDA00032774398600000318
为第x个PS候选点的缠绕差分干涉相位梯度，
Figure BDA00032774398600000319
为第x个PS候选点在第i幅干涉图中的模拟残余地形差分相位梯度，
Figure BDA00032774398600000320
为第x个PS候选点在第i幅干涉图中的线性形变差分干涉相位梯度，
Figure BDA0003277439860000041
为第x个PS候选点在第i幅干涉图中的周期形变差分干涉相位梯度，其差分形式可表示为
Figure BDA0003277439860000042
Figure BDA0003277439860000043
Figure BDA0003277439860000044
其中，bi与ti分别表示第i幅干涉图的垂直基线及观测时间，Δh为残余高程， r表示雷达到目标的距离，v目标的线性形变速率，p表示周期幅度，θ与λ分别表示雷达入射角与雷达波长，t0表示周期运动初始偏差常数，如上所述的步骤 3.4中，第x个PS候选点的缠绕差分干涉相位梯度
Figure BDA0003277439860000045
的真实相位梯度
Figure BDA0003277439860000046
Figure BDA0003277439860000047
Figure BDA0003277439860000048
为残余相位梯度时间解缠的缠绕倍数，<·>-π，π表示相位缠绕算子，
基于扩展最小费用流方法，在时间维度对真实相位梯度进行相位解缠，并获得
Figure BDA0003277439860000049
的估计值
Figure BDA00032774398600000410
根据
Figure BDA00032774398600000411
的估计值
Figure BDA00032774398600000412
在时间维根据三角边展开整数线性方程：
Figure BDA00032774398600000413
联立三角网所有三角形的整数线性方程构成超定方程组，主干涉图先验
Figure BDA00032774398600000414
其中mast表示主干涉图编号，mast∈i，整数线性规划估计残余相位梯度时间解缠的缠绕倍数
Figure BDA00032774398600000415
的估计值
Figure BDA00032774398600000416
进而获得真实相位梯度
Figure BDA00032774398600000417
的相位梯度时间解缠值
Figure BDA00032774398600000418
如上所述的步骤3.6中，依据步骤3.4获得的相位梯度时间解缠值
Figure BDA00032774398600000419
采用稀疏最小费用流方法，在时间维度解缠步骤3.4中真实相位梯度和对应的相位梯度时间解缠值的差值
Figure BDA0003277439860000051
的倍数
Figure BDA0003277439860000052
进而获得解缠后的真实相位梯度
Figure BDA0003277439860000053
即
Figure BDA0003277439860000054
采用基于稀疏最小费用流方法的优化模型：
Figure BDA0003277439860000055
Figure BDA0003277439860000056
其中α、β、ξ为步骤3.5方位-距离平面上约束三角网的三条边，求解基于稀疏最小费用流方法的优化模型即可获得真实相位梯度在时间与空间维的解缠值。
本发明与现有技术相比，具有以下优点：
1、高密度。在相同PS-InSAR数据处理流程下，本发明的基于多成分时间相干模型获得的监测点位密度比基于二维线性模型提升了20％。
2、高精度。在相同PS-InSAR数据处理流程下，本发明的基于多成分时间相干模型获得0.5米以内的高程不确定度的像素比传统的基于二维的时间相干模型的占比多15.2％；本发明的基于多成分时间相干模型获得0.1毫米/年以内的平均形变速率不确定度的像素比传统的基于二维的时间相干模型的占比多 11.32％；此外，本发明的基于多成分时间相干模型的PS-InSAR方法获得了更准确的时序形变。
3、高效性。尽管比传统的二维时间相干模型增加了一维估计量，但对于相同的PS候选点集，本发明的基于多成分时间相干模型的方法的计算耗时比传统的基于二维的时间相干模型的方法只低9.5％。
附图说明
图1：基于多成分时间相干模型的PS-InSAR测量方法的流程图；
图2：基于多成分时间相干模型的PS-InSAR测量方法的三维相位解缠流程图；
图3：月平均温度与周期运动模型拟合曲线图；
图4：PS-InSAR估计结果。(a)和(b)分别为本发明估计的高程[米]图、平均形变速率[毫米/年]图，(c)、(d)、(e)分别为本发明估计的高程[米]、平均形变速率[毫米/年]、周期幅度[毫米]图；
图5：不确定度的统计结果。(a)和(b)分别为本发明估计的高程、平均形变速率的不确定度，(c)、(d)、(e)分别为本发明估计的高程、平均形变速率、周期幅度的不确定度；
图6：两种方法在共同PS点的时序形变估计结果。(a)和(b)分别为传统的基于二维的时间相干模型以及本发明的基于多成分相干模型的PS-InSAR方法估计的时序形变结果。
具体实施方式
为了便于本领域普通技术人员理解和实施本发明，下面结合实施例对本发明作进一步的详细描述，应当理解，此处所描述的实施示例仅用于说明和解释本发明，并不用于限定本发明。
实施例：
如图1所示，一种基于多成分时间相干模型的PS-InSAR高程与形变测量方法，包括以下步骤：
步骤1，基于差分干涉相位堆栈数据，利用幅度离差准则
Figure BDA0003277439860000061
初步筛选差分干涉相位堆栈数据中的PS候选点，其中DA为幅度离差指数，
Figure BDA0003277439860000062
为设置的指数阈值；
步骤2，根据月平均温度拟合周期运动初始偏差常数；
步骤3，将步骤2估计的周期运动初始偏差常数代入多成分时间相干模型，采用基于多成分时间相干模型的三维相位解缠方法，对步骤1筛选的PS候选点的差分干涉相位进行三维相位解缠，并计算时间相干系数、残余相位梯度、平均形变速率梯度、周期形变幅度梯度；
步骤4，去除PS候选点的差分干涉相位中的残余地形相位、线性形变相位、周期形变相位等模拟成分后，利用PS候选点的差分干涉相位中大气与轨道相位、 SAR系统热噪声相位、以及非模拟形变相位的不同时空特性，采用高通时间、低通空间滤波分离大气与轨道相位；
步骤5，将解缠的PS候选点的差分干涉相位减去大气与轨道相位，线性拟合残余地形、平均形变信息以及总的时序形变，进而获得SAR观测地形的高程与形变结果；
步骤6，将解缠的差分干涉相位减去大气与轨道相位、残余地形相位、以及时序形变相位，获得噪声相位，进而估计PS候选点的相干系数。
步骤7，若步骤6计算的相干系数大于相干阈值或迭代次数超过所设最大更新次数，则输出步骤5中地形的高程与形变结果；反之，剔除低于相干阈值的 PS候选点，返回步骤3，重复第3-6步骤。
如上所述的步骤3中，采用图2中基于多成分时间相干模型的三维相位解缠方法，对步骤1筛选的PS候选点的差分干涉相位进行三维相位解缠，并计算时间相干系数、残余相位梯度、平均形变速率梯度、周期形变幅度梯度的方法为：
步骤3.1，根据各个PS候选点在SAR图像的行列位置，构建空间三角网；
步骤3.2，对空间三角网边连接的两个PS候选点的相位差分，最大化多成分时间相干模型：
Figure BDA0003277439860000071
获得任意相邻PS候选点(x1与x2)的时间相干系数
Figure BDA0003277439860000072
以及残余相位梯度
Figure BDA0003277439860000073
平均形变速率梯度
Figure BDA0003277439860000074
周期形变幅度梯度
Figure BDA0003277439860000075
x1∈x，x2∈x。
式(1)中
Figure BDA0003277439860000076
为任意相邻PS候选点的差分算子，M表示干涉图幅数，j表示虚数单位，
Figure BDA0003277439860000077
表示第x个PS候选点在第i幅干涉图中的差分干涉相位，ΔφH,x,i、
Figure BDA0003277439860000078
与
Figure BDA0003277439860000079
分别表示第x个PS候选点在第i幅干涉图中的模拟残余地形、线性形变干涉相位以及周期形变干涉相位，
Figure BDA00032774398600000710
为第x个PS候选点的缠绕差分干涉相位梯度，
Figure BDA00032774398600000711
为第x个PS候选点在第i幅干涉图中的模拟残余地形差分相位，
Figure BDA00032774398600000712
为第x个PS候选点在第i幅干涉图中的线性形变差分干涉相位梯度，
Figure BDA00032774398600000713
为第x个PS候选点在第i幅干涉图中的周期形变差分干涉相位梯度，其差分形式可表示为
Figure BDA0003277439860000081
其中，bi与ti分别表示第i幅干涉图的垂直基线及观测时间，Δh为残余高程， r表示雷达到目标的距离，v目标的线性形变速率，p表示周期幅度，θ与λ分别表示雷达入射角与雷达波长，t0表示周期运动初始偏差常数。该优化问题可在 CPU并行及GPU计算策略下用网格搜索方法快速求解。
步骤3.3，在基线-时间二维平面构建约束三角网，约束主要由时间与基线约束构成，剔除基线-时间三角网中基线与时间大于设定阈值的三角形；
步骤3.4，假设第x个PS候选点的缠绕差分干涉相位梯度
Figure BDA0003277439860000082
的真实相位梯度
Figure BDA0003277439860000083
而
Figure BDA0003277439860000084
表示模拟相位与残余缠绕相位之和，
Figure BDA0003277439860000085
为残余相位梯度时间解缠的缠绕倍数，其中<·>-π，π表示相位缠绕算子。基于扩展最小费用流方法，在时间维度对真实相位梯度进行相位解缠，并获得
Figure BDA0003277439860000086
的估计值
Figure BDA0003277439860000087
基于扩展最小费用流方法的模型为：
Figure BDA0003277439860000088
Figure BDA0003277439860000089
其中，
Figure BDA00032774398600000810
表示以最小化{}内函数为目标估计
Figure BDA00032774398600000811
为自变量，round[·]表示整数近似算子，Z为整数，
Figure BDA00032774398600000812
表示任意第i1幅与第i2幅干涉图的差分算子，i1和i2∈i，
Figure BDA0003277439860000091
α、β、ξ为步骤3.3中PS候选点在基线-时间平面上约束三角网的三角形的三条边，Δα、Δβ、Δξ分别表示三角形三条边上PS候选点的差分算子。
根据式(3)
Figure BDA0003277439860000092
的估计值
Figure BDA0003277439860000093
在时间维根据三角边展开整数线性方程：
Figure BDA0003277439860000094
联立三角网所有三角形的线性方程(4)构成超定方程组，由主干涉图先验
Figure BDA0003277439860000095
其中mast表示主干涉图编号，mast∈i，整数线性规划估计残余相位梯度时间解缠的缠绕倍数
Figure BDA0003277439860000096
的估计值
Figure BDA0003277439860000097
进而获得真实相位梯度
Figure BDA0003277439860000098
的相位梯度时间解缠值
Figure BDA0003277439860000099
步骤3.5，设置低相干阈值，依据相干阈值约束，剔除步骤3.2计算的相干系数小于低相干阈值对应的三角网边，在方位-距离二维平面构建高相干约束的最大连通三角网；
步骤3.6，依据步骤3.4获得的相位梯度时间解缠值
Figure BDA00032774398600000910
采用稀疏最小费用流方法，在时间维度解缠步骤3.4中真实相位梯度和对应的相位梯度时间解缠值的差值
Figure BDA00032774398600000911
的倍数
Figure BDA00032774398600000912
进而获得解缠后的真实相位梯度
Figure BDA00032774398600000913
即
Figure BDA00032774398600000914
采用基于稀疏最小费用流方法的优化模型：
Figure BDA00032774398600000915
其中α、β、ξ为步骤3.5方位-距离平面上约束三角网的三条边。求解基于稀疏最小费用流方法的优化模型(5)即可获得真实相位梯度在时间与空间维的解缠值。
步骤3.7，对步骤3.6获得的解缠后的真实相位梯度
Figure BDA0003277439860000101
沿步骤3.5获得的最大连通三角网作为最优路径进行积分，从而获得缠绕的第x个PS候选点在第i 幅干涉图中的差分干涉相位
Figure BDA0003277439860000102
的解缠相位值ψx,i。根据式(5)获得的真实相位梯度在时间与空间维的解缠值，在方位-距离平面展开三角网每条边的表达式
Figure BDA0003277439860000103
联立所有边的表达式组成边的方程组，将步骤3.2中的时间相干系数作为边的表达式的权重，重加权最小二乘求解任意第x个PS候选点在第i幅干涉图中的差分干涉相位
Figure BDA0003277439860000104
的解缠相位值ψx,i。
通过TerraSAR-X/TanDEM-X获取的31幅首都国际机场区域的真实数据，验证本发明一种基于多成分时间相干模型的PS-InSAR高程与形变测量方法的有效性。31幅TerraSAR-X/TanDEM-X图像获取时间、基线参数及观测时刻月平均温度信息如表1所示，利用月平均温度与周期运动模型拟合，如图3所示，可估计出多成分时间相干模型的初始偏差常数t0＝-0.483年。
表1为31幅TerraSAR-X/TanDEM-X图像获取时间与基线参数
Figure BDA0003277439860000105
Figure BDA0003277439860000111
将本发明的处理结果与传统的基于线性时间相干模型的PS-InSAR方法的结果对比，说明本发明所提方法的优势。表2、3中分别列出了基于两种相干模型在三次迭代中筛选的PS点数目及三次迭代的计算耗时。从表2可以看出，在初始PS候选点集同为343040个的情况下，最终多成分时间相干模型筛选的高相干点数比线性时间相干模型提升了20％。从表3可以看出，在第一迭代、相同的 PS候选点集下，尽管比二维线性模型增加了一维估计量，基于多成分时间相干模型的方法的计算耗时比二维线性模型的方法只低9.5％。
表2为三次迭代后的PS候选点数目
Figure BDA0003277439860000112
表3为三次迭代计算耗时
Figure BDA0003277439860000113
在相同PS-InSAR数据处理流程下，从图4基于两种时间相干模型的PS-InSAR估计结果可以看出，两种时间相干模型下获得的高程与平均形变速率总体分布趋势相同，本发明提供了额外的周期形变幅度。计算高程、平均形变速率、周期幅度成分的不确定度，统计结果如图5所示，可以看出，本发明的基于多成分时间相干模型获得0.5米以内的高程不确定度的像素比传统的基于二维的线性时间相干模型的占比多15.2％；本发明的基于多成分时间相干模型获得0.1 毫米/年以内的平均形变速率不确定度的像素比传统的基于二维的线性时间相干模型的占比多11.32％。图6画出了两种方法在任选的共同PS点的时序形变曲线图，从图中可以看出，传统的基于二维的线性时间相干模型的PS-InSAR方法估计的第6、7个时段的形变量与模型存在较大偏差，而本发明的基于多成分时间相干模型的PS-InSAR方法获得了更准确的时序形变。
以上仅是本发明的优选实施方式，本发明的保护范围并不仅局限于上述实施例，凡属于本发明思路下的技术方案均属于本发明的保护范围。应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理前提下的若干改进和润饰，应视为本发明的保护范围。

Claims (5)

1.基于多成分时间相干模型的PS-InSAR测量方法，其特征在于，还包括以下步骤：
步骤1，基于差分干涉相位堆栈数据，利用幅度离差准则初步筛选差分干涉相位堆栈数据中的PS候选点；
步骤2，根据月平均温度拟合周期运动初始偏差常数；
步骤3，将步骤2估计的周期运动初始偏差常数代入多成分时间相干模型，采用基于多成分时间相干模型的三维相位解缠方法，对步骤1筛选的PS候选点的差分干涉相位进行三维相位解缠，并计算时间相干系数、残余相位梯度、平均形变速率梯度、周期形变幅度梯度；
步骤4，去除PS候选点的差分干涉相位中的残余地形相位、线性形变相位、周期形变相位等模拟成分后，利用PS候选点的差分干涉相位中大气与轨道相位、SAR系统热噪声相位、以及非模拟形变相位的不同时空特性，采用高通时间、低通空间滤波分离大气与轨道相位；
步骤5，将解缠的PS候选点的差分干涉相位减去大气与轨道相位，线性拟合残余地形、平均形变信息以及总的时序形变，进而获得SAR观测地形的高程与形变结果；
步骤6，将解缠的差分干涉相位减去大气与轨道相位、残余地形相位、以及时序形变相位，获得噪声相位，进而估计PS候选点的相干系数；
步骤7，若步骤6计算的相干系数大于相干阈值或迭代次数超过所设最大更新次数，则输出步骤5中地形的高程与形变结果；反之，剔除低于相干阈值的PS候选点，返回步骤3，重复第3-6步骤；
所述的步骤3包括以下步骤：
步骤3.1，根据各个PS候选点在SAR图像的行列位置，构建空间三角网；
步骤3.2，对空间三角网边连接的两个PS候选点的相位差分，最大化多成分时间相干模型：
步骤3.3，在基线-时间二维平面构建约束三角网，约束主要由时间与基线约束构成，剔除基线-时间三角网中基线与时间大于设定阈值的三角形；
步骤3.4，计算第x个PS候选点的缠绕差分干涉相位梯度
Figure FDA0003907853360000011
的真实相位梯度
Figure FDA0003907853360000021
的相位梯度时间解缠值
Figure FDA0003907853360000022
步骤3.5，设置低相干阈值，依据相干阈值约束，剔除步骤3.2计算的相干系数小于低相干阈值对应的三角网边，在方位-距离二维平面构建高相干约束的最大连通三角网；
步骤3.6，依据步骤3.4获得的相位梯度时间解缠值
Figure FDA0003907853360000023
采用稀疏最小费用流方法，在时间维度解缠步骤3.4中真实相位梯度和对应的相位梯度时间解缠值的差值
Figure FDA0003907853360000024
的倍数
Figure FDA0003907853360000025
进而获得解缠后的真实相位梯度
Figure FDA0003907853360000026
采用基于稀疏最小费用流方法的优化模型获得真实相位梯度在时间与空间维的解缠值；
步骤3.7，对步骤3.6获得的解缠后的真实相位梯度
Figure FDA0003907853360000027
沿步骤3.5获得的最大连通三角网作为最优路径进行积分，从而获得缠绕的第x个PS候选点在第i幅干涉图中的差分干涉相位
Figure FDA0003907853360000028
的解缠相位值ψx,i，根据步骤3.6获得的解缠后的真实相位梯度
Figure FDA0003907853360000029
估计值，在方位-距离平面展开三角网每条边的表达式，联立三角网所有边的表达式组成边的方程组，将步骤3.2中的时间相干系数作为边的表达式的权重，重加权最小二乘求解任意第x个PS候选点在第i幅干涉图中的差分干涉相位
Figure FDA00039078533600000210
的解缠相位值ψx,i。
2.根据权利要求1所述的基于多成分时间相干模型的PS-InSAR测量方法，其特征在于，所述的步骤1中幅度离差准则为
Figure FDA00039078533600000211
其中DA为幅度离差指数，
Figure FDA00039078533600000212
为设置的指数阈值。
3.根据权利要求1所述的基于多成分时间相干模型的PS-InSAR测量方法，其特征在于，所述的步骤3.2中最大化多成分时间相干模型基于以下公式：
Figure FDA00039078533600000213
获得任意相邻PS候选点的时间相干系数
Figure FDA00039078533600000214
以及残余相位梯度
Figure FDA00039078533600000215
平均形变速率梯度
Figure FDA00039078533600000216
周期形变幅度梯度
Figure FDA00039078533600000217
x1∈x，x2∈x，
Figure FDA00039078533600000218
为任意相邻PS候选点的差分算子，M表示干涉图幅数，j表示虚数单位，
Figure FDA0003907853360000031
表示第x个PS候选点在第i幅干涉图中的差分干涉相位，△φH,x,i、
Figure FDA0003907853360000032
与
Figure FDA0003907853360000033
分别表示第x个PS候选点在第i幅干涉图中的模拟残余地形、线性形变干涉相位以及周期形变干涉相位，
Figure FDA0003907853360000034
为第x个PS候选点的缠绕差分干涉相位梯度，
Figure FDA0003907853360000035
为第x个PS候选点在第i幅干涉图中的模拟残余地形差分相位梯度，
Figure FDA0003907853360000036
为第x个PS候选点在第i幅干涉图中的线性形变差分干涉相位梯度，
Figure FDA0003907853360000037
为第x个PS候选点在第i幅干涉图中的周期形变差分干涉相位梯度，其差分形式可表示为：
Figure FDA0003907853360000038
Figure FDA0003907853360000039
Figure FDA00039078533600000310
其中，bi与ti分别表示第i幅干涉图的垂直基线及观测时间，△h为残余高程，r表示雷达到目标的距离，v目标的线性形变速率，p表示周期幅度，θ与λ分别表示雷达入射角与雷达波长，t0表示周期运动初始偏差常数。
4.根据权利要求3所述的基于多成分时间相干模型的PS-InSAR测量方法，其特征在于，所述的步骤3.4中，第x个PS候选点的缠绕差分干涉相位梯度
Figure FDA00039078533600000311
的真实相位梯度
Figure FDA00039078533600000312
Figure FDA00039078533600000313
Figure FDA00039078533600000314
为残余相位梯度时间解缠的缠绕倍数，<·>-π，π表示相位缠绕算子，基于扩展最小费用流方法，在时间维度对真实相位梯度进行相位解缠，并获得
Figure FDA00039078533600000315
的估计值
Figure FDA00039078533600000316
根据
Figure FDA00039078533600000317
的估计值
Figure FDA00039078533600000318
在时间维根据三角边展开整数线性方程：
Figure FDA00039078533600000319
联立三角网所有三角形的整数线性方程构成超定方程组，主干涉图先验
Figure FDA0003907853360000041
其中，
Figure FDA0003907853360000042
为任意相邻PS候选点的差分算子，
Figure FDA0003907853360000043
表示任意第i1幅与第i2幅干涉图的差分算子，i1和i2∈i，mast表示主干涉图编号，mast∈i，整数线性规划估计残余相位梯度时间解缠的缠绕倍数
Figure FDA0003907853360000044
的估计值
Figure FDA0003907853360000045
进而获得真实相位梯度
Figure FDA0003907853360000046
的相位梯度时间解缠值
Figure FDA0003907853360000047
5.根据权利要求4所述的基于多成分时间相干模型的PS-InSAR测量方法，其特征在于，所述的步骤3.6中，依据步骤3.4获得的相位梯度时间解缠值
Figure FDA0003907853360000048
采用稀疏最小费用流方法，在时间维度解缠步骤3.4中真实相位梯度和对应的相位梯度时间解缠值的差值
Figure FDA0003907853360000049
的倍数
Figure FDA00039078533600000410
进而获得解缠后的真实相位梯度
Figure FDA00039078533600000411
即
Figure FDA00039078533600000412
采用基于稀疏最小费用流方法的优化模型：
Figure FDA00039078533600000413
Figure FDA00039078533600000414
其中，△α、△β、△ξ分别表示三角形三条边上PS候选点的差分算子，α、β、ξ为步骤3.5方位-距离平面上约束三角网的三条边，求解基于稀疏最小费用流方法的优化模型即可获得真实相位梯度在时间与空间维的解缠值。
CN202111121640.4A 2021-09-24 2021-09-24 基于多成分时间相干模型的PS-InSAR测量方法 Active CN113866765B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111121640.4A CN113866765B (zh) 2021-09-24 2021-09-24 基于多成分时间相干模型的PS-InSAR测量方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111121640.4A CN113866765B (zh) 2021-09-24 2021-09-24 基于多成分时间相干模型的PS-InSAR测量方法

Publications (2)

Publication Number Publication Date
CN113866765A CN113866765A (zh) 2021-12-31
CN113866765B true CN113866765B (zh) 2022-12-13

Family

ID=78993904

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111121640.4A Active CN113866765B (zh) 2021-09-24 2021-09-24 基于多成分时间相干模型的PS-InSAR测量方法

Country Status (1)

Country Link
CN (1) CN113866765B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115267774A (zh) * 2022-06-14 2022-11-01 深圳大学 一种城区多时相InSAR相位解缠方法、终端及存储介质
CN114966692B (zh) * 2022-07-19 2022-11-08 之江实验室 基于Transformer的InSAR技术冻土区多变量时序形变预测方法及装置

Citations (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6583751B1 (en) * 1999-05-25 2003-06-24 Politecnico Di Milano Process for radar measurements of the movement of city areas and landsliding zones
CN103822598A (zh) * 2014-02-26 2014-05-28 北京理工大学 地基sar在时间去相关严重区域的形变监测方法
CN104091064A (zh) * 2014-07-02 2014-10-08 北京航空航天大学 基于优化解空间搜索法的PS-DInSAR地表形变测量参数估计方法
CN106772342A (zh) * 2017-01-11 2017-05-31 西南石油大学 一种适用于大梯度地表沉降监测的时序差分雷达干涉方法
CN108627832A (zh) * 2018-05-11 2018-10-09 电子科技大学 一种基于多时序sar图像提取输电通道地表形变的方法
CN109031301A (zh) * 2018-09-26 2018-12-18 云南电网有限责任公司电力科学研究院 基于PSInSAR技术的山区地形形变提取方法
CN109541593A (zh) * 2018-10-30 2019-03-29 北京航空航天大学 一种改进的最小费用流InSAR相位解缠方法
CN110109112A (zh) * 2019-04-30 2019-08-09 成都理工大学 一种基于InSAR的填海地区机场形变监测方法
CN110412574A (zh) * 2019-09-05 2019-11-05 河海大学 一种时空相干性增强的分布式目标InSAR时序处理方法和装置
CN112986993A (zh) * 2021-02-07 2021-06-18 同济大学 一种基于空间约束的InSAR形变监测方法

Patent Citations (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6583751B1 (en) * 1999-05-25 2003-06-24 Politecnico Di Milano Process for radar measurements of the movement of city areas and landsliding zones
CN103822598A (zh) * 2014-02-26 2014-05-28 北京理工大学 地基sar在时间去相关严重区域的形变监测方法
CN104091064A (zh) * 2014-07-02 2014-10-08 北京航空航天大学 基于优化解空间搜索法的PS-DInSAR地表形变测量参数估计方法
CN106772342A (zh) * 2017-01-11 2017-05-31 西南石油大学 一种适用于大梯度地表沉降监测的时序差分雷达干涉方法
CN108627832A (zh) * 2018-05-11 2018-10-09 电子科技大学 一种基于多时序sar图像提取输电通道地表形变的方法
CN109031301A (zh) * 2018-09-26 2018-12-18 云南电网有限责任公司电力科学研究院 基于PSInSAR技术的山区地形形变提取方法
CN109541593A (zh) * 2018-10-30 2019-03-29 北京航空航天大学 一种改进的最小费用流InSAR相位解缠方法
CN110109112A (zh) * 2019-04-30 2019-08-09 成都理工大学 一种基于InSAR的填海地区机场形变监测方法
CN110412574A (zh) * 2019-09-05 2019-11-05 河海大学 一种时空相干性增强的分布式目标InSAR时序处理方法和装置
CN112986993A (zh) * 2021-02-07 2021-06-18 同济大学 一种基于空间约束的InSAR形变监测方法

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
A triangle-oriented Spatial–Temporal phase unwrapping algorithm based on irrotational constraints for time-series InSAR;Rui Li等;《IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING》;20191231;第57卷(第12期);第10263-10275页 *
New advances of the extended minimum cost flow phase unwrapping algorithm for SBAS-DInSAR analysis at full spatial resolution;Antonio Pepe等;《IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING》;20111031;第49卷(第10期);第4062-4079页 *
PSInSAR在地面沉降监测中的研究分析;刘欢欢;《中国优秀硕士学位论文全文数据库信息科技辑(月刊)》;20100815(第8期);第I136-404页 *
PS-InSAR技术在城市地表沉降监测的应用探索;胡文;《2019年12月建筑科技与管理学术交流会论文集》;20191220;第95-100页 *

Also Published As

Publication number Publication date
CN113866765A (zh) 2021-12-31

Similar Documents

Publication Publication Date Title
CN113866765B (zh) 基于多成分时间相干模型的PS-InSAR测量方法
Kolecka et al. Assessment of the accuracy of SRTM C-and X-Band high mountain elevation data: A case study of the Polish Tatra Mountains
CN111059998B (zh) 一种基于高分辨率的时序InSAR形变监测方法及系统
CN105954747A (zh) 基于电网不良地质体三维形变监测的塔基稳定性分析方法
CN109031301A (zh) 基于PSInSAR技术的山区地形形变提取方法
CN104123464A (zh) 一种高分辨率InSAR时序分析反演地物高程与地面沉降量的方法
CN111208512A (zh) 一种基于视频合成孔径雷达的干涉测量方法
Booth et al. Multi-year, three-dimensional landslide surface deformation from repeat lidar and response to precipitation: Mill Gulch earthflow, California
CN110456346A (zh) 一种基于InSAR技术的输电铁塔倾斜监测方法
CN105824022A (zh) 一种电网不良地质体三维形变监测方法
CN108663678B (zh) 基于混合整数优化模型的多基线InSAR相位解缠算法
CN108919266A (zh) 一种基于PSInSAR技术的桥梁安全预警方法
CN116338607B (zh) 时间域和空间域两步式InSAR对流层延迟矫正方法
CN116299455A (zh) 基于PSInSAR与SqueeSAR的设施形变分析方法
CN112685819A (zh) 一种用于大坝及滑坡变形gb-sar监测的数据后处理方法及系统
CN113238228B (zh) 基于水准约束的三维地表形变获取方法、系统及装置
Feng et al. Automatic selection of permanent scatterers-based GCPs for refinement and reflattening in InSAR DEM generation
CN116977580A (zh) 一种基于机载LiDAR的山区大比例尺DEM制作方法
CN115546264A (zh) 一种星载InSAR图像精细配准与立体测量方法
Deo et al. Evaluation of interferometric SAR DEMs generated using TanDEM-X data
Pérez-Falls et al. Land subsidence in Villahermosa Tabasco Mexico, using radar interferometry
Ferretti et al. Monitoring terrain deformations using multi-temporal SAR images
CN114114258A (zh) 一种基于SBAS-InSAR技术的输电铁塔形变监测方法
Gao et al. Accuracy comparison and analysis of interpolation methods in DEM generation with 3D laser point cloud data
Bagheri et al. Urban TanDEM-X Raw DEM Fusion Based ON TV-L1 and Huber Models

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
GR01 Patent grant
GR01 Patent grant