CN108919262A - Dem辅助强度相关的冰川表面运动三维矢量反演方法 - Google Patents

Dem辅助强度相关的冰川表面运动三维矢量反演方法 Download PDF

Info

Publication number
CN108919262A
CN108919262A CN201810391224.8A CN201810391224A CN108919262A CN 108919262 A CN108919262 A CN 108919262A CN 201810391224 A CN201810391224 A CN 201810391224A CN 108919262 A CN108919262 A CN 108919262A
Authority
CN
China
Prior art keywords
displacement
glacier
sar
dem
satellite
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
CN201810391224.8A
Other languages
English (en)
Other versions
CN108919262B (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.)
China Aero Geophysical Survey & Remote Sensing Center For Land And Resources
Original Assignee
China Aero Geophysical Survey & Remote Sensing Center For Land And Resources
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 China Aero Geophysical Survey & Remote Sensing Center For Land And Resources filed Critical China Aero Geophysical Survey & Remote Sensing Center For Land And Resources
Priority to CN201810391224.8A priority Critical patent/CN108919262B/zh
Publication of CN108919262A publication Critical patent/CN108919262A/zh
Application granted granted Critical
Publication of CN108919262B publication Critical patent/CN108919262B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

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
    • 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
    • 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/9029—SAR image post-processing techniques specially adapted for moving target detection within a single SAR image or within multiple SAR images taken at the same time
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
    • G06T17/05—Geographic models

Landscapes

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

Abstract

本发明公开了一种DEM辅助强度相关的冰川表面运动三维矢量反演方法，涉及冰川运动监测技术领域。通过基于高精度DEM分别将升降轨SAR影像对高精度配准，然后利用星载SAR幅度信息应用偏移量跟踪技术获得了冰川区四个不同方向的位移观测结果，掩膜阴影、叠掩和低于偏移量跟踪相关性阈值的区域，并依据SAR影像成像几何关系建立起四个方向形变观测结果与东西、南北和垂直方向三个冰川三维运动分量之间的函数方程组，并计算该三维运动分量，进而得出冰川表面真实流动速率和流动方向。实现了仅利用星载SAR反演冰川表面真实三维运动，并在地形陡峭以及SAR影像空间基线较大时降低了偏移量误差，提高了偏移量跟踪结果的精度。

Description

