CN112797886A - 面向缠绕相位的InSAR时序三维形变监测方法 - Google Patents

面向缠绕相位的InSAR时序三维形变监测方法 Download PDF

Info

Publication number
CN112797886A
CN112797886A CN202110111179.8A CN202110111179A CN112797886A CN 112797886 A CN112797886 A CN 112797886A CN 202110111179 A CN202110111179 A CN 202110111179A CN 112797886 A CN112797886 A CN 112797886A
Authority
CN
China
Prior art keywords
insar
time sequence
dimensional
target arc
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.)
Granted
Application number
CN202110111179.8A
Other languages
English (en)
Other versions
CN112797886B (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.)
Central South University
Original Assignee
Central South University
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 Central South University filed Critical Central South University
Priority to CN202110111179.8A priority Critical patent/CN112797886B/zh
Publication of CN112797886A publication Critical patent/CN112797886A/zh
Application granted granted Critical
Publication of CN112797886B publication Critical patent/CN112797886B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01B—MEASURING LENGTH, THICKNESS OR SIMILAR LINEAR DIMENSIONS; MEASURING ANGLES; MEASURING AREAS; MEASURING IRREGULARITIES OF SURFACES OR CONTOURS
    • G01B7/00—Measuring arrangements characterised by the use of electric or magnetic techniques
    • G01B7/16—Measuring arrangements characterised by the use of electric or magnetic techniques for measuring the deformation in a solid, e.g. by resistance strain gauge
    • 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
    • G01S13/9023—SAR image post-processing techniques combined with interferometric techniques
    • 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
    • G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
    • G01S7/02—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00
    • G01S7/41—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00 using analysis of echo signal for target characterisation; Target signature; Target cross-section

Landscapes

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

Abstract

本发明属于干涉合成孔径雷达技术领域，具体涉及面向缠绕相位的InSAR时序三维形变监测方法。包括：收集多轨道时序SAR影像，并获取干涉对数据集和时序高相干点数据集，并构建狄洛尼三角网，基于地表应力应变模型建立多轨道时序InSAR缠绕相位与时序三维地表形变梯度之间的函数关系；采用迭代加权最小二乘方法求解时序三维地表形变梯度，得到目标弧段时序三维地表形变，通过空间积分获得每个高相干点上的时序三维地表形变，该方法避免了地表形变求解过程中干涉图解缠步骤，并降低了粗差及误差较大的观测值对未知参数解算的影响，InSAR三维地表形变监测过程高效、准确。

Description

