一种基于二维观测的地表三维变形insar测量方法2t3d-insar
技术领域
1.本发明涉及卫星遥感测量地表变形的领域，具体涉及一种基于二维观测的地表三维变形insar监测方法。
背景技术：2.合成孔径雷达差分干涉测量（dinsar，differential interferometry synthetic aperture radar，也泛称insar）是一种高精度对地微波观测技术，利用同一地区的两幅具有相干性的sar（synthetic aperture radar）影像的重复观测形成干涉变形，根据图像像元的相位信息获取地表形变信息（图1），其测距原理类似高中物理课双缝干涉实验，利用条纹特性来测量缝与板间的距离。insar技术在测量地面高程和地表形变方面具有范围广、精度高、全覆盖、可回溯、全天时、全天候和非接触式测量的特点，广泛应用于地质灾害、地震和火山形变、冰川运动监测、城市沉降监测、矿山形变监测和大型建筑的形变监测等。
3.但是，insar（干涉合成孔径雷达）技术只能测量卫星与地表目标点间相对距离的变化（图1），即卫星视线向的变形，无法像gnss、全站仪等传统地面测量方法一样获得空间的x、y、z三维变形。
4.为了解决这个问题，前人通过三个方向的insar观测分解来求取地表x、y、z三维变形，但问题是sar卫星只能从升轨、降轨两个水平方位进行insar观测（图2），由于sar卫星属于太阳同步卫星，各种卫星之间观测的水平角度差距不大（一般
±5°
），观测的变形方向同质化严重，空间几何差异小，引入第三种insar观测分解出来的三维变形误差较大，甚至无法分解。
技术实现要素：5.基于上述技术问题，本发明提出一种基于二维观测的地表三维变形insar监测方法—2t3d-insar，其仅用升轨、降轨2个方向的观测，结合地形条件，求取地表三维变形。
6.本发明所采用的技术解决方案如下：一种基于二维观测的地表三维变形insar监测方法，包括以下步骤：s1：对监测目标进行insar升、降轨两个方向的观测，准备原始观测数据和参数数据，分别采用升轨和降轨各至少两期sar数据，应用dinsar、stacking-insar、sbas-insar或ps-insar技术获取监测对象的雷达视线向变形量或变形速率，升、降轨观测值分别为v
as
和v
des
，并在sar数据的参数文件中查找雷达入射角θ和轨道方位角α值。
7.s2：利用dem计算观测区每个像素点的坡度β和坡向λ，作为后续的地表联动变形假设模型的参数。
8.s3：建立insar实测变形与理论模型相结合的二维至三维空间转化的公式，实际应用中习惯采用垂直、东西、南北三个方向的变形量来表征地表的三维变形，分别设为v
u，
向上为正，v
e，
东为正，v
n，
北为正，根据雷达观测的空间几何特征得到视线向和三维形变量转换
关系为第一公式：θa升轨变形入射角，θd降轨变形入射角，αa升轨轨道方位角，αd降轨轨道方位角；设：insar变形观测向量l=[v
as
，v
des
] t
三维变形向量d= [vu，ve，vn]
t
转化参数向量b =第一公式中由于b不是方阵，不可逆，也就是根据二个雷达视线向变形向量l是无法求出空间三维变形向量d，因此需要其他方法的辅助。
[0009]
地表的岩土体是连续介质，具有空间上的三维应变关联，理论上任何一个方向的连续变形都会关联到其他两个方向，斜坡垂向变形是沿坡度的正弦（sin）分量，与水平分量是正切关系（坡度的tan值，也就是坡比），此时三维变形的关系表示为第二公式：-vu=ve·
se+vn·
snvu垂直变形量，ve东西方向变形量，vn南北方向变形量，se东西方向坡比，sn南北方向坡比，其中se和sn是由坡度（β）和坡向（λ）计算；将第一公式和第二公式结合，可建立二维insar变形观测求取三维空间变形的第三公式；通过逆矩阵运算就可求得三维空间变形向量d，涉及公式为第四公式：d=b-1
l。
[0010]
s4：将v
as
和v
des
的形变图像信息，利用dem分析计算得到的β和λ图像，sar原始参数θ和α图像（升、降轨），共8个图重采样到相同分辨率的图像并使各同一地理位置像元对应。然后带入第四公式就可以计算每个像元的垂直vu、东西ve和南北vn三维变形值。
[0011]
特别说明的是本发明中涉及的insar卫星技术和测量参数、数字地面模型dem参数均为本领域技术人员的公知常识，因此不再对此类基本知识进行展开介绍。
[0012]
本发明的有益技术效果是：地表变形是空间三维变形。但是，现有insar技术单一轨道观测仅能观测卫星与地面间距离变化的一维变形；升、降轨观测仅能观测二维变形不足以反算三维变形；再增加观
测方向受轨道限制方向同质化严重，解算的误差大甚至无法解算。本发明是采用2个方位的insar观测，结合滑坡地表变形平行流动关系模型，建立空间变形分解矩阵，然后计算求取滑坡3维变形。与传统需要3个方位insar变形分解比较，仅需升、降轨2个方位的insar观测量，以及利用数字高程模型dem计算观测区每个像素点的坡度（β）和坡向（λ），现今dem数据来源广泛，坡度、坡向计算也成熟。整个发明理论清晰、逻辑严谨、形式整齐，应用时所需数据量小、计算效率高、容错性强，对基于insar观测的地表三维变形求取具有重要用途。
附图说明
[0013]
图1 insar测量地表与卫星间一维变形（视线向）原理示意图；图2 insar升、降轨两个方向测量地表变形方向示意图；图3 为本发明（2t3d-insar）的技术流程图；图4用insar测量的四川金沙江沿岸滑坡案例的升轨、降轨变形图；图5 用于地表流动变形假设公式案例区的坡度和坡向图；图6 利用本发明计算的案例区的三维变形图（箭头代表水平方向，点代表垂直方向）图7 本发明计算的变形值与gnss实测值变形比较图。
具体实施方式
[0014]
现在将参考附图更全面地描述示例实施例。然而，示例实施例能够以多种形式实施，且不应被理解为限于在此阐述的实施例；相反，提供这些实施例使得本技术将全面和完整，并将示例实施例的构思全面地传达给本领域的技术人员。在图中相同的附图标记表示相同或类似的部分，因而将省略对它们的重复描述。
[0015]
所描述的特征、结构或特性可以以任何合适的方式结合在一个或更多实施例中。在下面的描述中，提供许多具体细节从而给出对本公开的实施例的充分理解。
[0016]
附图中所示的流程图仅是示例性说明，不是必须包括所有的内容和操作/步骤，也不是必须按所描述的顺序执行。例如，有的操作/步骤还可以分解，而有的操作/步骤可以合并或部分合并，因此实际执行的顺序有可能根据实际情况改变。
[0017]
本技术的说明书和权利要求书及上述附图中的术语“第一”、“第二”等是用于区别不同对象，而不是用于描述特定顺序。此外，术语“包括”和“具有”以及它们任何变形，意图在于覆盖不排他的包含。例如包含了一系列步骤或单元的过程、方法、参数、结果没有限定于已列出的步骤或单元，而是可选地还包括没有列出的步骤或单元，或可选地还包括对于这些过程、方法、参数、结果固有的其他步骤或单元。
[0018]
s1：采用欧空局哨兵1号卫星的2019.12.04—2020.12.22期间的升轨sar数据，2019.12.06—2020.12.24期间的降轨sar数据，对危岩体进行变形观测，获得升轨和降轨两个方向的变形观测量，在sar数据的参数文件中查找雷达入射角和方位角。图4为升轨和降轨两个方向观测获得的insar变形量图。
[0019]
s2：采用15米分辨率dem在arcgis软件中通过空间分析模块计算获得，案例区的坡度和坡向值。图5中的a图和b图分别为计算得到的坡度和坡向图。
[0020]
s3：直接引用发明中推导的第三公式和第四公式。
[0021]
s4：将哨兵1号卫星中的升、降轨观测图分别提取v
as
和v
des
，通过坡度和坡向图提取坡度角β和坡向角λ，在哨兵1号sar数据的原始参数文件中查找θ和α值，然后带入第三公式，在通过第四公式的矩阵逆计算得到三维变形量，最后在arcgis中进行三维显示，如图6所示。案例设有6个地面监测站，与其比较该发明测量的精度，可达mm级（图7）。
技术特征：1.一种基于二维观测的地表三维变形insar监测方法，其特征在于包括以下步骤：s1：对监测目标进行insar升、降轨两个方向的观测，准备原始观测数据和卫星参数数据；s2：利用dem计算观测区每个像素点的坡度β和坡向λ，作为后续的地表联动变形假设模型的参数；s3：建立insar实测变形与理论模型相结合的二维至三维空间转化的公式；s4：带入insar升降轨观测值、地貌参数、升降轨sar原始参数计算空间三维变形量。2.根据权利要求1所述的一种基于二维观测的地表三维变形insar监测方法，其特征在于，步骤s1中：分别采用升轨和降轨各至少两期sar数据，应用dinsar、stacking-insar、sbas-insar或ps-insar技术获取监测对象的雷达视线向变形量或变形速率，升、降轨观测值分别为v
as
和v
des
，并在sar数据的参数文件中查找雷达入射角θ和轨道方位角α值。3.根据权利要求1所述的一种基于二维观测的滑坡三维变形insar监测方法，其特征在于，步骤s3中：建立insar实测变形与理论模型相结合的二维至三维空间转化的公式，实际应用中习惯采用垂直、东西、南北三个方向的变形量来表征地表的三维变形，分别设为v
u，
向上为正，v
e，
东为正，v
n，
北为正，根据雷达观测的空间几何特征得到视线向和三维形变量转换关系为第一公式：θ
a
升轨变形入射角，θ
d
降轨变形入射角，α
a
升轨轨道方位角，α
d
降轨轨道方位角；设：insar变形观测向量l=[v
as
，v
des
] t
三维变形向量d= [v
u
，v
e
，v
n
]
t
转化参数向量b =第一公式中由于b不是方阵，不可逆，也就是根据二个雷达视线向变形向量l是无法求出空间三维变形向量d，因此需要其他方法的辅助；地表的岩土体是连续介质，具有空间上的三维应变关联，理论上任何一个方向的连续变形都会关联到其他两个方向，斜坡垂向变形是沿坡度的正弦（sin）分量，与水平分量是正切关系（坡度的tan值，也就是坡比），此时三维变形的关系表示为第二公式：-v
u
=v
e
·
s
e
+v
n
·
s
n
。4.v
u
垂直变形量，v
e
东西方向变形量，v
n
南北方向变形量，s
e
东西方向坡比，s
n
南北方向坡比，其中s
e
和s
n
是由坡度（β）和坡向（λ）计算；将第一公式和第二公式结合，可建立二维insar变形观测求取三维空间变形的第三公式；
通过逆矩阵运算就可求得三维空间变形向量d，涉及公式为第四公式：d=b-1
l。5.根据权利要求1所述的一种基于二维观测的滑坡三维变形insar监测方法，其特征在于，步骤s4中：将v
as
和v
des
的形变图像信息，利用dem分析计算得到的β和λ图像，sar原始参数θ和α图像（升、降轨），共8个图重采样到相同分辨率的图像并使各同一地理位置像元对应，然后带入公式4就可以计算每个像元的垂直v
u
、东西v
e
和南北v
n
三维变形值。
技术总结本发明公开了一种基于二维观测的地表三维变形InSAR测量方法——2T3D-InSAR方法，属于卫星遥感测量地表变形的领域。InSAR（Interferematic Synthetic Aperture Radar）技术利用卫星雷达两次回波的相位差计算地表变形。本发明采用2个方向的InSAR观测，结合地表变形平行流假设模型，建立空间变形分解矩阵，然后求取地表3维变形。解决了现有InSAR技术只能测量一个方向（卫星视线向）变形的问题，扩展了InSAR的应用领域；且与传统三维变形分解方法比较，本方法仅需2个方向的InSAR观测量，计算效率高、数据获取容易、容错性强。容错性强。容错性强。
技术研发人员：姚鑫
受保护的技术使用者：姚鑫
技术研发日：2022.01.06
技术公布日：2022/7/5