DEM辅助强度相关的冰川表面运动三维矢量反演方法
技术领域
本发明涉及冰川运动监测技术领域，尤其涉及一种DEM辅助强度相关的冰川表面运动三维矢量反演方法。
背景技术
山地冰川是冰川的重要组成部分，是反映全球气候变化的敏感的记录器和指示器。冰川运动是冰川的重要特征之一，能够引起泥石流、冰川湖溃决等严重威胁人类生产生活的自然灾害，反演冰川表面真实流动状况对预测冰川流动和灾害发生有重要的参考意义。但是山地冰川一般位于自然环境恶劣的区域，难以进行规律地、大面积的实地监测，遥感的出现为低成本、高频率、大面积监测冰川变化提供了可能。
合成孔径雷达作为冰川研究的重要工具，可以全天时全天候收取SAR影像，利用SAR影像相位信息的DInSAR技术、MAI技术和利用SAR影像幅度信息的偏移量跟踪技术是如今测量冰川表面移动速度的三种重要技术。DInSAR技术和MAI技术可以获得厘米级甚至亚厘米级精度的冰川表面位移，但是仅可得到LOS向的一维形变或方位向的位移，且两者均受限于SAR影像之间的相干性。偏移量跟踪技术利用SAR影像幅度信息的配准偏移量来估计冰川表面的LOS向和方位向的二维形变场，能够更好的抵抗相位失相关的影响，每个方向形变监测精度约为所使用SAR影像分辨率的1/10到1/30，以3m分辨率的COSMO-SkyMed数据为例，其结果精度可以达到10cm甚至更高。受限于雷达数据的限制，早期一些学者提出了假设冰川流动方向平行于冰川表面、结合DEM反演出冰川三维运动场的方法，取得了一些成果，但是该方法是在假设的基础上建立起来，计算出的冰川三维运动场和冰川真实运动仍有差别。且由于山地冰川常位于地形起伏剧烈的区域,在SAR影像对空间基线较大时,应用传统偏移量跟踪技术获取的冰川LOS向的精度往往较低。
随着更高分辨率的雷达卫星尤其是卫星星座的发射，在同一时段(相隔不超过两天)获取同一地区的升降轨SAR数据已经成为可能，利用至少三个方向的SAR观测结果计算冰川运动的真实情况也得以实现。现有的山地冰川的三维流速反演多集中于融合InSAR技术和偏移量跟踪技术的结果，而冰川区域易失相干给InSAR技术的实施带来了限制。
发明内容
本发明的目的在于提供一种DEM辅助强度相关的冰川表面运动三维矢量反演方法，从而解决现有技术中存在的前述问题。
为了实现上述目的，本发明采用的技术方案如下：
一种DEM辅助强度相关的冰川表面运动三维矢量反演方法，包括如下步骤：
S1，获取高精度高分辨率的DEM数据和升降轨星载SAR数据，其中，DEM数据的高精度是指水平精度和高程精度优于现有的免费的30米分辨率的SRTM，DEM数据的高分辨率是指分辨率优于现有的免费的30米分辨率的SRTM，升降轨星载SAR数据包含一对升轨主辅SAR影像和一对降轨主辅SAR影像，且覆盖同一区域、时间周期分布类似；
S2，利用高精度高分辨率的DEM数据，辅助升降轨星载SAR数据进行精确配准，得到精确配准的升降轨星载SAR数据；
S3，基于精确配准的升降轨星载SAR数据的幅度信息，利用偏移量跟踪技术获取冰川升降轨雷达坐标系下的方位向位移和视线向位移四个不同方向的位移观测结果，并将得到的位移观测结果地理编码到地图坐标系下；
S4，将升降轨叠掩阴影区和升降轨四个不同方向的位移观测结果共同覆盖范围之外的区域进行掩膜，得到掩膜区；
S5，对掩膜区之外区域的有效像元，根据SAR卫星成像几何关系，分别建立四个不同方向的位移观测结果与冰川表面在东西向、南北向和垂向位移分量之间的函数关系，并利用最小二乘法计算出冰川表面在东西向、南北向和垂向的位移分量；
S6，合成东西向和南北向的位移分量，得到平面上冰川表面的位移分量；合成平面上冰川表面的位移分量与垂向的位移分量，得到冰川表面真实位移大小和方向。
优选地，S1中，所述DEM数据由资源三号立体像对、TerraSAR-X/TanDEM-X双星系统或WorldView-3立体像对提供。
优选地，S1中，所述升降轨星载SAR数据为单一SAR卫星数据或不同SAR卫星的联合数据，所述升降轨星载SAR数据由COSMO-SkyMed雷达星座或TerraSAR-X/TanDEM-X雷达卫星提供。
优选地，S2包括如下步骤：
S201，基于距离-多普勒方程建立起DEM像元与主辅SAR影像中像元之间的对应关系，计算得到主辅SAR影像对应像元的配准偏移量；
S202，基于S201得到的配准偏移量，将主辅SAR影像采样到同一雷达坐标系统中；
S203，在同一雷达坐标系统中，基于影像强度互相关方法获取主影像和辅影像的配准多项式，并依据多项式对辅影像进行重采样，实现主辅SAR影像的精确配准。
优选地，S201包括如下步骤：
S2011，SAR影像像元X(l,)用方位向坐标l和距离向坐标p描述，则用下式
其中：S和Vs分别为成像时刻传感器的位置以及速度矢量，c为光速，τ为单程距离向时间，Re为地球赤道半径，h为目标相对于椭球的高度，Rp＝(1-f)(Re+h)，f为地球的扁率；
计算地面目标点P(x,y,z)的方位向和距离向成像时间，得到地面目标点P在SAR影像上的位置(l,p)，即得到地面目标点P与SAR影像中目标点X(l,p)之间的对应关系；
S2012，用高精度DEM的三维信息表示地面目标点P，依据S2011的方法得到DEM像元P0与主辅SAR影像对应像元Xm和Xs之间的对应关系，
S2013，基于S2012得到的对应关系，计算得到主辅SAR影像对应像元Xm和Xs的配准偏移量。
优选地，S3中，所述升降轨雷达坐标系下的方位向位移和视线向位移四个不同方向的位移观测结果，包括升轨视线向形变值升轨方位向形变值降轨视线向形变值和降轨方位向形变值
优选地，S5中，所述函数关系为：
其中，为升轨视线向形变值，为升轨方位向形变值，为降轨视线向形变值，为降轨方位向形变值；
θa为升轨参考影像每个像元的雷达波入射角，αa为升轨参考影像的卫星轨道方位角，θd为降轨参考影像每个像元的雷达波入射角，αd为降轨参考影像的卫星轨道方位角；
dE为冰川表面真实位移在东向的位移矢量，dN为冰川表面真实位移在北向的位移矢量，dU为冰川表面真实位移在上向的位移矢量；
所述利用最小二乘法计算出冰川表面的三维运动分量，具体为：
令
X＝[dUdNdE]T
则有：
D＝CX
则利用最小二乘法求解出冰川表面流动速度在东向、北向、上向的位移分量。
优选地，S6具体为：
根据以下公式计算冰川表面真实位移大小和表示出冰川表面真实位移方向：
其中，
dS为冰川表面真实位移大小；
∠A为冰川表面真实位移方向在平面上的投影与正北方向的夹角，以北向角度为0，顺时针旋转角度增大，角度范围为0到360°；
∠B为冰川表面真实位移方向与水平面之间的夹角，水平面以下夹角为负，水平面以上夹角为正，位移方向在水平面上时夹角为0，角度范围-90到+90°。
本发明的有益效果是：本发明实施例提供的DEM辅助强度相关的冰川表面运动三维矢量反演方法，通过选取覆盖研究区冰川的高精度DEM和具有相似时间周期的升降轨SAR影像，并基于高精度DEM分别将升降轨SAR影像对高精度配准，然后利用星载SAR幅度信息应用偏移量跟踪技术获得了冰川区升轨视线向、升轨方位向、降轨视线向、降轨方位向四个不同方向的位移观测结果，掩膜阴影、叠掩和低于偏移量跟踪相关性阈值的区域，保留可靠的四个方向的观测结果的共同区域，依据SAR影像成像几何关系建立起四个方向形变观测结果与东西、南北和垂直方向三个冰川三维运动分量之间的函数方程组，利用最小二乘法解算出冰川表面在东西、南北和垂向上的三维运动分量，进而得出冰川表面真实流动速率和流动方向。本发明实现了仅利用星载SAR幅度信息反演出冰川表面真实三维运动，有效的克服了以往需要假设冰川表面流动方向平行于冰川表面或融合InSAR结果与偏移量跟踪结果的反演方法易受冰川大梯度位移、失相干影响的限制；并且该技术在地形陡峭以及SAR影像空间基线较大时能有效降低地形引起的偏移量误差，提高影像配准精度，进而提高了偏移量跟踪技术结果的精度。此外，该技术要求升降轨星载SAR数据具有相似的时间周期，使得监测结果更加真实；使用高精度的DEM作为辅助，大大提高了SAR影像间的配准精度、地理编码精度和三维解算所需要的局部入射角的精度，降低了三维求解结果的误差。该技术能监测出冰川表面的真实运动分布，对预测冰川运动趋势有重要的意义，同时，也为分析冰川表面运动机理、检验有关冰川流动理论提供可靠的数据支持。
附图说明
图1是本发明提供的反演方法流程示意图；
图2是地理编码后的升轨星载SAR偏移量跟踪方位向位移结果图，m/d；
图3是地理编码后的降轨星载SAR偏移量跟踪方位向位移结果图，m/d；
图4是地理编码后的升轨星载SAR偏移量跟踪视线向位移结果图，m/d；
图5是地理编码后的降轨星载SAR偏移量跟踪视线向位移结果图，m/d；
图6是降轨星载SAR影像对的阴影叠掩图；
图7是不可靠区域的掩膜范围，其中，黑色代表掩膜范围；
图8是SAR卫星成像几何关系示意图；
图9是应用本发明提供的方法得到的冰川表面真实位移在东向的运动速率图；
图10是应用本发明提供的方法得到的冰川表面真实位移在北向的运动速率图；
图11是应用本发明提供的方法得到的冰川表面真实位移在上向的运动速率图；
图12是冰川表面平面上的运动方向和三维真实流速大小分布示意图；
图13是冰川表面真实流动方向与水平面之间的夹角。
具体实施方式
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施方式仅仅用以解释本发明，并不用于限定本发明。
如图1所示，本发明实施例提供了一种DEM辅助强度相关的冰川表面运动三维矢量反演方法，包括如下步骤：
S1，获取高精度高分辨率的DEM数据和升降轨星载SAR数据，其中，DEM数据的高精度是指水平精度和高程精度优于现有的免费的30米分辨率的SRTM，DEM数据的高分辨率是指分辨率优于现有的免费的30米分辨率的SRTM，升降轨星载SAR数据包含一对升轨主辅SAR影像和一对降轨主辅SAR影像，且覆盖同一区域、时间周期分布类似；
S2，利用高精度高分辨率的DEM数据，辅助升降轨星载SAR数据进行精确配准，得到精确配准的升降轨星载SAR数据；
S3，基于精确配准的升降轨星载SAR数据的幅度信息，利用偏移量跟踪技术获取冰川升降轨雷达坐标系下的方位向位移和视线向位移四个不同方向的位移观测结果，并将得到的位移观测结果地理编码到地图坐标系下；
S4，将升降轨叠掩阴影区和升降轨四个不同方向的位移观测结果共同覆盖范围之外的区域进行掩膜，得到掩膜区；
S5，对掩膜区之外区域的有效像元，根据SAR卫星成像几何关系，分别建立四个不同方向的位移观测结果与冰川表面在东西向、南北向和垂向位移分量之间的函数关系，并利用最小二乘法计算出冰川表面在东西向、南北向和垂向的位移分量；
S6，合成东西向和南北向的位移分量，得到平面上冰川表面的位移分量；合成平面上冰川表面的位移分量与垂向的位移分量，得到冰川表面真实位移大小和方向。
上述方法中，基于精确配准的星载SAR数据的幅度信息，利用偏移量跟踪技术能直接得到LOS向和方位向位移，结合升降轨数据得到的四个方向位移观测结果，从而直接反演出冰川表面的真实流动结果，虽不如DInSAR技术和MAI技术的精度高，但在SAR影像之间不能保持良好相干性的情况下，利用偏移量跟踪技术可以获得形变结果，具有较好的鲁棒性。
另外，本发明中，偏移量跟踪技术是基于高精度DEM辅助配准的星载SAR数据的幅度信息，所以，对于山地冰川地形陡峭以及SAR数据对空间基线较大的情况下，能够获得可靠的结果，减小了地形以及空间基线长度给偏移量跟踪技术结果带来的限制，有利于更好开展山地冰川表面运动监测。
所以，采用本发明实施例提供的方法，利用高分辨率SAR数据，配合高精度DEM，利用偏移量跟踪技术能够在冰川面积较小，地势较为复杂的山区获得可靠的冰川表面运动的三维矢量，最终得到冰川表面的真实流动结果。解决了现有技术中，需基于冰川表面流动方向平行于冰川表面的假设，造成计算结果与冰川真实运动有差别，不能很好的监测冰川运动的问题，以及需融合InSAR结果与偏移量跟踪结果的反演方法，造成易受山地冰川大梯度位移、失相干影响的问题。
其中，S1中，所述DEM数据可由资源三号立体像对、TerraSAR-X/TanDEM-X双星系统或WorldView-3立体像对提供。
S1中，所述升降轨星载SAR数据可以为单一SAR卫星数据或不同SAR卫星的联合数据，所述升降轨星载SAR数据可由COSMO-SkyMed雷达星座或TerraSAR-X/TanDEM-X雷达卫星提供。
本实施例中，高精度DEM要求为精度和分辨率高于现有的免费的SRTM，资源三号立体像对、TerraSAR-X/TanDEM-X双星系统、WorldView-3立体像对等可为研究区提供高分辨率高精度的DEM数据；星载SAR数据要求包含一对升轨SAR重轨影像和一对降轨SAR重轨影像，可为单一SAR卫星数据或联合不同SAR卫星数据，影像覆盖范围均应包含冰川区，且升降轨SAR影像的时间周期应尽量的一致，以保证时间段内冰川具有相同的运动速度，以COSMO-SkyMed雷达星座和TerraSAR-X/TanDEM-X雷达卫星为代表的高分辨率雷达卫星系统可为获取覆盖同一区域、时间周期分布类似的升降轨星载SAR数据提供数据源。
本发明的一个优选实施例中，S2可以包括如下步骤：
S201，基于距离-多普勒方程建立起DEM像元与主辅SAR影像中像元之间的对应关系，计算得到主辅SAR影像对应像元的配准偏移量；
S202，基于S201得到的配准偏移量，将主辅SAR影像采样到同一雷达坐标系统中；
S203，在同一雷达坐标系统中，基于影像强度互相关方法获取主影像和辅影像的配准多项式，并依据多项式对辅影像进行重采样，实现主辅SAR影像的精确配准。
其中，S201可以包括如下步骤：
S2011，用方位向坐标l和距离向坐标p描述SAR影像像元X(l,p)，则用下式计算地面目标点P(x,y,z)的方位向和距离向成像时间，得到地面目标点P在SAR影像上的位置(l,p)，即得到地面目标点P与SAR影像中目标点X(l,p)之间的对应关系：
其中：S和Vs分别为成像时刻传感器的位置以及速度矢量，c为光速，τ为单程距离向时间，Re为地球赤道半径，h为目标相对于椭球的高度，Rp＝(1-f)(Re+h)，f为地球的扁率；
S2012，用高精度DEM的三维信息表示地面目标点P，依据S2011的方法得到DEM像元P0与主辅SAR影像对应像元Xm和Xs之间的对应关系，
S2013，基于S2012得到的对应关系，计算得到主辅SAR影像对应像元Xm和Xs的配准偏移量。
上述方法可使用现有的基于lookup table的SAR影像配准方法实现。
本发明实施例中，S3中，所述升降轨雷达坐标系下的方位向位移和视线向位移四个不同方向的位移观测结果，包括升轨视线向形变值升轨方位向形变值降轨视线向形变值和降轨方位向形变值
本实施例中，基于对已配准星载SAR数据的幅度信息，应用偏移量跟踪技术可获取其雷达坐标系下的方位向形变和视线向形变，则针对覆盖同一区域的升降轨两对SAR数据，可获得研究区四个不同方向的位移观测结果，将雷达坐标系下的结果地理编码到同一地图坐标系下，比如WGS84坐标系，则对一般的地面目标点，都包括升轨视线向形变值升轨方位向形变值降轨视线向形变值降轨方位向形变值四个方向的形变结果。
本发明实施例的S4中，由于升轨SAR和降轨SAR数据卫星轨道方位角的差异，导致SAR数据的阴影和叠掩的分布有所不同，所以，将升轨和降轨SAR偏移量跟踪结果的叠掩和阴影区进行掩膜。由于偏移量跟踪技术仅针对高于相关性阈值的区域进行计算，升轨SAR和降轨SAR的位移值区域也有所差别，所以，将两者低于阈值的区域取并集，然后进行掩膜。
本发明实施例的S5中，冰川表面真实位移可以分解为东向dE、北向dN和上向dU三个方向的位移矢量，其中，视线向位移是上述三个位移矢量在卫星视线向的投影之和，方位向位移是上述三个位移矢量在卫星方位向的投影之和，所以，根据SAR卫星的成像几何关系，可以得到所述函数关系为：
其中，为升轨视线向形变值，为升轨方位向形变值，为降轨视线向形变值，为降轨方位向形变值；
θa为升轨参考影像每个像元的雷达波入射角，αa为升轨参考影像的卫星轨道方位角，θd为降轨参考影像每个像元的雷达波入射角，αd为降轨参考影像的卫星轨道方位角；
dE为冰川表面真实位移在东向的位移矢量，dN为冰川表面真实位移在北向的位移矢量，dU为冰川表面真实位移在上向的位移矢量；
对于形变区域内的每个有效像元，都可以构建出以上四个方程；
令
X＝[dUdNdE]T
则：
D＝CX
则：可以利用最小二乘法求解出冰川表面流动速度在东向、北向、上向的位移分量。
在本发明的一个优选实施例中，S6具体可以为：
根据以下公式计算冰川表面真实位移大小和表示出冰川表面真实位移方向：
其中，
dS为冰川表面真实位移大小；
∠A为冰川表面真实位移方向在平面上的投影与正北方向的夹角，以北向角度为0，顺时针旋转角度增大，角度范围为0到360°；
∠B为冰川表面真实位移方向与水平面之间的夹角，水平面以下夹角为负，水平面以上夹角为正，位移方向在水平面上时夹角为0，角度范围-90到+90°。
具体实施例
以估算西藏那曲地区嘉黎县依嘎冰川表面运动三维矢量为例，说明本发明在实际应用时的具体操作步骤。如图1所示，本发明实施例可以按照如下具体步骤进行实施：
步骤1：选取高精度DEM和升降轨星载SAR数据。
COSMO-SkyMed雷达卫星是由4颗X波段的高分辨率雷达卫星组成的星座，可快速获取同一地区时间分布类似的升降轨SAR影像对，选取覆盖依嘎冰川的升降轨COSMO-SkyMed数据，为HImage成像模式获取的3m分辨率的影像，极化方式为HH，升降轨数据时间差仅有1天12个小时11分钟，可认为两个时间周期内冰川表面速度一致，具体参数详见表1。高精度DEM选取由TerraSAR-X/TanDEM-X生成的分辨率为5m的WorldDEMTM高程结果。
表1选用COSMO-SkyMed影像主要参数表
步骤2：基于高精度DEM辅助的星载SAR数据的精确配准，基于已配准星载SAR数据的幅度信息，利用偏移量跟踪技术获取冰川升降轨雷达坐标系下的方位向位移和视线向位移，并将结果地理编码到地图坐标系下。
具体实现过程可以为：基于距离-多普勒方程建立起DEM像元与主辅雷达图像中像元之间的对应关系，求出主辅影像对应像元的配准偏移量；基于此配准偏移量将主辅影像采样到同一雷达坐标系统中，进一步使用影像强度互相关方法获取主影像和重采样从影像的配准多项式，并依据多项式对辅影像进行重采样，从而实现主辅SAR影像精确配准。基于已配准的升降轨COSMO-SkyMed数据应用偏移量跟踪技术获得依嘎冰川升降轨雷达坐标系下的方位向位移和视线向位移共四个不同方向的位移观测结果，并将结果地理编码到WGS84坐标系下，对于冰川上每一个像元都有升轨视线向形变值(如图2所示)、升轨方位向形变值(如图3所示)、降轨视线向形变值(如图4所示)、降轨方位向形变值(如图5所示)四个方向的形变结果，并裁剪出包含依嘎冰川的共同区域。
步骤3：掩膜升降轨叠掩阴影区和升降轨四个位移结果共同覆盖范围之外的区域。
基于高精度DEM和COSMO-SkyMed影像的参数信息分别生成共同区域升轨和降轨WGS84坐标系下的阴影叠掩图(如图6所示)，将升轨和降轨的阴影叠掩范围取并集进行掩膜。应用偏移量跟踪技术时，需设定一个阈值(0.1)，仅高于阈值的区域有位移结果，因此升轨SAR影像对和降轨SAR影像对的方位向和视线向结果有不同的位移范围，将两对数据低于阈值(0.1)的区域取并集进行掩膜。将叠掩阴影掩膜范围和低于阈值的掩膜范围进行合并，得到最终的掩膜图(如图7所示)，除掩膜区之外的共同区域进行后续计算。
步骤4：根据SAR成像几何关系(如图8所示)建立四个不同方向的偏移量跟踪结果和冰川表面在东西向、南北向和垂向上的三维运动分量的函数关系，并利用最小二乘法解算出冰川表面的三维运动分量。
具体的实施过程可以为：将依嘎冰川四个方向的观测结果依据步骤3的掩膜图范围进行掩膜。将依嘎冰川表面真实位移分解为东dE、北dN和垂向dU上三个方向的位移矢量，可知视线向位移是三个位移矢量在卫星视线向的投影之和，方位向位移为三个位移矢量在卫星方位向的投影之和，根据SAR卫星的成像几何关系(如图8所示)可以得到：
其中，θa为升轨参考影像每个像元的雷达波入射角，αa为升轨参考影像的卫星轨道方位角，θd为降轨参考影像每个像元的雷达波入射角，αd为降轨参考影像的卫星轨道方位角。θa和θd可根据高精度DEM及卫星参数求解出，αa和αd可从卫星参数信息中查找到，由此对于冰川上的每一个有效像元，都可以构建出以上四个方程。
令
X＝[dUdNdE]T
则有
D＝CX
利用最小二乘法就可以求解出依嘎冰川表面流动速度在东、北、上三个方向的速度分量，见图9-11所示。
步骤5：合成东西向和南北向的位移分量，计算出平面上冰川表面位移大小和方向，合成冰川表面的平面位移与垂向位移，得出冰川表面真实运动速率和真实流动方向。
具体可以按照如下方法实施：依嘎冰川表面真实流动方向用其流动方向在平面上的投影与正北方向的夹角∠A(以北向角度为0，顺时针旋转，角度增大，角度范围为0到360度，结果见图12所示)和真实流动方向与水平面之间的夹角∠B(水平面以下夹角为负，水平面以上夹角为正，流动方向在水平面上时夹角为0，角度范围-90到+90度，结果见图13所示)来表示。冰川表面真实流速大小用ds来表示。
根据以上公式，可求出真实流速大小和表示出冰川表面的真实流动方向。
经过验证，东西、南北和垂向上速度分量的均方根误差分别为0.83、1.58和0.33cm/d，合成速度均方根误差约为1.8cm/d，远低于冰川表面的日均运动速度，所以，本发明实施例提供的反演方法是可行的。
通过采用本发明公开的上述技术方案，得到了如下有益的效果：本发明实施例提供的DEM辅助强度相关的冰川表面运动三维矢量反演方法，通过选取覆盖研究区冰川的高精度DEM和具有相似时间周期的升降轨SAR影像，并基于高精度DEM分别将升降轨SAR影像对高精度配准，然后利用星载SAR幅度信息应用偏移量跟踪技术获得了冰川区升轨视线向、升轨方位向、降轨视线向、降轨方位向四个不同方向的位移观测结果，掩膜阴影、叠掩和低于偏移量跟踪相关性阈值的区域，保留可靠的四个方向的观测结果的共同区域，依据SAR影像成像几何关系建立起四个方向形变观测结果与东西、南北和垂直方向三个冰川三维运动分量之间的函数方程组，利用最小二乘法解算出冰川表面在东西、南北和垂向上的三维运动分量，进而得出冰川表面真实流动速率和流动方向。本发明实现了仅利用星载SAR幅度信息反演出冰川表面真实三维运动，有效的克服了以往需要假设冰川表面流动方向平行于冰川表面或融合InSAR结果与偏移量跟踪结果的反演方法易受冰川大梯度位移、失相干影响的限制；并且该技术在地形陡峭以及SAR影像空间基线较大时能有效降低地形引起的偏移量误差，提高影像配准精度，进而提高了偏移量跟踪技术结果的精度。此外，该技术要求升降轨星载SAR数据具有相似的时间周期，使得监测结果更加真实；使用高精度的DEM作为辅助，大大提高了SAR影像间的配准精度、地理编码精度和三维解算所需要的局部入射角的精度，降低了三维求解结果的误差。该技术能监测出冰川表面的真实运动分布，对预测冰川运动趋势有重要的意义，同时，也为分析冰川表面运动机理、检验有关冰川流动理论提供可靠的数据支持。
以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视本发明的保护范围。