面向缠绕相位的InSAR时序三维形变监测方法
技术领域
本发明属于干涉合成孔径雷达技术领域，具体涉及面向缠绕相位的InSAR时序三维形变监测方法。
背景技术
合成孔径雷达干涉(Interferometric Synthetic Aperture Radar,SAR，InSAR)技术利用两景SAR影像可快速获取大范围、高分辨率的地表形变结果。然而，单个轨道的InSAR数据仅能获取真实三维形变沿雷达卫星视线向的一维投影，难以准确揭示真实地表运动特征，极易造成错判、误判。随着SAR卫星的不断增多，同一研究区域可获取多轨道/不同技术得到的不同方向InSAR时序形变观测数据，为实现真实时序三维形变的反演计算提供了契机。由于原始InSAR观测数据是通过2景单视复数影像(Single Look Complex,SLC)共轭相乘得到，因此，代表地表形变数据的InSAR相位观测值全部介于[-π，π]之间，即InSAR相位是缠绕的。在选定一个稳定参考点的情况下，将[-π，π]之间的InSAR缠绕相位恢复到实数域的过程称之为解缠。目前的InSAR时序三维形变监测方法均需要在InSAR干涉图正确解缠的前提下，才能得到可靠的三维形变结果。
然而，InSAR干涉相位解缠本身是一个秩亏的问题，该过程对失相干噪声较为敏感，且极易出现解缠误差。随着技术的发展，SAR影像的空间分辨率和时间采样率越来越高。海量SAR数据为InSAR时序三维形变监测提供数据基础的同时，也使得InSAR相位解缠更具挑战性。对于普通用户而言，受计算机硬件设备的限制，甚至无法直接对大范围、长时序的InSAR干涉图进行相位解缠操作。
发明内容
基于此，本发明针对传统方法在进行InSAR三维地表形变监测时相位解缠困难的问题，基于地表应力应变模型，通过利用迭代加权最小二乘方法求解三维形变梯度的方式，实现InSAR三维地表形变监测。
本发明提供了一种面向缠绕相位的InSAR时序三维形变监测方法，具体包括：
S1:获取多轨道时序SAR影像，将同一轨道SAR影像进行配准与差分干涉，获得预设时空基线阈值的多轨道InSAR干涉对数据集；
S2:采用振幅离差法获取多轨道SAR影像中的时序高相干点数据集；并基于所述高相干点构建狄洛尼三角网，经预处理后获得多个目标弧段；
S3:基于地表应力应变模型，建立目标弧段处预设范围内所有弧段上的InSAR相位差与目标弧段上时序三维地表变形梯度的函数关系；
S4:采用迭代加权最小二乘方法计算目标弧段时序三维地表变形梯度，并根据所述目标弧段所对应的坐标增量，获得目标弧段所对应的三维地表形变；
S5:获取所有目标弧段的三维地表形变，对每个时刻所有弧段上的三维地表形变依次进行空间积分，获得所有高相干点上的时序三维地表形变结果。
进一步的，所述预设时空基线阈值分别为100天和200m。
进一步的，所述步骤S2中的采用振幅离差法获取多轨道SAR影像中的时序高相干点数据集步骤具体包括：
计算所述SAR影像像素点的后向散射强度的标准差σA与平均值mA的比值
Figure BDA0002918935860000031
选取DA<0.4的像素点作为高相干点，获取高相干点数据集。
进一步的，所述步骤S2中的预处理具体为：
基于构建的狄洛尼三角网和高相干点集，获取所述狄洛尼三角网中每一个弧段中点的坐标和所有高相干点的边界多边形；
将中点位于边界多边形外的对应弧段剔除，剩余的有效弧段即为目标弧段。
进一步的，所述步骤S3具体包括：
将目标弧段预设范围内的高相干点处的InSAR观测相位与目标弧段的起点处的InSAR观测相位进行做差，获得目标弧段预设范围内的高相干点与目标弧段起点的InSAR相位差；
基于地表应力应变模型和卫星的方位角及入射角，获得单个InSAR干涉图中，目标弧段预设范围内的InSAR相位差与目标弧段处三维形变梯度之间的函数关系Bj,sm；
根据所述干涉对数据集与SAR影像的构网关系，获取干涉对相位与时序相位之间的系数矩阵Bsbl，j；
根据所述目标弧段预设范围内的高相干点与目标弧段起点的InSAR相位差、单个干涉图预设范围内弧段上的InSAR相位差与目标弧段处三维形变梯度之间的函数关系Bj,sm、以及干涉对相位与时序相位之间的系数矩阵Bsbl，j，构建在多个轨道InSAR时序干涉图中，预设范围内弧段上的InSAR相位差与时序三维形变梯度之间的函数关系为：
ΔLts＝Bts·lts
所述ΔLts＝[(ΔL1,ts)T,(ΔL2,ts)T,...,(ΔLj,ts)T,...,(ΔLJ,ts)T]T，ΔLts＝[(ΔL1，ts)T，(ΔL2，ts)T，...，(ΔLj，ts)T，...，(ΔLJ，ts)T]T,Bts＝[(B1，ts)T，(B2，ts)T，...，(Bj，ts)T，...，(BJ，ts)T]T，
Figure BDA0002918935860000041
代表第j个轨道中第i个InSAR干涉对所对应的Kj，i个Pk与目标弧段起点P0组成的弧段上的InSAR相位差
Figure BDA0002918935860000042
Figure BDA0002918935860000043
代表系数矩阵，
Figure BDA0002918935860000044
为克罗内克积符号，
Figure BDA0002918935860000045
代表时序三维地表形变梯度，
Figure BDA0002918935860000046
Figure BDA0002918935860000047
为第s个SAR影像获取时刻对应的累积三维地表形变梯度，下标ts代表时序。
进一步的，所述基于地表应力应变模型和卫星的方位角及入射角，获得单个InSAR干涉图中，预设范围内弧段上的InSAR相位差与目标弧段处三维形变梯度之间的函数关系Bj,sm，具体包括：
所述地表应力应变模型为dk＝H×Δxk+d0，所述dk、d0为预设范围内高相干点处和目标弧段起点的三维形变，H为目标弧段起点的三维形变梯度矩阵，Δxk＝xk-x0为预设范围内高相干点处和目标弧段起点的三维坐标之差；
所述预设范围内高相干点处和目标弧段起点的InSAR观测数据分别为
Figure BDA0002918935860000048
所对应的卫星方位角和入射角分别为
Figure BDA0002918935860000049
和
Figure BDA00029189358600000410
Figure BDA00029189358600000411
和
Figure BDA00029189358600000412
根据SAR卫星成像几何可得函数关系：
Figure BDA00029189358600000413
式中，
Figure BDA00029189358600000414
当
Figure BDA00029189358600000415
为InSAR视线向变形观测值时，
Figure BDA00029189358600000416
当
Figure BDA0002918935860000051
为InSAR方位向形变观测值时，
Figure BDA0002918935860000052
根据所述地表应力应变模型和卫星方位角和入射角，在忽略预设范围内高相干点之间的成像几何差异的前提下(即
Figure BDA0002918935860000053
)，可得
Figure BDA0002918935860000054
Figure BDA0002918935860000055
所述
Figure BDA0002918935860000056
进一步的，所述步骤S4具体包括：
获取观测向量ΔLts的初始权矩阵
Figure BDA0002918935860000057
令
Figure BDA0002918935860000058
Figure BDA0002918935860000059
基于最小二乘准则和矩阵伪逆得到未知参数向量lts的平差值
Figure BDA00029189358600000510
Figure BDA00029189358600000511
Figure BDA00029189358600000512
代表矩阵
Figure BDA00029189358600000513
的伪逆矩阵；
根据迭代加权最小二乘方法对初始权矩阵进行优化更新权矩阵，获得平差值
Figure BDA00029189358600000514
根据目标弧段在东西向、南北向和垂直向的坐标增量Δxe，Δxn，Δxu，获得目标弧段对应的累积三维地表形变：
Figure BDA00029189358600000515
进一步的，所述根据迭代加权最小二乘方法对初始权矩阵进行优化和更新的步骤具体包括：
根据平差值
Figure BDA00029189358600000516
可得观测值向量ΔLts的残差
Figure BDA00029189358600000517
根据观测值残差得到更新后的权矩阵
Figure BDA00029189358600000518
°代表哈达玛运算，G为对角矩阵，G的第z个对角线元素gz可基于以下权函数确定：
Figure BDA0002918935860000061
式中z＝1，2，…，Z,Z代表观测值向量ΔLts中元素的个数，uz代表标准化残差，Δrz代表Δrts中的第z个元素，Δrm为绝对值Δrts的中位值，hz为方阵Bts((Bts)TBts)+(Bts)T的第z个对角线元素，c是一个校准常数；
利用
Figure BDA0002918935860000062
代替
Figure BDA0002918935860000063
重新计算未知参数lts的平差值
Figure BDA0002918935860000064
迭代至
Figure BDA0002918935860000065
Figure BDA0002918935860000066
∈表示一个事先设定的阈值。
进一步的，所述步骤S5具体包括：
获取高相干点上的相位与弧段上的相位差之间的系数矩阵Bp2a；所有目标弧段每个时刻的观测向量Δd以及对应Δd的权矩阵WΔd，基于加权最小二乘准则获得高相干点时序三维地表形变d：
d＝((Bp2a)TWΔdBp2a)-1(Bp2a)TWΔdΔd。
进一步的，所述预设范围为目标弧段中心周围1km×1km范围。
有益效果：
本发明基于地表应力应变模型建立了邻近点InSAR相位差与中心弧段处时序三维地表形变梯度之间的函数模型；然后利用迭代加权最小二乘方法进行时序三维地表形变梯度的求解；根据形变梯度及目标弧段的坐标增量即可得到目标弧段上的时序三维地表形变；最后，对所有弧段上的形变结果进行空间积分，即得到每一个高相干点处的时序三维地表形变，该方法通过求解形变梯度的方式，有效避免了InSAR干涉图的解缠步骤；通过迭代加权最小二乘方法求解，降低了粗差及误差较大的观测值对未知参数解算的影响。整个方法流程均是基于最小二乘准则进行解算，操作简单，通俗易懂，进一步丰富了InSAR三维地表形变监测方法体系。
应当理解的是，以上的一般描述和后文的细节描述仅是示例性和解释性的，并不能限制本公开。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1为本发明实施例提供的一种面向缠绕相位的InSAR时序三维形变监测方法的流程图；
图2为本发明实施例提供的面向缠绕相位的InSAR时序三维形变监测方法模拟实验所用的降轨右视缠绕InSAR干涉图；
图3为本发明实施例提供的面向缠绕相位的InSAR时序三维形变监测方法模拟实验不同方法得到的三维时序形变对比图。
具体实施方式
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。
如图1所示，在本发明实施例中，提出了面向缠绕相位的InSAR时序三维形变监测方法，具体包括以下步骤：
步骤S1，获取多轨道时序SAR影像，将同一轨道SAR影像进行配准与差分干涉，获得预设时空基线阈值的多轨道InSAR干涉对数据集。
在本发明实施例中，获取得到特定研究区域的多轨道时序SAR影像，比如获取J个轨道成像几何差异明显的时序SAR数据集，基于第j(j＝1，2，..，J)个轨道Sj个SAR影像，首先将Sj个SAR影像与该轨道获取的第一个时刻的SAR影像进行配准，然后基于预设的时空基线阈值，对配准后的Sj个SAR影像进行差分干涉处理，共生成了Ij个InSAR干涉对，所有轨道共包含
Figure BDA0002918935860000081
Figure BDA0002918935860000082
个SAR影像和
Figure BDA0002918935860000083
个InSAR干涉对，时空基线阈值一般根据研究区域地表稳定性、SAR卫星波长等因素凭经验确定，一般情况下，时空基线阈值可分别设置为100天和200米。
步骤S2，采用振幅离差法获取多轨道SAR影像中的时序高相干点数据集；并基于所述高相干点构建狄洛尼三角网，经预处理后获得多个目标弧段。
在本发明实施例中，振幅离差法是经典永久散射体InSAR(Permanent ScatterersInSAR，PSInSAR)技术筛选高相干点的方法，该方法首先计算某一像素点所对应的时序SAR影像后向散射强度的标准差σA与平均值mA的比值
Figure BDA0002918935860000084
然后将DA<0.25对应的像素点作为高相干点。本发明中，为了获取更为密集的高相干候选点，取DA<0.4作为振幅离差法选取高相干点数据集的条件。
在本发明实施例中，根据获得的不同轨道SAR数据选取的高相干点整合成一个数据集，共包含
Figure BDA0002918935860000085
个高相干点，进而构建出1个完整的狄洛尼三角网，根据构建的狄洛尼三角网和高相干点集，获取了网中每一个弧段中点的坐标和所有高相干点的边界多边形，进而将中点位于边界多边形外的对应弧段剔除，网中剩余的有效弧段数为
Figure BDA0002918935860000086
个。
步骤S3，基于地表应力应变模型，建立目标弧段处预设范围内所有弧段上的InSAR相位差与目标弧段上时序三维地表变形梯度的函数关系。
在本发明实施例中，所述地表应力应变模型描述的是地表邻近点三维地表形变之间的物理力学关系，假如地表邻近两点P0，Pk的三维坐标和三维形变分别为
Figure BDA0002918935860000091
和
Figure BDA0002918935860000092
则基于地表应力应变模型可得：
dk＝H·Δxk+d0 (1)
其中
Figure BDA0002918935860000093
e，n，u分别代表东西向、南北向和垂直向。H表示点P0处三维形变梯度矩阵，可表示为
Figure BDA0002918935860000094
其中
Figure BDA0002918935860000095
代表求偏导，
Figure BDA0002918935860000096
表示dir1方向的地表形变在dir2方向的形变梯度(dir1和dir2代表e，n，u)，例如
Figure BDA0002918935860000097
代表东西向形变在南北方向的形变梯度。
在本发明实施例中，在P0，Pk处的第j个轨道所对应的观测数据分别为
Figure BDA0002918935860000098
所对应的卫星方位角和入射角分别为
Figure BDA0002918935860000099
和
Figure BDA00029189358600000910
Figure BDA00029189358600000911
和
Figure BDA00029189358600000912
则根据SAR卫星成像几何可得：
Figure BDA00029189358600000913
Figure BDA00029189358600000914
当
Figure BDA00029189358600000915
为InSAR视线向变形观测值时，
Figure BDA00029189358600000916
当
Figure BDA0002918935860000101
为InSAR方位向形变观测值时，
Figure BDA0002918935860000102
由于P0，Pk两点距离较近，则可忽略两点间方位角和入射角的差异，即可认为
Figure BDA0002918935860000103
联合式(1)-(3)可得：
Figure BDA0002918935860000104
Figure BDA0002918935860000105
Figure BDA0002918935860000106
式(5)中的
Figure BDA0002918935860000107
即为第j个轨道InSAR观测值在P0，Pk两点组成的弧段上的相位差。
取目标弧段中心周围1km×1km范围内的InSAR观测值与目标弧段处的时序三维地表形变梯度建立函数关系。假设1km×1km范围内的第j个轨道的高相干点共有Kj+1个，设目标弧段的起点为P0，则基于剩余的Kj个Pk，即可得到在单个InSAR干涉图中，一定范围内弧段上的InSAR相位差与目标弧段处三维形变梯度之间的函数关系：
ΔLj＝Bj,sm·l (8)
式中
Figure BDA0002918935860000108
对于短时空基线InSAR干涉对中，空间上邻近两点的真实形变相位差可认为是介于[-π，π]之间，由于原始InSAR观测相位是缠绕的，因此ΔLj中可能包含2π或-2π的整周模糊度，因此对ΔLj进行如下操作之后才可以用于后续计算：
ΔLj＝mod(ΔLj+π，2π)-π (9)
mod(p1，p2)代表p1除以p2得到的余数。
根据步骤S1中第j个轨道生成的Ij个InSAR干涉对与所有S个SAR影像的构网关系，可得干涉对相位与时序相位之间的系数矩阵为
Figure BDA0002918935860000111
Bsbl，j矩阵的大小为Ij×(S-1)；Bsbl，j的每一列代表所有轨道SAR影像按时间排序后所对应的SAR影像，每一行对应一个InSAR干涉对，例如第一个干涉对是由第2个和第4个SAR影像组成，则Bsbl，，j矩阵中第一行的第2列为-1，第4列为1，其余元素均为0。
联合式(8)和(10)，在考虑InSAR时序数据的基础上，点P0的1km×1km范围内的InSAR观测值与P0处InSAR观测值之差ΔLj，ts和时序三维地表形变梯度lts之间的函数关系为：
ΔLj，ts＝Bj，ts·lts (11)
其中
Figure BDA0002918935860000112
Figure BDA0002918935860000113
代表第j个轨道中第u个InSAR干涉对所对应的Kj，i个Pk与目标弧段起点P0处的InSAR观测相位之差
Figure BDA0002918935860000114
Figure BDA0002918935860000115
Figure BDA0002918935860000116
代表系数矩阵，
Figure BDA0002918935860000117
为克罗内克积符号，
Figure BDA0002918935860000118
代表时序三维地表形变梯度，
Figure BDA0002918935860000119
为第s个SAR影像获取时刻对应的累积三维地表形变梯度，下标ts代表时序。
考虑到共有J个轨道的InSAR数据，则根据式(11)可得，在多个轨道InSAR时序干涉图中，一定范围内弧段上的InSAR相位差与时序三维形变梯度之间的函数关系
ΔLts＝Bts·lts (12)
式中ΔLts＝[(ΔL1，ts)T，(ΔL2，ts)T，…，(ΔLj，ts)T，…，(ΔLJ，ts)T]T，Bts＝[(B1，ts)T，(B2，ts)T，…，(Bj，ts)T，…，(BJ，ts)T]T。
步骤S4，采用迭代加权最小二乘方法计算目标弧段时序三维地表变形梯度，并根据所述目标弧段所对应的坐标增量，获得目标弧段所对应的三维地表形变。
在本发明实施例中，求解式(12)中的未知参数lts之前，需事先确定观测向量ΔLts的初始权矩阵
Figure BDA0002918935860000121
所述
Figure BDA0002918935860000122
为对角矩阵，对角线元素
Figure BDA0002918935860000123
代表相应观测值的精度，
Figure BDA0002918935860000124
越大代表精度越高，本发明中以点Pk处的相干性作为相应观测值
Figure BDA0002918935860000125
的初始权重
Figure BDA0002918935860000126
确定初始权矩阵
Figure BDA0002918935860000127
之后，可令
Figure BDA0002918935860000128
Figure BDA0002918935860000129
本发明则基于最小二乘准则和矩阵伪逆得到未知参数向量lts的平差值
Figure BDA00029189358600001210
Figure BDA00029189358600001211
式中
Figure BDA00029189358600001212
代表矩阵
Figure BDA00029189358600001213
的伪逆矩阵。由于初始权矩阵
Figure BDA00029189358600001214
往往难以准确表征观测值向量的精度水平，且观测值向量ΔLts中可能包含2π模糊度，导致未知参数向量的求解不是最优的。本发明则引入迭代加权最小二乘方法，通过不断优化权矩阵数值，以提高未知参数的求解精度，具体包括：
根据平差值
Figure BDA00029189358600001215
可得观测值向量ΔLts的残差Δrts为：
Figure BDA00029189358600001216
可根据观测值残差得到更新后的权矩阵
Figure BDA00029189358600001217
Figure BDA00029189358600001218
式中°代表哈达玛运算，G为对角矩阵，其对角线元素可认为是降权因子。其中，G的第z个对角线元素gz可基于以下权函数确定
Figure BDA00029189358600001219
Figure BDA0002918935860000131
式中z＝1，2，…，Z,Z代表观测值向量ΔLts中元素的个数，uz代表标准化残差，Δrz代表Δrts中的第z个元素，Δrm为绝对值Δrts的中位值，hz为方阵Bts((Bts)TBts)+(Bts)T的第z个对角线元素，c是一个校准常数，在95％的置信度下可取值为4.685.利用
Figure BDA0002918935860000132
代替
Figure BDA0002918935860000133
重新计算未知参数lts的平差值
Figure BDA0002918935860000134
迭代式(13)-(15)直至
Figure BDA0002918935860000135
∈表示一个事先设定的阈值。
根据式(13)得到的未知参数向量的平差值，即可求解目标弧段上的时序三维地表形变；以第s个SAR影像获取时刻为例，弧段上对应的累积三维地表形变[Δde，s，Δdn，s，Δdu，s]T为：
Figure BDA0002918935860000136
式中Δxe，Δxn，Δxu表示目标弧段在东西向、南北向和垂直向的坐标增量。
根据观测值残差Δrts可得到目标弧段的时序相干性γ：
Figure BDA0002918935860000137
目标弧段上的时序相干性γ可在后续对弧段进行空间积分的时候作为权重因子，以降低噪声对空间积分的影响。
重复上述步骤3、4，直至得到所有目标弧段所对应的时序三维地表形变。
步骤S5，获取所有目标弧段的三维地表形变，对每个时刻所有弧段上的三维地表形变依次进行空间积分，获得所有高相干点上的时序三维地表形变结果。
在本发明实施例中，选定某一稳定点为参考，对每个时刻所有弧段上的三维地表形变依次进行空间积分，即可得到剩余所有
Figure BDA0002918935860000147
个高相干点上的时序三维地表形变结果，具体为：
高相干点上的相位与弧段上的相位差之间的系数矩阵Bp2a可表示为
Figure BDA0002918935860000141
其中，矩阵Bp2a的大小为
Figure BDA0002918935860000142
每一行代表一个弧段，每一列代表一个高相干点，其中，选取的参考点不包含在内。例如第一个弧段是由参考点和第1个高相干点组成，则Bp2a矩阵中第一行的第1列为1。第2个弧段是由第1个和第2个高相干点组成，则Bp2a矩阵中第2行的第1列为-1，第2列为1，其余元素均为0，以此类推，即可得到矩阵Bp2a。
以某一时刻的东西向形变为例，所有弧段上在该时刻所对应的东西向形变组成的观测向量为
Figure BDA0002918935860000143
Figure BDA0002918935860000144
代表第a个弧段上的东西向形变。对应Δde的权矩阵WΔd为对角矩阵，矩阵大小为
Figure BDA0002918935860000145
第a个对角线元素即为第a个弧度对应的时序相干性γa。进而基于加权最小二乘准则即可得到
Figure BDA0002918935860000146
个高相干点上对应的东西向形变de：
de＝((Bp2a)TWΔdBp2a)-1(Bp2a)TWΔdΔde (22)
其他时刻和其他维度进行空间积分的方法按照式(22)的方法进行空间积分，获得所有高相干点上的时序三维地表形变结果。
本发明通过获取多轨道时序SAR影像，分别经配准与差分干涉和振幅离差法获取多轨道InSAR干涉对数据集和高相干点数据集，基于高相干点数据集构建狄洛尼三角网，经预处理获得有效弧段，基于地表应力应变模型建立了邻近点InSAR相位差与目标弧段处时序三维地表形变梯度之间的函数模型；然后利用迭代加权最小二乘方法进行时序三维地表形变梯度的求解；根据形变梯度及目标弧段的坐标增量即可得到目标弧段上的时序三维地表形变，最后对所有弧段上的形变结果进行空间积分，获得高相干点处的时序三维地表形变，该方法通过求解形变梯度的方式，有效避免了InSAR干涉图的解缠步骤；通过迭代加权最小二乘方法求解，降低了粗差及误差较大的观测值对未知参数解算的影响。整个方法流程均是基于最小二乘准则进行解算，操作简单，获得了一种高效、准确的InSAR三维地表形变监测方法。
以下通过模拟实验进行进一步说明：
①根据现有SAR卫星的特征，模拟得到升降轨左右视四个轨道数据的成像几何参数及时空基线配置；
②在一定区域(图像大小200×200)模拟生成因地下流体变化导致的三维地表形变数据，并以此数据作为时序三维地表形变的空间特征，假定时序线性形变，模拟生成对应4个轨道(升轨左视、升轨右视、降轨左视和降轨右视)的InSAR干涉相位；
③假设干涉相位的相干性在[0.7,0.9]之间均匀分布，则可模拟生成失相干噪声，并加到干涉相位上，同时将加有噪声的干涉相位进行缠绕，得到模拟数据所用的InSAR缠绕相位观测值。如图2所示，为降轨右视InSAR缠绕干涉图。
本发明将步骤③中缠绕之前的InSAR干涉相位作为观测值，利用多维SBAS方法进行时序三维地表形变解算；如图3所示，为图2干涉图20180101_20180429中白色三角形处，模拟数据、本发明方法和传统多维SBAS方法得到的时序三维形变对比图，可以看出，本文发明方法在InSAR缠绕相位的情况下，可以准确地得到时序三维形变。同时，由于考虑了邻近点形变之间的物理力学关系，本发明方法比传统多维SBAS方法得到的三维形变结果更为精确。
以上所述实施例仅表达了本发明的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对本发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。因此，本发明专利的保护范围应以所附权利要求为准。
本领域技术人员在考虑说明书及实践这里公开的发明后，将容易想到本公开的其它实施方案。本申请旨在涵盖本公开的任何变型、用途或者适应性变化，这些变型、用途或者适应性变化遵循本公开的一般性原理并包括本公开未公开的本技术领域中的公知常识或惯用技术手段。说明书和实施例仅被视为示例性的，本公开的真正范围和精神由权利要求指出。
应该理解的是，虽然本发明各实施例的流程图中的各个步骤按照箭头的指示依次显示，但是这些步骤并不是必然按照箭头指示的顺序依次执行。除非本文中有明确的说明，这些步骤的执行并没有严格的顺序限制，这些步骤可以以其它的顺序执行。而且，各实施例中的至少一部分步骤可以包括多个子步骤或者多个阶段，这些子步骤或者阶段并不必然是在同一时刻执行完成，而是可以在不同的时刻执行，这些子步骤或者阶段的执行顺序也不必然是依次进行，而是可以与其它步骤或者其它步骤的子步骤或者阶段的至少一部分轮流或者交替地执行。
本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流程，是可以通过计算机程序来指令相关的硬件来完成，所述的程序可存储于一非易失性计算机可读取存储介质中，该程序在执行时，可包括如上述各方法的实施例的流程。其中，本申请所提供的各实施例中所使用的对存储器、存储、数据库或其它介质的任何引用，均可包括非易失性和/或易失性存储器。非易失性存储器可包括只读存储器(ROM)、可编程ROM(PROM)、电可编程ROM(EPROM)、电可擦除可编程ROM(EEPROM)或闪存。易失性存储器可包括随机存取存储器(RAM)或者外部高速缓冲存储器。作为说明而非局限，RAM以多种形式可得，诸如静态RAM(SRAM)、动态RAM(DRAM)、同步DRAM(SDRAM)、双数据率SDRAM(DDRSDRAM)、增强型SDRAM(ESDRAM)、同步链路(Synchlink)DRAM(SLDRAM)、存储器总线(Rambus)直接RAM(RDRAM)、直接存储器总线动态RAM(DRDRAM)、以及存储器总线动态RAM(RDRAM)等。
以上所述实施例的各技术特征可以进行任意的组合，为使描述简洁，未对上述实施例中的各个技术特征所有可能的组合都进行描述，然而，只要这些技术特征的组合不存在矛盾，都应当认为是本说明书记载的范围。

