CN114236541A - 基于Sentinel-1卫星SAR图像的大面积地表三维形变计算方法 - Google Patents

基于Sentinel-1卫星SAR图像的大面积地表三维形变计算方法 Download PDF

Info

Publication number
CN114236541A
CN114236541A CN202111489543.0A CN202111489543A CN114236541A CN 114236541 A CN114236541 A CN 114236541A CN 202111489543 A CN202111489543 A CN 202111489543A CN 114236541 A CN114236541 A CN 114236541A
Authority
CN
China
Prior art keywords
deformation
area
dimensional deformation
imaging
image
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
CN202111489543.0A
Other languages
English (en)
Other versions
CN114236541B (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.)
University of Electronic Science and Technology of China
Original Assignee
University of Electronic Science and Technology of China
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 University of Electronic Science and Technology of China filed Critical University of Electronic Science and Technology of China
Priority to CN202111489543.0A priority Critical patent/CN114236541B/zh
Publication of CN114236541A publication Critical patent/CN114236541A/zh
Application granted granted Critical
Publication of CN114236541B publication Critical patent/CN114236541B/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
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • Radar, Positioning & Navigation (AREA)
  • General Physics & Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Mathematical Analysis (AREA)
  • Computational Mathematics (AREA)
  • Algebra (AREA)
  • Mathematical Optimization (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Pure & Applied Mathematics (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Electromagnetism (AREA)
  • Image Processing (AREA)
  • Radar Systems Or Details Thereof (AREA)

Abstract

本发明公开了一种基于Sentinel‑1卫星SAR图像的大面积地表三维形变计算方法，该方法通过公开免费的对地遥感观测Sentinel卫星数据网站获取SAR图像数据文件进行处理，判断监测地区的SAR图像轨道类型，然后采取对应的地表三维形变计算方法，解算得到监测区域的地表三维形变序列。本发明支持大面积地表形变监测，解决高时效性、低成本地表三维形变序列计算问题。相对于传统方法，本发明方法可根据监测地区的不同轨道类型，给出不同的地表三维形变计算方法，加之利用该地区的所有可用SAR图像，提高地表三维形变计算的时效性。

Description

基于Sentinel-1卫星SAR图像的大面积地表三维形变计算 方法
技术领域
本发明涉及星载SAR成像处理的地表形变监测技术领域，尤其涉及一种基于Sentinel-1卫星SAR图像的大面积地表三维形变计算方法。
背景技术
地质灾害对我国影响巨大，据自然资源部统计，2020年全国共发生地质灾害7840起，共造成139人死亡(失踪)、58人受伤，直接经济损失50.2亿元。地质灾害类型主要以滑坡、崩塌和泥石流为主：其中滑坡4810起、崩塌1797起和泥石流899起。这些地质灾害孕育与发生均在该区域地表形变有直接体现。因此，监测地表形变数据进而预测地质灾害风险发生一直以来都是地质灾害监测预警研究的热点领域之一。星载SAR(合成孔径雷达)具有全天候、大范围、持续对地观测特点，被广泛应用于地表形变监测领域。在地质灾害监测中，由于地表三维形变可以真实全面地反映地质灾害的地表变化情况，地表三维形变监测对于地质灾害(如滑坡、泥石流、崩塌等)监测预警有着非常重要的作用。
在地表三维形变监测方面，由于传统星载SAR监测存在视线向(line of sight,LOS)模糊的问题。国内外研究者大都使用多轨道SAR图像解算地表三维形变，如通过使用SBAS-InSAR技术获取不同轨道SAR图像的LOS向形变，再将其进行投影计算得到地表三维形变。部分专家也提出一种从两平行轨道SAR图像解算地表三维移动量的方法，该方法突破了使用三个交叉轨道SAR图像才能计算地表三维形变的限制，即仅使用两平行轨道SAR图像即可解算得到地表三维形变。
此外，还有研究者提出使用一个轨道SAR图像计算三维形变序列的方法。该方法通过SBAS-InSAR技术获取视线向形变序列，再通过偏移量追踪技术获得方位向和距离向形变序列，进而解算出监测区域的地表三维形变序列。但是由于偏移量追踪技术受限于SAR图像的分辨率，该方法计算结果精度较低，难以满足地表形变监测的应用要求。因此该方法一般只用于地震等场景的监测。
上述这些方法在解算地表三维形变过程中均局限于某一特定的轨道类型。当监测区域范围较大时，该区域存在多种轨道类型的SAR图像，传统方法难以满足该场景的地表形变监测需求。此外，监测区域的SAR图像来源于多个轨道，每个轨道的SAR图像成像时间不同，上述方法只能获取其公共时间上的三维形变，难以获取整个时间段上的形变。此外，每个卫星的重访周期较长(如单个Sentinel-1卫星的重访时间为12天)，这些都影响地表三维形变计算的时效性。
本发明将解决现有地表三维形变计算方法适用性不强、计算时效性不足等问题，目标是实现一种低成本的大面积范围地表三维形变计算方法，可应用于地质灾害(如滑坡、泥石流等)监测预警。
针对地质灾害监测过程难以获取大面积且时效性高的地表三维形变序列问题，本发明采用偏移量追踪技术(offset-tracking，OT)、差分干涉测量小基线集时序分析(SmallBaseline Subset InSAR,SBAS-InSAR)技术对Sentinel-1卫星SAR图像进行处理以获取监测区域的三维地表形变数据，再基于滑动窗口处理方式，计算得到地表的三维形变序列，从而为地质灾害监测预警提供重要指标数据。本发明提出的地表三维形变计算方法可以适用于大面积地表形变监测处理，并解决多轨道SAR图像因成像时间不同导致地表形变计算时效性不足等问题。
发明内容
现有SAR图像解算地表三维形变技术只能针对单一的卫星轨道类型进行处理。当需要监测大面积地表形变时，观测区域存在多种卫星轨道类型SAR图像，传统的SAR图像三维形变计算方法存在一定局限性。
除此之外，现有SAR图像的地表三维形变计算方法大多都需要至少两个轨道的SAR图像。因不同轨道的SAR卫星观测成像时间不同，现有方法只能解算公共时间段上的三维形变。即便已经获取到了一个轨道最新的SAR图像，但若没有相同时段的其它轨道SAR图像数据，仍无法解算得到完整的三维形变序列，这会降低三维形变序列解算的时效性。例如，覆盖监测地区X的卫星在A轨道上有两幅不同的SAR图像A1、A2，在B轨道上有两幅不同的SAR图像B1和B2，且四幅图像的成像时间分别为
Figure BDA0003398642600000021
和
Figure BDA0003398642600000022
其中
Figure BDA0003398642600000023
由于现有技术计算形变需要两幅SAR图像必须来源于同一成像区域，因此只能计算出A轨道上
Figure BDA0003398642600000024
到
Figure BDA0003398642600000025
的LOS向形变、距离向形变或方位向形变和B轨道上
Figure BDA0003398642600000026
到
Figure BDA0003398642600000027
的LOS向形变、距离向形变或方位向形变。在监测区域X上进行多轨道SAR图像解算地表三维形变过程中，由于除
Figure BDA0003398642600000028
到
Figure BDA0003398642600000029
外的其他时间段没有同时具备两轨道的图像，也就无法使用这种方法解算三维形变。而类似
Figure BDA00033986426000000210
到
Figure BDA00033986426000000211
这一段同时具备所有轨道计算出形变值的时间段为公共时间段。
另一方面，当进行大范围地表形变监测时，从商业卫星获取SAR图像数据会有较高成本，为了降低成本，本发明方法采用Sentinle-1卫星网站提供的免费SAR图像数据进行地表形变监测。
针对上述问题，本发明提出一种基于Sentinel-1 A/B卫星多轨道类型SAR图像的地表三维形变计算方法。该方法支持大面积地表形变监测，解决高时效性、低成本地表三维形变序列计算问题。相对于传统方法，本发明方法可根据监测地区的不同轨道类型，给出不同的地表三维形变计算方法，加之利用该地区的所有可用SAR图像，提高地表三维形变计算的时效性。
本发明提出了一种基于Sentinel-1卫星SAR图像的大面积地表三维形变计算方法。该方法通过公开免费的对地遥感观测Sentinel卫星数据网站(https://search.asf.alaska.edu/#/)获取SAR图像数据文件进行处理，判断监测地区的SAR图像轨道类型，然后采取对应的地表三维形变计算方法，解算得到监测区域的地表三维形变序列。
1.选定监测地点
从公开免费的对地遥感观测Sentinel卫星数据网站(https://search.asf.alaska.edu/#/)选定待监测地表三维形变的监测地区，配置数据参数，准备下载该待监测地区的SAR图像文件。
2.定时下载SAR图像
SAR图像包括SLC模式和GRD模式的图像文件，每天定时从公开免费的Sentinel卫星数据网站(https://search.asf.alaska.edu/#/)获取待监测地区的单视复数图像文件(Single Look Complex，SLC)和地距多视图像文件(Ground Range Detected，GRD)。在后期数据处理中，使用SBAS-InSAR技术对SLC图像文件处理，获取视线向形变；使用OT技术对GRD图像文件处理，获取方位向形变或距离向形变。
3.判断轨道类型
在下载了待监测地区的SAR图像后，可以判断其轨道类型。在地球的不同地区，Sentin el-1卫星经过的频次和成像角度有所不同，即轨道数量与轨道之间关系(包括平行和交叉)不同。现有方法仅针对特定轨道类型进行地表三维形变计算，而没有考虑不同地区有不同轨道类型的处理情况。针对不同地区的轨道覆盖情况，可以将不同地区的轨道覆盖情况归类四种情况，分别为一个轨道覆盖类型的地区、两个交叉轨道覆盖的地区、两个平行轨道覆盖的地区和三个或更多轨道覆盖地区。
4.对应解算方案
本发明针对不同轨道类型的地区提出了不同的地表三维形变解算方案。
a.针对一个轨道覆盖类型的地区，其地表三维形变解算方案为：
S表示卫星，dx，dy，dz分别代表地表的南北方向、东西方向以及垂直地面方向上的形变。dLOS，ddistance，dazimu分别代表视线向、距离向(水平面上垂直于卫星飞行方向)以及方位向(卫星飞行方向)上的形变，且均由卫星S多次经过当前地区对地观测形成的多幅SAR图像计算得来。本发明使用SBAS-InSAR技术计算卫星S视线向的形变量dLOS，再使用OT技术计算出方位向形变dazimuth和距离向形变ddistance。然后再将视线向形变dLOS、方位向形变dazimuth和距离向形变ddistance投影到包含东西方向、南北方向和垂直方向的地表三维坐标系，即可解算出地表的南北方向、东西方向以及垂直地面方向上的形变。
b.针对两个交叉轨道覆盖的地区，其地表三维形变解算方案为：
S1，S2表示两颗轨道交叉的卫星，通过SBAS-InSAR技术分别计算出卫星S1和卫星S2的视线向形变dLOS1和dLOS2，再使用OT技术计算出距离向形变ddistance。将dLOS1、dLOS2和ddistance分别投影到包含东西方向、南北方向和垂直方向的地表三维坐标系，即可解算出地表的南北方向、东西方向以及垂直地面方向上的形变。本发明仅使用距离向形变而非方位向形变的原因是：OT技术受分辨率的影响较大，而Sentinel-1系列卫星用于地表形变监测的IW模式距离向分辨率更高，因此选用距离向形变会有更高的精度。
c.针对两个平行轨道覆盖的地区，其地表三维形变解算方案为：
S1，S2表示两颗轨道平行的卫星，通过SBAS-InSAR技术分别计算出卫星S1和卫星S2的视线向形变dLOS1和dLOS2，再使用OT技术计算出方位向形变dazimut。将dLOS1、dLOS2和dazimuth分别投影到包含东西方向、南北方向和垂直方向的地表三维坐标系，即可解算出地表的南北方向、东西方向以及垂直地面方向上的形变。在本方案中，仅使用方位向形变而非距离向形变的原因是：Sentinel-1卫星提供的SAR图像上的每一个像素点都可以认为是在垂直于卫星飞行方向的条件下获得的，视线向形变的方向与距离向形变的方向平行，因此两个平行轨道卫星的视线向形变和距离向形变处在同一个平面上，无法投影到东西、南北和垂直方向的地表三维坐标系上。
d.对于三个或更多轨道覆盖地区，其地表三维形变解算方案为：
S1，S2表示两颗轨道平行的卫星，S3表示与S1，S2轨道交叉的卫星。通过SBAS-InSAR技术计算出卫星S1、卫星S2和卫星S3各自的视线向形变dLOS、dLOS2和dLOS3，将dLOS1、dLOS2和dLOS3分别投影到东西方向、南北方向和垂直方向即可解算出地表的南北方向、东西方向以及垂直地面方向上的形变。
5.解算形变序列
在不同轨道类型的地表三维形变解算方案基础上，本发明提出了一种基于滑动窗口机制的处理方法，其目标是在提高时效性的同时下尽力提高精度。
对于一个地表形变监测地区，设其包含n个成像区域，ar为成像区域的编号，且ar＝1,2,…,n。
Figure BDA0003398642600000051
表示成像区域ar的第j幅SAR图像，
Figure BDA0003398642600000052
表示成像区域ar的第j幅SAR图像的获取时间。采集的地表形变监测地区SAR图像序列P如下：
Figure BDA0003398642600000053
为了计算整个图像序列P上的形变，同时能够获取尽量高精度的形变值，本发明方法采用滑动窗口机制。窗口大小k为1到n之间的整数。整体思想是：令k取遍1到n的所有值，当k＝n0时，将窗口从图像序列P起始点
Figure BDA0003398642600000054
开始向后移动，每次移动一幅SAR图像，直到
Figure BDA0003398642600000055
被纳入进窗口后结束，共移动n-n0次。在移动过程中，每次移动均根据当前窗口内部包含的成像区域的轨道类型选用对应的地表三维形变解算方案，将计算出的形变按照时间顺序组织起来，覆盖上一次遍历的结果，逐步得到最终的地表三维形变序列。
达到的效果是：当滑动窗口大小为1时，对所有的成像区域进行遍历，计算得到了所有图像获取时间上的形变序列，相较于传统方法只能获取公共时间上的形变序列，显著提高了时效性。随着滑动窗口变大，每次参与计算的成像区域增多，可获得的LOS向(视线向)形变数量增多。由于LOS向形变的精度要高于距离向形变和方位向形变的精度，计算出的地表三维形变的精度也不断提高。
先设窗口k＝1，代表每次只在该待监测地区的一个成像区域上计算。从
Figure BDA0003398642600000056
开始，认为该待监测地区为一轨道覆盖地区。对成像区域1的所有图像
Figure BDA0003398642600000057
中的任意两幅SAR图像(即图像
Figure BDA0003398642600000058
和图像
Figure BDA0003398642600000059
a取遍1到n)使用一个轨道覆盖地区的地表三维形变解算方案，得到成像区域1的所有图像的获取时间
Figure BDA00033986426000000510
上的形变序列L1，向后移动窗口，对成像区域2的任意两幅SAR图像使用一个轨道覆盖地区的地表三维形变解算方案，得到
Figure BDA00033986426000000511
上的形变序列L2，以此类推，直到该待监测地区的成像区域n计算完成。
再设窗口k＝2，代表每次在该待监测地区的两个成像区域上计算地表三维形变。从
Figure BDA00033986426000000512
开始，窗口内部有
Figure BDA00033986426000000513
和
Figure BDA00033986426000000514
两幅SAR图像。如果窗口内部两幅SAR图像所属的两个成像区域属于同一轨道，则不进行操作。如果两个成像区域所属轨道相互平行，则使用两个平行轨道覆盖地区的地表三维形变解算方案计算成像区域1和2上的地表三维形变，得到
Figure BDA00033986426000000515
Figure BDA0003398642600000061
即公共时间段上的形变序列L1-2。如果两个成像区域所属轨道交叉，则使用两个交叉轨道覆盖地区的地表三维形变解算方案计算成像区域1和2上的地表三维形变，得到
Figure BDA0003398642600000062
即公共时间段上的形变序列L1-2。再移动窗口，以此类推，直到该待监测地区的成像区域n计算完成，采用此时得到的公共时间段上的形变更新k＝1时得到的形变序列在该公共时间段上的形变。
再令k＝3，代表每次在该待监测地区的三个成像区域上计算地表三维形变。从
Figure BDA0003398642600000063
开始，窗口内部有
Figure BDA0003398642600000064
和
Figure BDA0003398642600000065
三幅SAR图像。如果窗口内部三幅SAR图像所属的三个成像区域对应三个轨道两两平行，则不进行操作。如果三个成像区域所属轨道存在交叉，则使用多个轨道覆盖地区的地表三维形变解算方案计算成像区域1、2、3上的地表三维形变，得到
Figure BDA0003398642600000066
Figure BDA0003398642600000067
即公共时间段上的形变序列L1-3。再移动窗口，以此类推，直到该待监测地区的成像区域n计算完成，进一步采用此时得到的公共时间段上的形变更新k＝2时得到的形变序列在该公共时间段上的形变。
再依次令k＝4，5……n，按照上述思路逐步计算，得到最终的地表三维形变序列L。
本发明在使用Sentinel-1A/B卫星SAR图像数据解算地表三维形变的过程中，充分考虑到监测区域的多种轨道类型，对于不同轨道类型采用了相应的地表三维形变序列解算方案，使得所有Sentinel-1卫星覆盖地区都可以进行地表三维形变计算，大大提高了SBAS-InSAR技术和OT技术在地表三维形变监测领域的应用范围和适应性。
本发明支持不同轨道覆盖类型地区的地表三维形变解算。以三轨覆盖地区为例，在公共时间段内使用三个轨道上卫星的视线向形变计算。在非公共时间段内，使用两个轨道上卫星的视线向形变结合方位向或距离向形变进行地表三维形变解算。在非公共时间段内，还可仅使用一个轨道上卫星的视线向形变结合方位向形变和距离向形变进行地表三维形变解算。本发明方法实现了不同轨道覆盖地区的所有时间段均可解算出地表三维形变序列，提高了地表三维形变序列解算的时效性。同时本发明使用开源SAR图像数据实现了低成本的大面积地表三维形变计算处理。
附图说明
图1为本发明地表三维形变计算方法总体流程图；
图2为本发明不同地区的轨道覆盖情况的示意图；
图3为本发明一个轨道覆盖地区的地表三维形变解算方案示意图；
图4为本发明两轨道交叉覆盖地区的地表三维形变解算方案示意图；
图5为本发明两平行轨道覆盖地区的地表三维形变解算方案示意图；
图6为本发明多轨道覆盖地区的地表三维形变解算方案示意图；
图7为本发明SAR图像下载相关参数选择界面；
图8为本发明SAR图像下载脚本获取界面；
图9为本发明SAR图像查看界面；
图10为本发明A,B,C,D四个成像区域示意图。
具体实施方式
下面将结合附图对本发明实施例中的技术方案进行清楚、完整地描述。
如图1所示，本发明提出了一种基于Sentinel-1卫星SAR图像的大面积地表三维形变计算方法，该方法具体包括如下步骤：
1.选定监测地点
在提供公开SAR图像数据的ASF Data Search网站(https://search.asf.alaska.edu/#/)中，选取需要对地观测地表形变的监测地区。
2.定时下载SAR图像
在ASF Data Search网站中，上传待监测地区位置信息(如shp文件、geojson文件、kml文件之一)，限定下载图像的范围，点击“Filters”按钮，如图7所示：
选择需要下载SAR图像的起止日期，“File Type”选择SLC和GRD类型分别用于SBAS-InSAR和OT技术。“Beam Mode”选择IW模式，该模式主要用于陆地观测。“Direction”升降轨全选，图7中“Subtype”表示需要获取的卫星图像来源，当前代表同时获取Sentinel-1A/B卫星图像。在点击搜索按钮后，就可以检索出符合条件的SAR图像。点击“Queue”将SAR图像加入下载队列，即可在“Download”界面中获取到下载图像的python脚本，如图8所示，执行该脚本即可获取相关SAR图像。
3.判断轨道类型
在下载SAR图像前，可以先查看需要对地观测地表形变的待监测地区所有的SAR图像文件，显示结果如下图9所示：
在地图中，椭圆框所示区域包含多个成像区域。点击成像区域后会展示出图像的Frame编号和Path编号。当两个成像区域Path编号相同而Frame编号不同则代表这两幅图像属于同一轨道的不同成像区域。两个成像区域Path编号不同且两个成像区域的矩形平行时，则代表这两个轨道平行。当Path编号不同且成像区域矩形不平行时，则代表这两个轨道交叉。如果有三个成像区域Path编号不同，且对应的矩形框不两两平行，则说明这三个轨道交叉。
针对不同地区的轨道覆盖情况，可以将不同地区的轨道覆盖情况归类为图2所示的四种情况，这四种情况分别为一个轨道覆盖类型的地区、两个交叉轨道覆盖的地区、两个平行轨道覆盖的地区和三个或更多轨道覆盖地区。
图2中每一个方框都代表对待监测地区的一个成像。同一监测地区因卫星轨道和入射角度不同，会产生多个成像。在Sentinel卫星数据网站上可以查看每一幅图像的轨道信息等数据。Sentinel-1卫星在不同地区的成像情况有所不同，但它们都属于图2的四种类型之一，图2中a)表示一个轨道覆盖地区，b)表示两个交叉轨道覆盖地区，c)表示两个平行轨道覆盖地区，d)表示三轨(多个轨道)覆盖地区。在后续步骤中，给出针对每种类型的地表三维形变序列给出解算方案。
4.对应解算方案
在图3到图6所示的不同轨道图像解算方案中，将使用SBAS-InSAR技术和OT技术。针对一个轨道覆盖类型的地区，其地表三维形变解算方案如图3所示，针对两个交叉轨道覆盖的地区，其地表三维形变解算方案如图4所示，针对两个平行轨道覆盖的地区，其地表三维形变解算方案如图5所示，对于三个或更多轨道覆盖地区，其地表三维形变解算方案如图6所示。
SBAS-InSAR是在差分InSAR基础上发展起来的一种新时间序列InSAR分析方法，能够降低相位噪声和误差。
假定在时间t1至tM内获取同一地区的M幅SAR图像，然后根据干涉组合条件，在短基线距的条件下形成N幅干涉条纹图，且满足
Figure BDA0003398642600000081
对tA和tB(t1<tA<tB<tM)时刻两幅SAR图像生成的第nin(nin＝1,2,…,N)幅干涉条纹图，在去除平地及地形相位影响后，第nin幅干涉条纹图中的第x个像素的干涉相位可表示为：
Figure BDA0003398642600000082
式中，
Figure BDA0003398642600000083
为tA～tB对应的视线向形变，
Figure BDA0003398642600000084
为tA～tB对应的地形相位误差，
Figure BDA0003398642600000085
为tA～tB对应的大气相位误差，
Figure BDA0003398642600000086
为tA～tB对应的基线轨道引起的相位误差；
Figure BDA0003398642600000087
为tA～tB对应的噪声误差。
对N辐干涉条纹图进行相位解缠、去除平地效应、去除大气相位、轨道矫正后即可获取t1至tM的视线向形变dLOS。
OT技术则利用强度相关对两幅SAR图像逐像素配准，并从配准偏移量中解算地表位移。对于Sentinel-1卫星而言，在进行轨道矫正之后可以近似认为配准偏移量即为地表偏移量，再将其分解到方位向和距离向，就得到了方位向形变dazimuth和距离向形变ddistance。
令卫星的入射角为θ，轨道倾角为α，θi、αi分别第i个轨道对应卫星的入射角和轨道倾角。则有：
Figure BDA0003398642600000091
dU为垂直水平面向上的地表形变，dE为东西方向上的地表形变，dN为南北方向上的地表形变。
Figure BDA0003398642600000092
为使用SBAS-InSAR计算出的第i个轨道的视线向形变。对于OT技术计算出的方位向形变dazimuth和距离向形变ddistance有：
dazimut＝dEcos(αi-3π/2)+dNsin(αi-3π/2)
ddistance＝dEsin(αi-3π/2)+dNcos(αi-3π/2)
对于一个轨道覆盖的类型，解下方方程组即可得到dU、dE、dN：
dLOS＝dUcosθ-dEsinθsin(α-3π/2)-dNsinθcos(α-3π/2)
dazimuth＝dEcos(α-3π/2)+dNsin(α-3π/2)
ddistance＝dEsin(α-3π/2)+dNcos(α-3π/2)
对于两个轨道交叉覆盖的类型，解下方方程组即可得到dU、dE、dN：
Figure BDA0003398642600000093
Figure BDA0003398642600000094
ddistance＝dEsin(α-3π/2)+dNcos(α-3π/2)
对于两个轨道平行覆盖的类型，解下方方程组即可得到dU、dE、dN：
Figure BDA0003398642600000095
Figure BDA0003398642600000096
dazimuth＝dEcos(α-3π/2)+dNsin(α-3π/2)
对于多轨道覆盖的类型，以三轨道为例，解下方方程组即可得到dU、dE、dN：
Figure BDA0003398642600000097
Figure BDA0003398642600000098
Figure BDA0003398642600000101
5.解算形变序列
为便于说明，此处选择以需要对地观测地表形变的待监测地区中某四个成像区域A、B、C、D的情况进行说明，其中A、B成像区域轨道共线，A、B成像区域轨道与C和D成像区域所处轨道交叉，C、D成像区域轨道平行，如图10所示。假设A、B、C、D四个成像区域的成像时间均匀分布。
图中方框为对应成像区域，成像区域右侧边线为其所属轨道。观测的图像序列如下：
A1，B1，C1，D1，A2，B2，C2，D2
A1表示A的第一幅图像，
Figure BDA0003398642600000102
表示其获取时间，(A2表示A的第二幅图像，
Figure BDA0003398642600000103
表示其获取时间，B1表示B的第一幅图像，
Figure BDA0003398642600000104
表示其获取时间，B2表示B的第二幅图像，
Figure BDA0003398642600000105
表示其获取时间，C1表示C的第一幅图像，
Figure BDA0003398642600000106
表示其获取时间，C2表示C的第二幅图像，
Figure BDA0003398642600000107
表示其获取时间，D1表示D的第一幅图像，
Figure BDA0003398642600000108
表示其获取时间，D2表示D的第二幅图像，
Figure BDA0003398642600000109
表示其获取时间)最新图像为D2。设k代表滑动窗口大小(1≤k≤4，k为整数)。
令k＝1，从A1开始，对成像区域A使用一轨道覆盖类型的地表三维形变解算方案，即使用A1和A2获取
Figure BDA00033986426000001010
到
Figure BDA00033986426000001011
的地表三维形变d1(这里将该成像区域A的地表三维形变标记为d1，实质上包含东西向、南北向和垂直地面向上的地表三维形变)，再假定每次获取图像的时间间隔相同，就可以得到
Figure BDA00033986426000001012
到
Figure BDA00033986426000001013
之间的形变为d1/4，
Figure BDA00033986426000001014
到
Figure BDA00033986426000001015
之间的形变为d1/2，
Figure BDA00033986426000001016
到
Figure BDA00033986426000001017
之间的形变为3d1/4。移动窗口至B1，获取
Figure BDA00033986426000001018
到
Figure BDA00033986426000001019
的地表三维形变d2，
Figure BDA00033986426000001020
到
Figure BDA00033986426000001021
之间的形变为d1/4+d2。移动窗口至C1，获取
Figure BDA00033986426000001022
到
Figure BDA00033986426000001023
的地表三维形变d3,
Figure BDA00033986426000001024
到
Figure BDA00033986426000001025
之间的形变为d1/2+d3。移动窗口至D1，获取
Figure BDA00033986426000001026
到
Figure BDA00033986426000001027
的地表三维形变d4,
Figure BDA00033986426000001028
到
Figure BDA00033986426000001029
之间的形变为3d1/4+d4。则形变序列为(相对于时间
Figure BDA00033986426000001030
的形变)：
Figure BDA00033986426000001031
令k＝2，采用k＝2得到的公共时间段上的形变更新k＝1得到的形变序列在该公共时间段上的形变。从A1，B1开始，由于A、B共线，仍相当于一个轨道覆盖类型的地区，因此不进行操作。移动窗口至B1，C1，B、C所属轨道交叉，相当于两个交叉轨道覆盖类型的地区，对B1，C1和B2，C2获取
Figure BDA00033986426000001032
到
Figure BDA00033986426000001033
(
Figure BDA00033986426000001034
到
Figure BDA00033986426000001035
和
Figure BDA00033986426000001036
到
Figure BDA00033986426000001037
的公共时间段)的地表三维形变d5,
Figure BDA00033986426000001038
到
Figure BDA00033986426000001039
之间的形变为d1/2+d5。C、D轨道平行，相当于两个平行轨道覆盖类型的地区对C1，D1和C2，D2获取
Figure BDA0003398642600000111
到
Figure BDA0003398642600000112
(
Figure BDA0003398642600000113
到
Figure BDA0003398642600000114
和
Figure BDA0003398642600000115
到
Figure BDA0003398642600000116
的公共时间段)的地表三维形变d6，
Figure BDA0003398642600000117
到
Figure BDA0003398642600000118
之间的形变为3d1/4+d6，则形变序列变为：
Figure BDA0003398642600000119
令k＝3，采用k＝3得到的公共时间段上的形变进一步更新k＝2得到的形变序列在该公共时间段上的形变。从A1，B1，C1开始，由于A、B共线，相当于两个交叉轨道覆盖类型的地区，不进行操作，移动窗口到B1，C1，D1，相当于多个轨道覆盖类型的地区，对B1，C1，D1和B2，C2，D2获取
Figure BDA00033986426000001110
到
Figure BDA00033986426000001111
的地表三维形变d7，
Figure BDA00033986426000001112
到
Figure BDA00033986426000001113
之间的形变为3d1/4+d7，则形变序列变为：
Figure BDA00033986426000001114
令k＝4，形变序列的更新规则和之前一样。由于此时仍然只包含三个不同方向的轨道，因此和k＝3时等同。最终形变序列为：
Figure BDA00033986426000001115
相较于传统方法，由于其只能获取整个图像序列的公共时间段，即
Figure BDA00033986426000001116
到
Figure BDA00033986426000001117
上的地表三维形变，本发明的方法获取了所有时间上的形变序列，提高了形变数据的时效性。
本发明使用开源免费的Sentinel-1 A/B双星座卫星SAR图像数据进行地表三维形变序列的解算，这在进行大面积地表形变监测中，可极大降低使用商业卫星数据成本。分析待监测地区的卫星轨道覆盖类型，对于不同轨道类型的地区使用不同的三维形变解算方案。将卫星轨道覆盖地区分为：一轨覆盖地区、两轨平行覆盖地区、两轨交叉覆盖地区以及多轨覆盖(三轨及以上)地区。对于一轨覆盖地区，使用SBAS-InSAR技术获取LOS向形变，使用OT技术获取距离向和方位向形变。对于两轨平行覆盖地区，使用SBAS-InSAR技术获取两个LOS向形变，使用OT技术获取方位向形变。对于两轨交叉覆盖地区，使用SBAS-InSAR技术获取两个LOS向形变，使用OT技术获取距离向形变。对于多轨道覆盖地区(三轨及以上)，使用SBAS-InSAR技术获取三个LOS向形变，在已经获取三个方向不共面的形变之后，对地表三维坐标进行投影解算，得到该地区的地表三维形变。将监测区域的多个SAR图像结合在一起进行地表三维形变序列解算。在地表形变序列解算过程中，针对不同轨道类型图像使用不同的解算方案，逐幅图像计算三维形变。同时使用滑动窗口机制，将其合并为形变序列，实现了公共时间段和非公共时间段上的三维形变序列计算。
采用上述技术方案，本发明针对大面积监测区域中存在多种轨道类型SAR图像，使用不同的地表三维形变解算方案，不再局限于某一种特殊轨道类型，这可提高三维形变计算方法在大面积地区地表形变监测的适应性。使用监测地区所有成像区域的SAR图像进行形变计算，保证了较高频率进行图像处理，为高时效性地解算三维形变序列提供支持；基于滑动窗口机制进行形变序列解算处理，保证了形变序列解算的高时效性，从而确保形变解算的精度：当滑动窗口大小为1时，计算得到所有已获取图像之间的形变序列；逐步增大滑动窗口尺寸，滑动窗口内部的成像区域数量增加，滑动窗口可支持更多的轨道类型图像处理，并可使用高精度的形变覆盖低精度的形变：使用多轨道覆盖地区的地表三维形变解算方案计算出的形变覆盖其他方案的计算结果；使用两交叉轨道覆盖地区的地表三维形变解算方案计算出的形变覆盖两平行轨道覆盖地区的地表三维形变解算方案、一个轨道覆盖地区的地表三维形变解算方案计算出的形变；使用两平行轨道覆盖地区的地表三维形变解算方案计算出的形变覆盖一个轨道覆盖地区的地表三维形变解算方案计算出的形变。
以上所述，仅为本发明的具体实施方式，本说明书中所公开的任一特征，除非特别叙述，均可被其他等效或具有类似目的的替代特征加以替换；所公开的所有特征、或所有方法或过程中的步骤，除了互相排斥的特征和/或步骤以外，均可以任何方式组合；本领域的技术人员根据本发明技术方案的技术特征所做出的任何非本质的添加、替换，均属于本发明的保护范围。

Claims (5)

1.一种基于Sentinel-1卫星SAR图像的大面积地表三维形变计算方法，其特征在于，该方法具体包括如下步骤：
步骤1：选定监测地点
从公开免费的对地遥感观测Sentinel卫星数据网站选定待监测地表三维形变的监测地区，配置数据参数，准备下载该待监测地区的SAR图像文件；
步骤2：定时下载SAR图像
SAR图像包括SLC模式和GRD模式的图像文件，每天定时从所述对地遥感观测Sentinel卫星数据网站获取该待监测地区的单视复数图像文件SLC和地距多视图像文件GRD，在后期数据处理中，使用差分干涉测量小基线集时序分析SBAS-InSAR技术对SLC图像文件处理，获取视线向形变；使用偏移量追踪OT技术对GRD图像文件处理，获取方位向形变或距离向形变；
步骤3：判断轨道类型
在下载了待监测地区的SAR图像后，判断轨道类型，针对不同地区的轨道覆盖情况，分为四种情况：一个轨道覆盖类型的地区、两个交叉轨道覆盖的地区、两个平行轨道覆盖的地区和三个或更多轨道覆盖地区；
步骤4：对应解算方案
针对不同轨道类型的地区给出了不同的地表三维形变解算方案：
a.针对一个轨道覆盖类型的地区，其地表三维形变解算方案为：
S表示卫星，dLOS，ddistance，dazimuth分别代表视线向、距离向以及方位向上的形变，使用SBAS-InSAR技术计算卫星S视线向的形变量dLOS，再使用OT技术计算出方位向形变dazimuth和距离向形变ddistance，然后再将视线向形变dLOS、方位向形变dazimuth和距离向形变ddistance分别投影到包含东西方向、南北方向和垂直方向的地表三维坐标系，便能够解算出地表的南北方向、东西方向以及垂直地面方向上的形变；
b.针对两个交叉轨道覆盖的地区，其地表三维形变解算方案为：
S1，S2表示两颗轨道交叉的卫星，通过SBAS-InSAR技术分别计算出卫星S1和卫星S2的视线向形变dLOS1和dLOS2，再使用OT技术计算出距离向形变ddistance，将dLOS1、dLOS2和ddistance分别投影到包含东西方向、南北方向和垂直方向的地表三维坐标系，便能够解算出地表的南北方向、东西方向以及垂直地面方向上的形变；
c.针对两个平行轨道覆盖的地区，其地表三维形变解算方案为：
S1，S2表示两颗轨道平行的卫星，通过SBAS-InSAR技术分别计算出卫星S1和卫星S2的视线向形变dLOS1和dLOS2，再使用OT技术计算出方位向形变dazimuth，将dLOS1、dLOS2和dazimuth分别投影到包含东西方向、南北方向和垂直方向的地表三维坐标系，便能够解算出地表的南北方向、东西方向以及垂直地面方向上的形变；
d.对于三个或更多轨道覆盖地区，其地表三维形变解算方案为：
S1，S2表示两颗轨道平行的卫星，S3表示与S1，S2轨道交叉的卫星，通过SBAS-InSAR技术计算出卫星S1、卫星S2和卫星S3各自的视线向形变dLOS1、dLOS2和dLOS3，将dLOS1、dLOS2和dLOS3分别投影到东西方向、南北方向和垂直方向便能够解算出地表的南北方向、东西方向以及垂直地面方向上的形变；
步骤5：解算形变序列
在不同轨道类型的地区的地表三维形变解算方案基础上，提出一种基于滑动窗口机制的处理方法，具体包括：
对于一个待监测地表三维形变的监测地区，设其包含n个成像区域，ar为成像区域的编号，且ar＝1，2，...，n，
Figure FDA0003398642590000021
表示成像区域ar的第j幅SAR图像，
Figure FDA0003398642590000022
表示成像区域ar的第j幅SAR图像的获取时间，采集的该待监测地表三维形变的监测地区的SAR图像序列P如下：
Figure FDA0003398642590000023
为了计算整个图像序列P上的形变，采用滑动窗口机制，设置窗口大小k为1到n之间的整数，整体思想是：令k取遍1到n的所有值，当k＝n0时，将窗口从图像序列P起始点
Figure FDA0003398642590000024
开始向后移动，每次移动一幅SAR图像，直到
Figure FDA0003398642590000025
被纳入进窗口后结束，共移动n-n0次，在移动过程中，每次移动均根据当前窗口内部包含的成像区域的轨道类型选用对应的地表三维形变解算方案，将计算出的形变按照时间顺序组织起来，覆盖上一次遍历的结果，逐步得到最终的地表三维形变序列；
先设窗口k＝1，代表每次只在该待监测地区的一个成像区域上计算，从
Figure FDA0003398642590000026
开始，认为该待监测地区为一个轨道覆盖类型的地区，对成像区域1的所有图像
Figure FDA0003398642590000027
中的任意两幅SAR图像使用一个轨道覆盖类型的地区的地表三维形变解算方案，得到成像区域1的所有图像的获取时间
Figure FDA0003398642590000028
上的形变序列L1，向后移动窗口，对成像区域2的任意两幅SAR图像使用一个轨道覆盖类型的地区的地表三维形变解算方案，得到
Figure FDA0003398642590000029
上的形变序列L2，以此类推，直到该待监测地区的成像区域n计算完成；
再设窗口k＝2，代表每次在该待监测地区的两个成像区域上计算地表三维形变，从
Figure FDA0003398642590000031
开始，窗口内部有
Figure FDA0003398642590000032
和
Figure FDA0003398642590000033
两幅SAR图像，如果窗口内部两幅SAR图像所属的两个成像区域属于同一轨道，则不进行操作；如果两个成像区域所属轨道相互平行，则使用两个平行轨道覆盖的地区的地表三维形变解算方案计算成像区域1和2上的地表三维形变，得到
Figure FDA0003398642590000034
Figure FDA0003398642590000035
即公共时间段上的形变序列L1-2；如果两个成像区域所属轨道交叉，则使用两个交叉轨道覆盖的地区的地表三维形变解算方案计算成像区域1和2上的地表三维形变，得到
Figure FDA0003398642590000036
即公共时间段上的形变序列L1-2；再移动窗口，以此类推，直到该待监测地区的成像区域n计算完成，采用本次得到的公共时间段
Figure FDA0003398642590000037
Figure FDA0003398642590000038
上的形变序列L1-2覆盖k＝1时得到的
Figure FDA0003398642590000039
上的形变序列；
再令k＝3，代表每次在该待监测地区的三个成像区域上计算地表三维形变，从
Figure FDA00033986425900000310
开始，窗口内部有
Figure FDA00033986425900000311
和
Figure FDA00033986425900000312
三幅SAR图像，如果窗口内部三幅SAR图像所属的三个成像区域对应三个轨道两两平行，则不进行操作；如果三个成像区域所属轨道存在交叉，则使用三个或更多轨道覆盖地区的地表三维形变解算方案计算成像区域1、2、3上的地表三维形变，得到
Figure FDA00033986425900000313
即公共时间段上的形变序列L1-3；再移动窗口，以此类推，直到该待监测地区的成像区域n计算完成，进一步采用本次得到的公共时间段
Figure FDA00033986425900000314
Figure FDA00033986425900000315
上的形变更新k＝2时得到的
Figure FDA00033986425900000316
上的形变序列；
再依次令k＝4，5......n，按照上述思路逐步计算，直到得到最终的地表三维形变序列L。
2.根据权利要求1所述的基于Sentinel-1卫星SAR图像的大面积地表三维形变计算方法，其特征在于，所述步骤2中定时下载SAR图像具体包括：在所述对地遥感观测Sentinel卫星数据网站中，设置搜索待监测地区SAR图像的位置信息，点击“Filters”按钮，选择需要下载SAR图像的起止日期，“File Type”选择SLC和GRD类型分别用于SBAS-InSAR和OT技术，“Beam Mode”选择IW模式，该模式主要用于陆地观测，“Direction”升降轨全选，“Subtype”表示需要获取的卫星图像来源，当前代表同时获取Sentinel-1 A/B卫星图像，在点击搜索按钮后，检索出符合条件的SAR图像，点击“Queue”将SAR图像加入下载队列，在“Download”界面中获取到下载图像的python脚本，执行该脚本获取相关SAR图像。
3.根据权利要求2所述的基于Sentinel-1卫星SAR图像的大面积地表三维形变计算方法，其特征在于，所述步骤3中判断轨道类型具体包括：在下载SAR图像前，先查看待监测地表三维形变的监测地区所有的SAR图像文件，在显示结果的地图中，点击成像区域后会展示出图像的Frame编号和Path编号，当两个成像区域Path编号相同而Frame编号不同则代表这两幅图像属于同一轨道的不同成像区域；当两个成像区域Path编号不同且两个成像区域的矩形框平行时，则代表这两个轨道平行；当Path编号不同且成像区域矩形框不平行时，则代表这两个轨道交叉；如果有三个成像区域Path编号不同，且对应的矩形框不两两平行，则说明这三个轨道交叉。
4.根据权利要求3所述的基于Sentinel-1卫星SAR图像的大面积地表三维形变计算方法，其特征在于，所述步骤4中不同的地表三维形变解算方案具体为：
采用SBAS-InSAR技术降低相位噪声和误差，假定在时间t1至tM内获取同一地区的M幅SAR图像，然后根据干涉组合条件，在短基线距的条件下形成N幅干涉条纹图，且满足
Figure FDA0003398642590000041
对tA和tB时刻两幅SAR图像生成的第nin幅干涉条纹图，在去除平地及地形相位影响后，第nin幅干涉条纹图中的第x个像素的干涉相位表示为：
Figure FDA0003398642590000042
式中，
Figure FDA0003398642590000043
为tA～tB对应的视线向形变，
Figure FDA0003398642590000044
为tA～tB对应的地形相位误差，
Figure FDA0003398642590000045
为tA～tB对应的大气相位误差，
Figure FDA0003398642590000046
为tA～tB对应的基线轨道引起的相位误差；
Figure FDA0003398642590000047
为tA～tB对应的噪声误差，且t1＜tA＜tB＜tM，nin＝1，2，…，N；
对N辐干涉条纹图进行相位解缠、去除平地效应、去除大气相位、轨道矫正后便能够获取t1至tM的视线向形变dLOS；
OT技术则利用强度相关对两幅SAR图像逐像素配准，并从配准偏移量中解算地表位移，对于Sentinel-1卫星而言，在进行轨道矫正之后认为配准偏移量即为地表偏移量，再将其分解到方位向和距离向，就得到了方位向形变dazimuth和距离向形变ddistance；
令卫星的入射角为θ，其轨道倾角为α；θi、αi分别第i个轨道对应卫星的入射角和轨道倾角，则有：
Figure FDA0003398642590000048
dU为垂直方向上的地表形变，dE为东西方向上的地表形变，dN为南北方向上的地表形变，
Figure FDA0003398642590000051
为使用SBAS-InSAR技术计算出的第i个轨道的视线向形变，对于OT技术计算出的方位向形变dazimu和距离向形变ddistance有：
dazimuth＝dEcos(αi-3π/2)+dNsin(αi-3π/2)
ddistance＝dEsin(αi-3π/2)+dNcos(αi-3π/2)
对于一个轨道覆盖类型的地区，解下列方程组便能够得到dU、dE、dN：
dLOS＝dUcosθ-dEsinθsin(α-3π/2)-dNsinθcos(α-3π/2)
dazimuth＝dEcos(α-3π/2)+dNsin(α-3π/2)
ddistance＝dEsin(α-3π/2)+dNcos(α-3π/2)
对于两个交叉轨道覆盖的地区，解下列方程组便能够得到dU、dE、dN：
Figure FDA0003398642590000052
Figure FDA0003398642590000053
ddistance＝dEsin(α-3π/2)+dNcos(α-3π/2)
对于两个平行轨道覆盖的地区，解下列方程组便能够得到dU、dE、dN：
Figure FDA0003398642590000054
Figure FDA0003398642590000055
dazimuth＝dEcos(α-3π/2)+dNsin(α-3π/2)
对于三个或更多轨道覆盖地区，以三轨道为例，解下列方程组便能够得到dU、dE、dN：
Figure FDA0003398642590000056
Figure FDA0003398642590000057
Figure FDA0003398642590000058
5.根据权利要求4所述的基于Sentinel-1卫星SAR图像的大面积地表三维形变计算方法，其特征在于，所述步骤5中选择待监测地区中某四个成像区域A、B、C、D的情况，解算形变序列的具体内容包括：所述四个成像区域A、B、C、D中A、B成像区域轨道共线，A、B成像区域轨道与C和D成像区域所处轨道交叉，C、D成像区域轨道平行，并假设A、B、C、D四个成像区域的成像时间均匀分布，此时观测的图像序列如下：
A1，B1，C1，D1，A2，B2，C2，D2
其中，A1表示A的第1幅图像，
Figure FDA0003398642590000059
表示其获取时间；A2表示A的第2幅图像，
Figure FDA00033986425900000510
表示其获取时间；B1表示B的第1幅图像，
Figure FDA0003398642590000061
表示其获取时间；B2表示B的第2幅图像，
Figure FDA0003398642590000062
表示其获取时间；C1表示C的第1幅图像，
Figure FDA0003398642590000063
表示其获取时间；C2表示C的第2幅图像，
Figure FDA0003398642590000064
表示其获取时间；D1表示D的第1幅图像，
Figure FDA0003398642590000065
表示其获取时间；D2表示D的第2幅图像，
Figure FDA0003398642590000066
表示其获取时间；最新图像为D2；设k代表滑动窗口大小，且1≤k≤4，k为整数；
令k＝1，从A1开始，对成像区域A使用一个轨道覆盖类型的地区的地表三维形变解算方案，即使用A1和A2获取
Figure FDA0003398642590000067
到
Figure FDA0003398642590000068
的地表三维形变d1，所述d1包含了东西方向、南北方向和垂直方向上的地表三维形变；再假定每次获取图像的时间间隔相同，得到
Figure FDA0003398642590000069
到
Figure FDA00033986425900000610
之间的形变为d1/4，
Figure FDA00033986425900000611
到
Figure FDA00033986425900000612
之间的形变为d1/2，
Figure FDA00033986425900000613
到
Figure FDA00033986425900000614
之间的形变为3d1/4；移动窗口至B1，获取
Figure FDA00033986425900000615
到
Figure FDA00033986425900000616
的地表三维形变d2，
Figure FDA00033986425900000617
到
Figure FDA00033986425900000618
之间的形变为d1/4+d2；移动窗口至C1，获取
Figure FDA00033986425900000619
到
Figure FDA00033986425900000620
的地表三维形变d3，
Figure FDA00033986425900000621
到
Figure FDA00033986425900000622
之间的形变为d1/2+d3；移动窗口至D1，获取
Figure FDA00033986425900000623
到
Figure FDA00033986425900000624
的地表三维形变d4，
Figure FDA00033986425900000625
到
Figure FDA00033986425900000626
之间的形变为3d1/4+d4，得到此时的地表三维形变序列L1；
然后，令k＝2，采用k＝2得到的公共时间段上的形变更新k＝1时得到的形变序列在该公共时间段上的形变；从A1，B1开始，由于A、B共线，仍相当于一个轨道覆盖类型的地区，因此不进行操作；移动窗口至B1，C1，B、C所属轨道交叉，相当于两个交叉轨道覆盖类型的地区，对B1，C1和B2，C2获取
Figure FDA00033986425900000627
到
Figure FDA00033986425900000628
公共时间段上的地表三维形变d5，则
Figure FDA00033986425900000629
到
Figure FDA00033986425900000630
之间的形变为d1/2+d5；由于C、D轨道平行，相当于两个平行轨道覆盖类型的地区，对C1，D1和C2，D2获取
Figure FDA00033986425900000631
到
Figure FDA00033986425900000632
公共时间段上的地表三维形变d6，则
Figure FDA00033986425900000633
到
Figure FDA00033986425900000634
之间的形变为3d1/4+d6，得到此时的地表三维形变序列L2；
再令k＝3，采用k＝3得到的公共时间段上的形变进一步更新k＝2时得到的形变序列在该公共时间段上的形变；从A1，B1，C1开始，由于A、B共线，相当于两个交叉轨道覆盖类型的地区，不进行操作；移动窗口到B1，C1，D1，相当于多个轨道覆盖类型的地区，对B1，C1，D1和B2，C2，D2获取
Figure FDA00033986425900000635
到
Figure FDA00033986425900000636
的公共时间段上的地表三维形变d7，则
Figure FDA00033986425900000637
到
Figure FDA00033986425900000638
之间的形变为3d1/4+d7，得到此时的地表三维形变序列L3；
最后令k＝4，形变序列的更新规则和之前一样，由于此时仍然只包含三个不同方向的轨道，因此和k＝3时等同，得到最终的地表三维形变序列L。
CN202111489543.0A 2021-12-08 2021-12-08 基于Sentinel-1卫星SAR图像的大面积地表三维形变计算方法 Active CN114236541B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111489543.0A CN114236541B (zh) 2021-12-08 2021-12-08 基于Sentinel-1卫星SAR图像的大面积地表三维形变计算方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111489543.0A CN114236541B (zh) 2021-12-08 2021-12-08 基于Sentinel-1卫星SAR图像的大面积地表三维形变计算方法

Publications (2)

Publication Number Publication Date
CN114236541A true CN114236541A (zh) 2022-03-25
CN114236541B CN114236541B (zh) 2023-05-16

Family

ID=80753810

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111489543.0A Active CN114236541B (zh) 2021-12-08 2021-12-08 基于Sentinel-1卫星SAR图像的大面积地表三维形变计算方法

Country Status (1)

Country Link
CN (1) CN114236541B (zh)

Citations (15)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106023157A (zh) * 2016-05-10 2016-10-12 电子科技大学 一种基于sar图像的山区地表微形变信息提取方法
CN106526590A (zh) * 2016-11-04 2017-03-22 山东科技大学 一种融合多源sar影像工矿区三维地表形变监测及解算方法
CN107144213A (zh) * 2017-06-29 2017-09-08 中南大学 基于sar强度影像的矿区大量级三维时序形变估计方法及装置
CN108627832A (zh) * 2018-05-11 2018-10-09 电子科技大学 一种基于多时序sar图像提取输电通道地表形变的方法
CN108983232A (zh) * 2018-06-07 2018-12-11 中南大学 一种基于邻轨数据的InSAR二维地表形变监测方法
CN110058237A (zh) * 2019-05-22 2019-07-26 中南大学 面向高分辨率SAR影像的InSAR点云融合及三维形变监测方法
CN110058234A (zh) * 2019-05-20 2019-07-26 太原理工大学 一种解算矿区地表沉降三维形变的方法
CN110161501A (zh) * 2019-05-24 2019-08-23 电子科技大学 一种多时序sar图像的目标区域地表起伏信息提取方法
CN110888130A (zh) * 2019-10-30 2020-03-17 华东师范大学 一种基于升降轨时序InSAR的煤矿区地表形变监测方法
CN111458709A (zh) * 2020-06-08 2020-07-28 河南大学 一种星载雷达广域地表二维形变场监测方法及装置
CN111474544A (zh) * 2020-03-04 2020-07-31 广东明源勘测设计有限公司 一种基于sar数据的滑坡形变监测及预警方法
CN111650579A (zh) * 2020-06-12 2020-09-11 中南大学 一种岩移参数自适应获取的InSAR矿区三维形变估计方法、装置及介质
CN112068136A (zh) * 2020-09-14 2020-12-11 广东省核工业地质局测绘院 一种基于幅度偏移量的方位向形变监测方法
CN112904337A (zh) * 2021-05-07 2021-06-04 北京东方至远科技股份有限公司 一种基于Offset Tracking技术的边坡形变时序监测方法
CN113091599A (zh) * 2021-04-06 2021-07-09 中国矿业大学 融合无人机dom和星载sar影像的地表三维形变提取方法

Patent Citations (15)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106023157A (zh) * 2016-05-10 2016-10-12 电子科技大学 一种基于sar图像的山区地表微形变信息提取方法
CN106526590A (zh) * 2016-11-04 2017-03-22 山东科技大学 一种融合多源sar影像工矿区三维地表形变监测及解算方法
CN107144213A (zh) * 2017-06-29 2017-09-08 中南大学 基于sar强度影像的矿区大量级三维时序形变估计方法及装置
CN108627832A (zh) * 2018-05-11 2018-10-09 电子科技大学 一种基于多时序sar图像提取输电通道地表形变的方法
CN108983232A (zh) * 2018-06-07 2018-12-11 中南大学 一种基于邻轨数据的InSAR二维地表形变监测方法
CN110058234A (zh) * 2019-05-20 2019-07-26 太原理工大学 一种解算矿区地表沉降三维形变的方法
CN110058237A (zh) * 2019-05-22 2019-07-26 中南大学 面向高分辨率SAR影像的InSAR点云融合及三维形变监测方法
CN110161501A (zh) * 2019-05-24 2019-08-23 电子科技大学 一种多时序sar图像的目标区域地表起伏信息提取方法
CN110888130A (zh) * 2019-10-30 2020-03-17 华东师范大学 一种基于升降轨时序InSAR的煤矿区地表形变监测方法
CN111474544A (zh) * 2020-03-04 2020-07-31 广东明源勘测设计有限公司 一种基于sar数据的滑坡形变监测及预警方法
CN111458709A (zh) * 2020-06-08 2020-07-28 河南大学 一种星载雷达广域地表二维形变场监测方法及装置
CN111650579A (zh) * 2020-06-12 2020-09-11 中南大学 一种岩移参数自适应获取的InSAR矿区三维形变估计方法、装置及介质
CN112068136A (zh) * 2020-09-14 2020-12-11 广东省核工业地质局测绘院 一种基于幅度偏移量的方位向形变监测方法
CN113091599A (zh) * 2021-04-06 2021-07-09 中国矿业大学 融合无人机dom和星载sar影像的地表三维形变提取方法
CN112904337A (zh) * 2021-05-07 2021-06-04 北京东方至远科技股份有限公司 一种基于Offset Tracking技术的边坡形变时序监测方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
马艳鸽 等: ""函数模型与随机模型双约束的GPS-INSAR数据融合方法建立三维形变场"", 《大地测量与地球动力学》 *

Also Published As

Publication number Publication date
CN114236541B (zh) 2023-05-16

Similar Documents

Publication Publication Date Title
Henriksen et al. Extracting accurate and precise topography from LROC narrow angle camera stereo observations
CN110111274B (zh) 一种星载推扫式光学传感器外方位元素定标方法
CN102968631B (zh) 山区多光谱遥感卫星影像的自动几何纠正与正射校正方法
CN108983232B (zh) 一种基于邻轨数据的InSAR二维地表形变监测方法
CN111458709B (zh) 一种星载雷达广域地表二维形变场监测方法及装置
CN102735216B (zh) Ccd立体相机三线阵影像数据平差处理方法
CN113358091B (zh) 一种利用三线阵立体卫星影像生产数字高程模型dem的方法
CN112284332B (zh) 基于高分辨率insar的高层建筑沉降监测结果三维定位方法
CN107705329A (zh) 基于几何约束的高分辨率光学卫星凝视影像配准方法
CN114791273B (zh) 一种针对滑坡的InSAR形变监测结果解释方法
CN109471104B (zh) 一种从两平行轨道sar数据中获取地表三维移动量的方法
CN109633639A (zh) Topsar干涉数据的高精度快速配准方法
Fraser et al. Precise georefrencing of long strips of ALOS imagery
Zhang et al. SAR mapping technology and its application in difficulty terrain area
CN115097450A (zh) 跨轨道高分三号sar偏移量大梯度滑坡形变估计方法
CN113238228B (zh) 基于水准约束的三维地表形变获取方法、系统及装置
CN108444451B (zh) 一种行星表面影像匹配方法与装置
CN112505696B (zh) 一种基于星载sar的露天矿边坡形变动态监测方法
CN116245757B (zh) 多模态数据的多场景通用性遥感影像云修复方法和系统
CN114236541A (zh) 基于Sentinel-1卫星SAR图像的大面积地表三维形变计算方法
CN116045920A (zh) 一种基于gf7号立体像对提取dem的方法
CN116299453A (zh) 星载sar非沿迹模式干涉高程测量方法
CN113532377B (zh) 一种利用高分七号激光测高数据辅助区域网平差的方法
Kirk et al. Ultrahigh resolution topographic mapping of Mars with HiRISE stereo images: Methods and first results
CN102735225B (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
GR01 Patent grant
GR01 Patent grant