Claims (8)

1.一种DEM辅助强度相关的冰川表面运动三维矢量反演方法，其特征在于，包括如下步骤：
S1，获取高精度高分辨率的DEM数据和升降轨星载SAR数据，其中，DEM数据的高精度是指水平精度和高程精度优于现有的免费的30米分辨率的SRTM，DEM数据的高分辨率是指分辨率优于现有的免费的30米分辨率的SRTM，升降轨星载SAR数据包含一对升轨主辅SAR影像和一对降轨主辅SAR影像，且覆盖同一区域、时间周期分布类似；
S2，利用高精度高分辨率的DEM数据，辅助升降轨星载SAR数据进行精确配准，得到精确配准的升降轨星载SAR数据；
S3，基于精确配准的升降轨星载SAR数据的幅度信息，利用偏移量跟踪技术获取冰川升降轨雷达坐标系下的方位向位移和视线向位移四个不同方向的位移观测结果，并将得到的位移观测结果地理编码到地图坐标系下；
S4，将升降轨叠掩阴影区和升降轨四个不同方向的位移观测结果共同覆盖范围之外的区域进行掩膜，得到掩膜区；
S5，对掩膜区之外区域的有效像元，根据SAR卫星成像几何关系，分别建立四个不同方向的位移观测结果与冰川表面在东西向、南北向和垂向位移分量之间的函数关系，并利用最小二乘法计算出冰川表面在东西向、南北向和垂向的位移分量；
S6，合成东西向和南北向的位移分量，得到平面上冰川表面的位移分量；合成平面上冰川表面的位移分量与垂向的位移分量，得到冰川表面真实位移大小和方向。
2.根据权利要求1所述的DEM辅助强度相关的冰川表面运动三维矢量反演方法，其特征在于，S1中，所述DEM数据由资源三号立体像对、TerraSAR-X/TanDEM-X双星系统或WorldView-3立体像对提供。
3.根据权利要求1所述的DEM辅助强度相关的冰川表面运动三维矢量反演方法，其特征在于，S1中，所述升降轨星载SAR数据为单一SAR卫星数据或不同SAR卫星的联合数据，所述升降轨星载SAR数据由COSMO-SkyMed雷达星座或TerraSAR-X/TanDEM-X雷达卫星提供。
4.根据权利要求1所述的DEM辅助强度相关的冰川表面运动三维矢量反演方法，其特征在于，S2包括如下步骤：
S201，基于距离-多普勒方程建立起DEM像元与主辅SAR影像中像元之间的对应关系，计算得到主辅SAR影像对应像元的配准偏移量；
S202，基于S201得到的配准偏移量，将主辅SAR影像采样到同一雷达坐标系统中；
S203，在同一雷达坐标系统中，基于影像强度互相关方法获取主影像和辅影像的配准多项式，并依据多项式对辅影像进行重采样，实现主辅SAR影像的精确配准。
5.根据权利要求4所述的DEM辅助强度相关的冰川表面运动三维矢量反演方法，其特征在于，S201包括如下步骤：
S2011，SAR影像像元X(l,p)用方位向坐标l和距离向坐标p描述，则用下式
其中：S和Vs分别为成像时刻传感器的位置以及速度矢量，c为光速，τ为单程距离向时间，Re为地球赤道半径，h为目标相对于椭球的高度，Rp＝(1-f)(Re+h)，f为地球的扁率；
计算地面目标点P(x,y,z)的方位向和距离向成像时间，得到地面目标点P在SAR影像上的位置(l,p)，即得到地面目标点P与SAR影像中目标点X(l,p)之间的对应关系；
S2012，用高精度DEM的三维信息表示地面目标点P，依据S2011的方法得到DEM像元P0与主辅SAR影像对应像元Xm和Xs之间的对应关系，
S2013，基于S2012得到的对应关系，计算得到主辅SAR影像对应像元Xm和Xs的配准偏移量。
6.根据权利要求1所述的DEM辅助强度相关的冰川表面运动三维矢量反演方法，其特征在于，S3中，所述升降轨雷达坐标系下的方位向位移和视线向位移四个不同方向的位移观测结果，包括升轨视线向形变值升轨方位向形变值降轨视线向形变值和降轨方位向形变值
7.根据权利要求6所述的DEM辅助强度相关的冰川表面运动三维矢量反演方法，其特征在于，S5中，所述函数关系为：
其中，为升轨视线向形变值，为升轨方位向形变值，为降轨视线向形变值，为降轨方位向形变值；
θa为升轨参考影像每个像元的雷达波入射角，αa为升轨参考影像的卫星轨道方位角，θd为降轨参考影像每个像元的雷达波入射角，αd为降轨参考影像的卫星轨道方位角；
dE为冰川表面真实位移在东向的位移矢量，dN为冰川表面真实位移在北向的位移矢量，dU为冰川表面真实位移在上向的位移矢量；
所述利用最小二乘法计算出冰川表面的三维运动分量，具体为：
令
X＝[dUdNdE]T
则有：
D＝CX
则利用最小二乘法求解出冰川表面流动速度在东向、北向、上向的位移分量。
8.根据权利要求7所述的DEM辅助强度相关的冰川表面运动三维矢量反演方法，其特征在于，S6具体为：
根据以下公式计算冰川表面真实位移大小和表示出冰川表面真实位移方向：
其中，
dS为冰川表面真实位移大小；
∠A为冰川表面真实位移方向在平面上的投影与正北方向的夹角，以北向角度为0，顺时针旋转角度增大，角度范围为0到360°；
∠B为冰川表面真实位移方向与水平面之间的夹角，水平面以下夹角为负，水平面以上夹角为正，位移方向在水平面上时夹角为0，角度范围-90到+90°。
CN201810391224.8A 2018-04-27 2018-04-27 Dem辅助强度相关的冰川表面运动三维矢量反演方法 Expired - Fee Related CN108919262B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810391224.8A CN108919262B (zh) 2018-04-27 2018-04-27 Dem辅助强度相关的冰川表面运动三维矢量反演方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810391224.8A CN108919262B (zh) 2018-04-27 2018-04-27 Dem辅助强度相关的冰川表面运动三维矢量反演方法