Claims (10)

1.一种面向缠绕相位的InSAR时序三维形变监测方法，其特征在于，具体包括：
S1：获取多轨道时序SAR影像，将同一轨道SAR影像进行配准与差分干涉，获得预设时空基线阈值的多轨道InSAR干涉对数据集；
S2：采用振幅离差法获取多轨道SAR影像中的时序高相干点数据集；并基于所述高相干点构建狄洛尼三角网，经预处理后获得多个目标弧段；
S3：基于地表应力应变模型，建立目标弧段处预设范围内所有弧段上的InSAR相位差与目标弧段上时序三维地表变形梯度的函数关系；
S4：采用迭代加权最小二乘方法计算目标弧段时序三维地表变形梯度，并根据所述目标弧段所对应的坐标增量，获得目标弧段所对应的三维地表形变；
S5：获取所有目标弧段的三维地表形变，对每个时刻所有弧段上的三维地表形变依次进行空间积分，获得所有高相干点上的时序三维地表形变结果。
2.根据权利要求1所述的面向缠绕相位的InSAR时序三维形变监测方法，其特征在于，所述预设时空基线阈值分别为100天和200m。
3.根据权利要求1所述的面向缠绕相位的InSAR时序三维形变监测方法，其特征在于，所述步骤S2中的采用振幅离差法获取多轨道SAR影像中的时序高相干点数据集步骤具体包括：
计算所述SAR影像像素点的后向散射强度的标准差σA与平均值mA的比值
Figure FDA0002918935850000011
选取DA＜0.4的像素点作为高相干点，获取高相干点数据集。
4.根据根据权利要求1所述的面向缠绕相位的InSAR时序三维形变监测方法，其特征在于，所述步骤S2中的预处理具体为：
基于构建的狄洛尼三角网和高相干点集，获取所述狄洛尼三角网中每一个弧段中点的坐标和所有高相干点的边界多边形；
将中点位于边界多边形外的对应弧段剔除，剩余的有效弧段即为目标弧段。
5.根据根据权利要求1所述的面向缠绕相位的InSAR时序三维形变监测方法，其特征在于，所述步骤S3具体包括：
将目标弧段预设范围内的高相干点处的InSAR观测相位与目标弧段的起点处的InSAR观测相位进行做差，获得目标弧段预设范围内的高相干点与目标弧段起点的InSAR相位差；
基于地表应力应变模型和卫星的方位角及入射角，获得单个InSAR干涉图中，目标弧段预设范围内的InSAR相位差与目标弧段处三维形变梯度之间的函数关系Bj，sm；
根据所述干涉对数据集与SAR影像的构网关系，获取干涉对相位与时序相位之间的系数矩阵Bsbl，j；
根据所述目标弧段预设范围内的高相干点与目标弧段起点的InSAR相位差、单个干涉图预设范围内弧段上的InSAR相位差与目标弧段处三维形变梯度之间的函数关系Bj，sm、以及干涉对相位与时序相位之间的系数矩阵Bsbl，j，构建在多个轨道InSAR时序干涉图中，预设范围内弧段上的InSAR相位差与时序三维形变梯度之间的函数关系为：
ΔLts＝Bts·lts
所述ΔLts＝[(ΔL1，ts)T，(ΔL2，ts)T，...，(ΔLj，ts)T，...，(ΔLJ，ts)T]T，Bts＝[(B1，ts)T，(B2，ts)T，...，(Bj，ts)T，...，(BJ，ts)T]T，
Figure FDA0002918935850000021
代表第j个轨道中第i个InSAR干涉对所对应的Kj，i个Pk与目标弧段起点P0组成的弧段上的InSAR相位差
Figure FDA0002918935850000022
Figure FDA0002918935850000023
代表系数矩阵，
Figure FDA0002918935850000024
为克罗内克积符号，
Figure FDA0002918935850000031
代表时序三维地表形变梯度，
Figure FDA0002918935850000032
Figure FDA0002918935850000033
为第s个SAR影像获取时刻对应的累积三维地表形变梯度，下标ts代表时序。
6.根据权利要求5所述的面向缠绕相位的InSAR时序三维形变监测方法，其特征在于，所述基于地表应力应变模型和卫星的方位角及入射角，获得单个InSAR干涉图中，预设范围内弧段上的InSAR相位差与目标弧段处三维形变梯度之间的函数关系Bj，sm，具体包括：
所述地表应力应变模型为dk＝H·Δxk+d0，所述dk、d0为预设范围内高相干点处和目标弧段起点的三维形变，H为目标弧段起点的三维形变梯度矩阵，Δxk＝xk-x0为预设范围内高相干点处和目标弧段起点的三维坐标之差；
所述预设范围内高相干点处和目标弧段起点的InSAR观测数据分别
Figure FDA0002918935850000034
所对应的卫星方位角和入射角分别为
Figure FDA0002918935850000035
和
Figure FDA0002918935850000036
和
Figure FDA0002918935850000037
根据SAR卫星成像几何可得函数关系：
Figure FDA0002918935850000038
式中，
Figure FDA0002918935850000039
当
Figure FDA00029189358500000310
为InSAR视线向形变观测值时，
Figure FDA00029189358500000311
当
Figure FDA00029189358500000312
为InSAR方位向形变观测值时，
Figure FDA00029189358500000313
根据所述地表应力应变模型和卫星方位角和入射角，在忽略预设范围内高相干点之间的成像几何差异的前提下(即
Figure FDA0002918935850000041
)，可得
Figure FDA0002918935850000042
Figure FDA0002918935850000043
所述
Figure FDA0002918935850000044
7.根据权利要求5所述的面向缠绕相位的InSAR时序三维形变监测方法，其特征在于，所述步骤S4具体包括：
获取观测向量ΔLts的初始权矩阵
Figure FDA0002918935850000045
令
Figure FDA0002918935850000046
Figure FDA0002918935850000047
基于最小二乘准则和矩阵伪逆得到未知参数向量lts的平差值
Figure FDA0002918935850000048
Figure FDA0002918935850000049
Figure FDA00029189358500000410
代表矩阵
Figure FDA00029189358500000411
的伪逆矩阵；
根据迭代加权最小二乘方法对初始权矩阵进行优化更新权矩阵，获得平差值
Figure FDA00029189358500000412
根据目标弧段在东西向、南北向和垂直向的坐标增量Δxe，Δxn，Δxu，获得目标弧段对应的累积三维地表形变：
Figure FDA00029189358500000413
8.根据权利要求7所述面向缠绕相位的InSAR时序三维形变监测方法，其特征在于，所述根据迭代加权最小二乘方法对初始权矩阵进行优化和更新的步骤具体包括：
根据平差值
Figure FDA00029189358500000414
可得观测值向量ΔLts的残差
Figure FDA00029189358500000415
根据观测值残差得到更新后的权矩阵
Figure FDA00029189358500000416
代表哈达玛运算，G为对角矩阵，G的第z个对角线元素gz可基于以下权函数确定：
Figure FDA00029189358500000417
式中z＝1，2，…，Z，Z代表观测值向量ΔLts中元素的个数，uz代表标准化残差，Δrz代表Δrts中的第z个元素，Δrm为绝对值Δrts的中位值，hz为方阵Bts((Bts)TBts)+(Bts)T的第z个对角线元素，c是一个校准常数；
利用
Figure FDA0002918935850000051
代替
Figure FDA0002918935850000052
重新计算未知参数lts的平差值
Figure FDA0002918935850000053
迭代至
Figure FDA0002918935850000054
∈表示一个事先设定的阈值。
9.根据权利要求1所述的面向缠绕相位的InSAR时序三维形变监测方法，其特征在于，所述步骤S5具体包括：
获取高相干点上的相位与弧段上的相位差之间的系数矩阵Bp2a；所有目标弧段每个时刻的观测向量Δd以及对应Δd的权矩阵WΔd，基于加权最小二乘准则获得高相干点时序三维地表形变d：
d＝((Bp2a)TWΔdBp2a)-1(Bp2a)TWΔdΔd。
10.根据权利要求1所述的面向缠绕相位的InSAR时序三维形变监测方法，其特征在于，所述预设范围为目标弧段中心周围1km×1km范围。
CN202110111179.8A 2021-01-27 2021-01-27 面向缠绕相位的InSAR时序三维形变监测方法 Active CN112797886B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110111179.8A CN112797886B (zh) 2021-01-27 2021-01-27 面向缠绕相位的InSAR时序三维形变监测方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110111179.8A CN112797886B (zh) 2021-01-27 2021-01-27 面向缠绕相位的InSAR时序三维形变监测方法

Publications (2)

Publication Number Publication Date
CN112797886A true CN112797886A (zh) 2021-05-14
CN112797886B CN112797886B (zh) 2022-04-22

Family

ID=75812100

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110111179.8A Active CN112797886B (zh) 2021-01-27 2021-01-27 面向缠绕相位的InSAR时序三维形变监测方法

Country Status (1)

Country Link
CN (1) CN112797886B (zh)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113567979A (zh) * 2021-06-03 2021-10-29 长安大学 基于模拟退火算法的多时相InSAR相位解缠绕方法
CN115993601A (zh) * 2023-03-22 2023-04-21 四川省公路规划勘察设计研究院有限公司 一种强盐渍土区域公路变形的时序InSAR监测方法
WO2023142205A1 (zh) * 2022-01-26 2023-08-03 中山大学 一种InSAR时序相位的优化方法及装置

Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5659318A (en) * 1996-05-31 1997-08-19 California Institute Of Technology Interferometric SAR processor for elevation
CN103698749A (zh) * 2013-12-31 2014-04-02 中国人民解放军国防科学技术大学 一种利用小数据集sar图像序列提取永久散射体的方法
CN104062660A (zh) * 2014-07-14 2014-09-24 中南大学 一种基于时域离散InSAR干涉对的矿区地表时序形变监测方法
CN104123470A (zh) * 2014-07-25 2014-10-29 首都师范大学 一种优化地面沉降监测网的方法
CN109738892A (zh) * 2019-01-24 2019-05-10 中南大学 一种矿区地表高时空分辨率三维形变估计方法
CN110888130A (zh) * 2019-10-30 2020-03-17 华东师范大学 一种基于升降轨时序InSAR的煤矿区地表形变监测方法
CN111398959A (zh) * 2020-04-07 2020-07-10 中南大学 基于地表应力应变模型的InSAR时序地表形变监测方法
CN111998766A (zh) * 2020-08-31 2020-11-27 同济大学 一种基于时序InSAR技术的地表形变反演方法
US20210011149A1 (en) * 2019-05-21 2021-01-14 Central South University InSAR and GNSS weighting method for three-dimensional surface deformation estimation