Publications (2)

Publication Number Publication Date
CN108919262A true CN108919262A (zh) 2018-11-30
CN108919262B CN108919262B (zh) 2019-04-16

Family

ID=64403793

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810391224.8A Expired - Fee Related CN108919262B (zh) 2018-04-27 2018-04-27 Dem辅助强度相关的冰川表面运动三维矢量反演方法

Country Status (1)

Country Link
CN (1) CN108919262B (zh)

Cited By (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109635713A (zh) * 2018-12-07 2019-04-16 云南大学 高原山地阴影区冰川识别方法
CN110017762A (zh) * 2019-03-01 2019-07-16 四川省地质工程勘察院 一种利用Offset-Tracking技术监测蠕变体形变的方法
CN110020404A (zh) * 2019-04-10 2019-07-16 自然资源部第二海洋研究所 一种角度约束的遥感反演流场的矢量数据处理方法
CN111310649A (zh) * 2020-02-13 2020-06-19 西南交通大学 一种无人机高分辨率影像对山岳冰川运动消融的提取方法
CN112179314A (zh) * 2020-09-25 2021-01-05 北京空间飞行器总体设计部 一种基于三维网格投影的多角度sar高程测量方法及系统
CN112395789A (zh) * 2020-10-23 2021-02-23 马培峰 一种耦合InSAR和数值模拟分析城区滑坡形变的方法
CN112729145A (zh) * 2020-12-21 2021-04-30 兰州大学 一种确定滑坡位移与库区水位波动关系的方法
RU2773277C2 (ru) * 2020-10-19 2022-06-01 Публичное акционерное общество "Ракетно-космическая корпорация "Энергия" имени С.П. Королёва" Способ определения с орбитального космического аппарата параметров движения объекта преимущественно смещающихся природных масс ледника и оползня
CN115291215A (zh) * 2022-09-29 2022-11-04 中国科学院空天信息创新研究院 基于升降轨sar卫星的长时序二维形变快速解算方法
CN115951354A (zh) * 2023-02-14 2023-04-11 中国铁道科学研究院集团有限公司铁道建筑研究所 一种融合升降轨的D-InSAR形变监测方法
CN116485857A (zh) * 2023-05-05 2023-07-25 中山大学 一种基于多源遥感数据的高时间分辨率冰川厚度反演方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105929398A (zh) * 2016-04-20 2016-09-07 中国电力工程顾问集团中南电力设计院有限公司 结合外控点的InSAR高精度高分辨率DEM获取方法
CN106066478A (zh) * 2016-05-27 2016-11-02 中国矿业大学 融合像元偏移跟踪和短基线集的矿区地表形变解算方法
CN106556834A (zh) * 2016-11-24 2017-04-05 首都师范大学 一种从两平行轨道sar数据集中精确提取地面垂直形变方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105929398A (zh) * 2016-04-20 2016-09-07 中国电力工程顾问集团中南电力设计院有限公司 结合外控点的InSAR高精度高分辨率DEM获取方法
CN106066478A (zh) * 2016-05-27 2016-11-02 中国矿业大学 融合像元偏移跟踪和短基线集的矿区地表形变解算方法
CN106556834A (zh) * 2016-11-24 2017-04-05 首都师范大学 一种从两平行轨道sar数据集中精确提取地面垂直形变方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
JUN LU: "Sentinel-1 Toolbox", 《OFFSET TRACKING TUTORIAL》 *
SHIYONG YAN ET AL.: "Mountain glacier displacement estimation using a DEM-assisted offset tracking method with ALOS/PALSAR data", 《REMOTE SENSING LETTERS》 *

Cited By (16)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109635713A (zh) * 2018-12-07 2019-04-16 云南大学 高原山地阴影区冰川识别方法
CN110017762A (zh) * 2019-03-01 2019-07-16 四川省地质工程勘察院 一种利用Offset-Tracking技术监测蠕变体形变的方法
CN110020404A (zh) * 2019-04-10 2019-07-16 自然资源部第二海洋研究所 一种角度约束的遥感反演流场的矢量数据处理方法
CN111310649A (zh) * 2020-02-13 2020-06-19 西南交通大学 一种无人机高分辨率影像对山岳冰川运动消融的提取方法
CN111310649B (zh) * 2020-02-13 2022-09-23 西南交通大学 一种无人机高分辨率影像对山岳冰川运动消融的提取方法
CN112179314A (zh) * 2020-09-25 2021-01-05 北京空间飞行器总体设计部 一种基于三维网格投影的多角度sar高程测量方法及系统
RU2773277C2 (ru) * 2020-10-19 2022-06-01 Публичное акционерное общество "Ракетно-космическая корпорация "Энергия" имени С.П. Королёва" Способ определения с орбитального космического аппарата параметров движения объекта преимущественно смещающихся природных масс ледника и оползня
CN112395789A (zh) * 2020-10-23 2021-02-23 马培峰 一种耦合InSAR和数值模拟分析城区滑坡形变的方法
CN112729145B (zh) * 2020-12-21 2022-03-11 兰州大学 一种确定滑坡位移与库区水位波动关系的方法
CN112729145A (zh) * 2020-12-21 2021-04-30 兰州大学 一种确定滑坡位移与库区水位波动关系的方法
CN115291215A (zh) * 2022-09-29 2022-11-04 中国科学院空天信息创新研究院 基于升降轨sar卫星的长时序二维形变快速解算方法
CN115291215B (zh) * 2022-09-29 2022-12-20 中国科学院空天信息创新研究院 基于升降轨sar卫星的长时序二维形变快速解算方法
CN115951354A (zh) * 2023-02-14 2023-04-11 中国铁道科学研究院集团有限公司铁道建筑研究所 一种融合升降轨的D-InSAR形变监测方法
CN115951354B (zh) * 2023-02-14 2023-06-06 中国铁道科学研究院集团有限公司铁道建筑研究所 一种融合升降轨的D-InSAR形变监测方法
CN116485857A (zh) * 2023-05-05 2023-07-25 中山大学 一种基于多源遥感数据的高时间分辨率冰川厚度反演方法
CN116485857B (zh) * 2023-05-05 2024-06-28 中山大学 一种基于多源遥感数据的高时间分辨率冰川厚度反演方法

Also Published As

Publication number Publication date
CN108919262B (zh) 2019-04-16

Similar Documents

Publication Publication Date Title
CN108919262B (zh) Dem辅助强度相关的冰川表面运动三维矢量反演方法
CN101770027B (zh) 基于InSAR与GPS数据融合的地表三维形变监测方法
Furgale et al. The Devon Island rover navigation dataset
Li Potential of high-resolution satellite imagery for national mapping products
Kääb et al. Coseismic displacements of the 14 November 2016 M w 7.8 Kaikoura, New Zealand, earthquake using the Planet optical cubesat constellation
CN109212522B (zh) 一种基于双基星载sar的高精度dem反演方法及装置
CN102968631B (zh) 山区多光谱遥感卫星影像的自动几何纠正与正射校正方法
CN102654576B (zh) 基于sar图像和dem数据的图像配准方法
US20160259044A1 (en) Three-dimensional positioning method
CN101876701B (zh) 一种侧视雷达遥感影像定位方法
CN111736152B (zh) 一种道路边坡稳定性监测方法及车载平台装置
CN113340191B (zh) 时间序列干涉sar的形变量测量方法及sar系统
Scheiber et al. Speckle tracking and interferometric processing of TerraSAR-X TOPS data for mapping nonstationary scenarios
CN108983232B (zh) 一种基于邻轨数据的InSAR二维地表形变监测方法
CN107144293A (zh) 一种视频卫星面阵相机的几何定标方法
CN102628942B (zh) 一种雷达影像双视向信息补偿方法
CN105352509A (zh) 地理信息时空约束下的无人机运动目标跟踪与定位方法
CN107918127A (zh) 一种基于车载InSAR的道路边坡形变检测系统及方法
CN102866393B (zh) 一种基于pos与dem数据的sar多普勒参数估计方法
Poli et al. SPOT-5/HRS stereo images orientation and automated DSM generation
CN110456346A (zh) 一种基于InSAR技术的输电铁塔倾斜监测方法
CN115469308B (zh) 多轨道InSAR震间形变速率场拼接方法、装置、设备及介质
De Oliveira et al. Assessment of radargrammetric DSMs from TerraSAR-X Stripmap images in a mountainous relief area of the Amazon region
Liu et al. Correction of positional errors and geometric distortions in topographic maps and DEMs using a rigorous SAR simulation technique
Recla et al. From Relative to Absolute Heights in SAR-based Single-Image Height Prediction

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
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20190416

CF01 Termination of patent right due to non-payment of annual fee