Patent Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5659318A (en) * 1996-05-31 1997-08-19 California Institute Of Technology Interferometric SAR processor for elevation
CN103698749A (zh) * 2013-12-31 2014-04-02 中国人民解放军国防科学技术大学 一种利用小数据集sar图像序列提取永久散射体的方法
CN104062660A (zh) * 2014-07-14 2014-09-24 中南大学 一种基于时域离散InSAR干涉对的矿区地表时序形变监测方法
CN104123470A (zh) * 2014-07-25 2014-10-29 首都师范大学 一种优化地面沉降监测网的方法
CN109738892A (zh) * 2019-01-24 2019-05-10 中南大学 一种矿区地表高时空分辨率三维形变估计方法
US20210011149A1 (en) * 2019-05-21 2021-01-14 Central South University InSAR and GNSS weighting method for three-dimensional surface deformation estimation
CN110888130A (zh) * 2019-10-30 2020-03-17 华东师范大学 一种基于升降轨时序InSAR的煤矿区地表形变监测方法
CN111398959A (zh) * 2020-04-07 2020-07-10 中南大学 基于地表应力应变模型的InSAR时序地表形变监测方法
CN111998766A (zh) * 2020-08-31 2020-11-27 同济大学 一种基于时序InSAR技术的地表形变反演方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
YANG, ZEFA; LI, ZHIWEI; ZHU, JIANJUN; 等: "Deriving time-series three-dimensional displacements of mining areas from a single-geometry InSAR dataset", 《JOURNAL OF GEODESY》 *
王心雨: "高分辨率数据用于西安市地面沉降的InSAR监测研究", 《中国优秀硕士学位论文全文数据库 (基础科学辑)》 *

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113567979A (zh) * 2021-06-03 2021-10-29 长安大学 基于模拟退火算法的多时相InSAR相位解缠绕方法
CN113567979B (zh) * 2021-06-03 2023-08-29 长安大学 基于模拟退火算法的多时相InSAR相位解缠绕方法
WO2023142205A1 (zh) * 2022-01-26 2023-08-03 中山大学 一种InSAR时序相位的优化方法及装置
CN115993601A (zh) * 2023-03-22 2023-04-21 四川省公路规划勘察设计研究院有限公司 一种强盐渍土区域公路变形的时序InSAR监测方法

Also Published As

Publication number Publication date
CN112797886B (zh) 2022-04-22

Similar Documents

Publication Publication Date Title
CN112797886B (zh) 面向缠绕相位的InSAR时序三维形变监测方法
Yunjun et al. Small baseline InSAR time series analysis: Unwrapping error correction and noise reduction
Yague-Martinez et al. Coregistration of interferometric stacks of Sentinel-1 TOPS data
Zhang et al. Modeling PSInSAR time series without phase unwrapping
Xu et al. A refined strategy for removing composite errors of SAR interferogram
CN107102333B (zh) 一种星载InSAR长短基线融合解缠方法
Fornaro et al. A null-space method for the phase unwrapping of multitemporal SAR interferometric stacks
Yang et al. Deriving time-series three-dimensional displacements of mining areas from a single-geometry InSAR dataset
WO2010000870A1 (en) Identification and analysis of persistent scatterers in series of sar images
CN111398959A (zh) 基于地表应力应变模型的InSAR时序地表形变监测方法
Xu et al. Time-series InSAR dynamic analysis with robust sequential adjustment
Zhang et al. Minimizing height effects in MTInSAR for deformation detection over built areas
Mirzaee et al. Non-linear phase linking using joined distributed and persistent scatterers
CN112363166B (zh) 一种基于可靠冗余网络的InSAR相位解缠方法和系统
Hellwich et al. Geocoding SAR interferograms by least squares adjustment
Zhang et al. Deformation rate estimation on changing landscapes using temporarily coherent point InSAR
Sosnovsky et al. A Method of Phase Unwrapping Algorithms Efficiency Analysis for InSAR Data Processing
Ajourlou et al. A new strategy for phase unwrapping in insar time series over areas with high deformation rate: Case study on the southern tehran subsidence
Costantini et al. SAR interferometric baseline calibration without need of phase unwrapping
Biswas et al. Investigation of ground deformation with PSInSAR approach in an unstable urban area Naples, Italy using X-band SAR images
CN110286374B (zh) 基于分形布朗运动的干涉sar影像仿真方法
Huang et al. Reference network construction for persistent scatterer detection in SAR tomography: Ant colony search algorithm (ACSA)
Capková et al. Detecting land deformation in the area of northern bohemia using insar stacks (preliminary results)
Dawson Satellite radar interferometry with application to the observation of surface deformation in Australia
Pepe et al. DEM correction and mean surface displacement rate retrieval from a stack of wrapped multi-temporal DInSAR interferograms

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