陆探一号(LT-1)InSAR处理方法

发布时间 2023-06-12 20:50:54作者: ENVI-IDL技术殿堂

L波段差分干涉SAR卫星(陆探一号)采用双星编队构型运行,具备双星绕飞与双星跟飞两种模式,利用干涉测高和差分形变测量技术,实现高精度、全天时、全天候地形测量、地表形变和地质灾害监测等任务。

SARscape5.6.2.1新出补丁初步支持陆探一号数据的导入,数据导入后就能进行后续的所有处理,包括InSAR处理。本文是在SARscape5.6.2.1下完成,包括数据导入、基线估算、干涉图生成、滤波及相干性计算、相位解缠、轨道精炼与重去平等。

由于本文档中使用的数据对成像间隔内没有明确的沉降现象发生,因此本文档主要说明InSAR/DInSAR各个处理步骤和结果分析。使用分步骤处理工具(/SARscape/Interferometry/Phase Processing),在实际应用中可使用效率更高的流程化处理工具(/SARscape/Interferometry/InSAR DEM Workflow,或者/SARscape/Interferometry/DInSAR Displacement Workflow)。

SARscape5.6下载与试用:https://envi.geoscene.cn/sar_license

补丁下载:https://www.sarmap.ch/5.6/2/1/preliminary/last.html  或者

https://pan.baidu.com/s/1Bapv7ucrfJQ2QXvark8JIw?pwd=envi

1数据导入

陆探一号搭载了先进的L波段多极化多通道SAR载荷,具备6种成像模式,最高分辨率3米,最大观测幅宽可达400公里。

如下表所示为陆探一号成像模式说明。

成像模式名称

分辨率/m

幅宽/km

极化方式

条带模式1

3

50

HH/VV

条带模式2

12

100

HH/VV

条带模式3

3

50

HH+HV/VV+VH

条带模式4

6

30

HH+HV+VH+VV

条带模式5

20~30

150~250

HH/VV

扫描模式

30

400

HH/VV

 

本文档使用的是陆探一号L1级产品,它是单视复数据(SLC),条带模式2(STRIP2)成像模式。

在处理之前,我们选择一套默认参数,打开Toolbox/SARscape/Preferences/Preferences specific,在打开的界面中,选择Load Preferences->VHR(better 10m and 30m),Cartographic Grid Size(m):12。

注:没有拿到详细的文件格式说明文档,这样对应不知道是否正确。如通过查询元数据文件(*meta.xml)中的<IncidenceAngle>:39.8度,距离向采样间隔< columnSpacing >:1.67米,最大制图分辨率为1.67/sin(39.8)=2.6。

(1)在Toolbox中,选择/SARscape/Import Data/SAR Spaceborne/Single Sensor/LuTan-1。

(2)在打开的面板中,

  • 数据输入面板(Input Files)
  • 输入文件(Input File List):输入过滤的.meta.xml文件
  • 参数设置面板(Parameters):主要参数(Principal Parameters)
  • 极化方式(Polarization):ALL,输出所有的极化数据,可以选择只输出同极化或者交叉极化的数据;
  • 对数据重命名(Rename the File Using Parameters):True。软件会自动在输入文件名的基础上增加几个标识字母,如增加“_HH_slc”。
  • 数据输出面板(Output Files)

输出文件(Output file list):自动读取ENVI默认的数据输出目录以及输入面板中的数据文件名。

注:1.如果要修改输出的路径,在右边单击文件夹图标选择输出文件夹目录。

    2、如果要修改输出的路径,在输出文件名右键选择Edit菜单。

(3)单击Exec按钮开始执行。

图:数据导入

生成的结果除了图像文件外,还包括Shapefile和KML格式的图像轮廓线。

2基线估算

对slc数据对进行基线估算,查看数据的基线情况,包括时间基线、空间基线等。

打开基线估算工具/SARscape/Interferometry/Interferometric Tools/Baseline Estimation。在Input Files面板中,Input Master File选择一个slc数据输入,Input Slave File选择另外一个slc数据输入,其他按照默认,点击Exec按钮,计算基线。

图:基线估算结果

该数据对基线估算结果如下:

  • 空间基线为114米(远小于临界基线23165米);
  • 时间基线为4天;
  • 高程的变化是53米(2 PI Ambiguity height (InSAR) (m));
  • 相位变化一个2π周期,代表着地表发生了118米的形变;
  • 立体像对1个像素移动模糊高程为422米;
  • 振幅跟踪1个像素移动模糊形变为666米;

该数据对比较适合DinSAR或者Amplitude Tracking处理。

3干涉图生成

(1)Toolbox中,选择/SARscape/Interferometry/Phase Processing/1 - Interferogram Generation,打开Interferometry Generation面板。

(2)在Interferometry Generation面板:

  • Input Files面板
    • Input Master File项,选择前一时相的SLC文件,这时会根据头文件中的信息和默认的制图分辨率自动计算出视数,点击确定。
    • Input Slave File项,选择后一时相的SLC文件
  • Optional Files面板:
    • Geometry GCP File:地理控制点文件,用于配准时候,这里默认不设置。
    • Coregistration File:配准文件,若有已有的配准文件可以在此输入,这里默认不设置。
    • Shift Parameters File:配准时生成的偏移参数文件,当Output Files面板中设置输出文件后,该文件会自动添加。这里默认不设置。
  • DEM/Cartographic System面板:输入DEM文件或投影信息。若是输入DEM数据,自动从输入的参考DEM中获取相应的信息,输出结果默认以DEM投影参数为准。如果不输入DEM数据,则需要设置Output Projection。这里点击,选择参考DEM文件是通过:/SARscape/Import Data/DEM Extraction/ALOS World 3D 30m 自动下载的DEM文件。
  • Parameters面板。主要参数(Principal Parameters):
    • Range Looks:4
    • Azimuth Looks:2
    • Grid Size for Suggested Looks:12
    • Compute Shift Parameters:True
    • Generate Coregistered SLC:False
    • Coregistration With DEM:True

图:Interferometry Generation面板

(3)Output Files面板,设置输出路径及根文件名,生成的相应结果都会自动在该文件名后添加后缀,如_dint为去平后的干涉图。

处理完成后,弹出Completed对话框,点击确定。在输出路径下生成以下结果:

–  _dint:去已知地形的干涉图(差分干涉图)

–  _int:干涉图

–  _master_pwr:主影像强度图

–  _slave_pwr:辅影像强度图

–  _par:文本文件保存的配准时的偏移量数据,

–  _sint:合成的相位

–  _srdem:斜距几何下的参考DEM

–  _par_orbit_off:估算轨道偏移时用到的点矢量数据,包括在方位向和距离向 的点的位置坐标、测量到的偏移量、计算出的线性偏移量。

–  _par_winCC_off:从强度数据的相差上估算交叉相关偏移量的点矢量数据。包含每个点上交叉相关的值(CC),范围是0-1。

–  _par_winCoh_off:从相位数据的相差上估算相干性的点矢量数据,包含信噪比(SNR)和相干性,相干性值的范围是0-1。

每个结果都相应地生成了TIF格式的快视图文件,便于直接查看。如下是_int和_dint的结果。

图:int结果

图:Dint结果

4 自适应滤波及相干性计算

对上一步去已知地形的干涉图(_dint)进行滤波,去掉由平地干涉引起的位相噪声。同时生成相干系数图(描述位相质量)和滤波后的主影像强度图。

(1)Toolbox中,双击/SARscape/Interferometry/Phase Processing/2 - Adaptive Filter and Coherence Generation打开Adaptive Filter and Coherence Generation面板,在Adaptive Filter and Coherence Generation面板:

  • Input Files面板
    • Interferogram File项,选择上一步生成的干涉图_dint;
    • Input Master File和Input Slave File两项,会自动填入上一步生成的主从影像的强度图_master_pwr和_slave_pwr。
  • Parameters面板:主要参数(Principal Parameters)
    • Coherence Generation:True
    • Adaptive Filter:True
    • Filtering Method:Goldstein
    • Coherence from Fint:True

提供了三种滤波方法:

  • Adaptive:这种方法适用于高分辨率的数据(如TerraSAR-X或COSMO-SkyMed)
  • Boxcar :使用局部干涉条纹的频率来优化滤波器,该方法尽可能的保留了微小的干涉条纹。
  • Goldstein:这种滤波方法的滤波器是可变的,提高了干涉条纹的清晰度、减少了由空间基线或时间基线引起的失相干的噪声。这种方法是最常用的方法。
  • Output Files面板,输出根名称按照默认。

(2)单击Exec执行。

图:干涉图滤波和相干性计算

处理完成后,弹出Completed对话框,点击确定。查看输出路径下的结果,这一步生成的结果有:

  • _fint:滤波后的干涉图
  • _cc:相干性系数图

在ENVI中显示了相干系数图_cc,点击查看像元值,相干性系数分布在0-1,值越大说明该区域的相干性越高,值越小,相干性越低。

图:滤波后的干涉图_fint

本数据对覆盖区域90%以上区域为山区且植被茂密,从CC图上看到这个数据对相干性非常好,可能得益于L波段成像。

图:相干系数图_cc

5相位解缠

将以2π为周期循环变化的相位解算到连续变化的相位。

(1)Toolbox中,双击/SARscape/Interferometry/Phase Processing/3 - Phase Unwrapping,打开Phase Unwrapping面板;

(2)在Phase Unwrapping面板

  • Input Files
    • Coherence File项,选择上一步生成的相干性结果_cc;
    • Interferogram File项,自动选择_fint结果。
  • Parameters项,主要参数Principal Parameters
    • 解缠方法(Unpwrapping Method Type):Region Growing
    • 解缠分解等级(Unwrapping Decomposition Level):1
    • 解缠相干性阈值(Unwrapping Cohrence Threshold):15

说明:

1)解缠方法(Unwrapping Method Type)提供了三种:

  • 区域增长法(Region Growing):选这种方法,则不要设置过高的相干性阈值(15-0.2是比较好的)以便留下足够的自由增长空间,相位突变部分在解缠后的图像上以解缠孤岛存在,这种方法降低了由相位突变引起的误差,
  • 最小费用流(Minimum Cost Flow):默认的解缠方法,当有大面积的低相干或是其他限制增长的因素而使解缠困难时,最小费用流算法可以取得比区域增长法更好的结果。这种方法采用正方形的格网,考虑了图像上所有的像元,对相干性小于阈值的像元做了掩膜处理。
  • Delaunay MCF:和最小费用流法的不同在于,这种方法不是考虑了图像上所有的像元,而是仅考虑了相干性大于阈值的部分,而且不是用正方形的格网而是用了德罗尼三角形格网。只有对相干性高的部分进行解缠,不受低相干像元的影响。对于有大量相干性低的地物存在的时候,如影像上存在大量水体、浓密植被等,推荐使用该方法,最小化相位突变的影响。

2)分解等级(Unwrapping Decomposition Level):该分解是为了用迭代的方法对数据做多视和疏采样:干涉图以一个较低的分辨率被解缠然后被重采样成原来的分辨率。使用了分解可减少解缠错误,提高处理效率。用户可以指定迭代的次数,每个迭代相当于3次采样过疏。这里可以填的数有:-1、0、1、2、3。其中,-1和0代表不执行分解,用原始的像素采样,当形变很大或是地形很陡峭的情况下(多相位/高密度分布),分解可引起交迭效应,建议设置为-1或0;1代表最小的分解等级。3:最大的分解等级。建议这个次数不要超过3。通常情况设置原始的像素采样(如-1)或者最小的分解等级(1)。

  • Output Files
    • Unwrapped Phase File项,输出解缠结果,路径和文件名自动添加。

(3)单击Exec,开始执行。

图:相位解缠

处理完成后,弹出Completed对话框,点击确定。自动加载解缠后的结果。

图:解缠结果图

6 轨道精炼和重去平

当输入精确的轨道信息,为了矫正相位偏移,一般需要进行轨道参数的修正,矫正的结果不会生成新的数据文件,而是将解缠图像头文件中的信息做了修正。

(1)在Toolbox中,选择/SARscape/Interferometry/Phase Processing/4 - Refinement and Re-flattening,打开轨道精炼和重去平面板,在Refinement and Re-flattening面板

  • 数据输入面板(Input Files)
    • Coherence File项,选择_cc图像输入;
    • Input Master File项,自动根据上一项添加的路径选择主影像的强度图输入_master_pwr;
    • Input Slave File项,自动添加从影像的强度图_slave_pwr;
    • Unwrapped Phase File项,自动添加解缠后的相位图_upha;
    • Synthetic Phase File项,自动添加合成相位图_sint;
    • Slant Range DEM File项,自动添加斜距DEM数据_srdem;
    • Refinement GCP File,控制点文件,点击,打开生成GCP流程化工具。

1)数据输入:在生成GCP的面板,输入数据_upha、dem文件、参考数据_fint会自动添加,在文件输入面板上点击Next 按钮。

图:选择控制点时的输入文件

2)在控制点选择面板,鼠标在图像中选择控制点,控制点一般选择相位相干性好的,尽量选择平地区域,避免相位突变的区域。控制点选好之后切换到Cartographic System项,默认为参考DEM的坐标系GEO-GLOBAL WGS84;切换到Export项,默认控制点文件输出的路径和文件名_upha_refinement_and_reflat_In_gcp.xml。点击Finish按钮。控制点文件生成并自动添加到面板中。

说明:控制点是用来估算轨道精炼时候的修正参数的,去除系统的几何误差(轨道误差),控制点的选择依据以下原则:

1)在差分干涉图(_dint或_fint)和解缠图(_upha)选择控制点,避免有地形相位没有去除的区域和变化的区域(干涉条纹密集区域);

3)选择相干性高的区域;

4)远离形变区;

5)避免解缠错误的区域,不能位于解缠错误的相位跃变上(phase jump),如相位孤岛等;

通俗的讲:GCP要远离形变区,认为形变为0,陡峭的地形区域和有残余地形相位区域,当地形起伏大的山区,最好是选择山谷底部的平地区域。

选择好GCP之后,点击Finish,自动设置输出GCP点的文件名为:*_upha_refinement_and_reflat_gcp.xml。

  • 可选文件面板(Optional Files):干涉图作为参考文件,自动输入_fint。
  • DEM/Cartographic System面板:设置参考地形,输入参考DEM文件。
  • 参数设置面板(Parameters):
    • 轨道精炼方法(Refinement Method):Polynomial Refinement
    • Refinement Res Phase Poly Degree:3(如果控制点数量不够3次多项式,会自动降低多项式次数)
  • 数据输出面板(Output Files):默认数据输出的路径及文件名,自动添加_reflat后缀。

 轨道精炼的三种方法:

  • 自动优化(Automatic Refinement):这种方法是首先根据输入的控制点估算轨道形态,如果“Achievale RMS”大于阈值、或者空间基线的绝对值小于Minimum Baseline、或者“Final RMS”大于阈值、RMS Ratio大于阈值、GCP点的个数小于7、会自动切换到轨道优化方法。
  • 轨道优化(Orbital Refinement):这种方法是根据控制点估算轨道参数,这种方法唯一的前提条件就是控制点个数大于7。
  • 多项式优化(Polynomial Refinement):这种方法是从解缠后的相位中估算相位斜坡,不考虑轨道形态,这种方法控制点个数必须与多项式次数所需要的控制点个数对应(如3次多项式需要的控制点个数至少是10),否则多项式次数会自动降低。

(2)单击Exec执行

处理结束后,在Completed对话框单击end。自动打开轨道精炼计算出的偏移量面板,内容如下:

此外,还生成重去平后的一系列结果,结果文件带有_reflat 标识。浏览轨道精炼后的干涉图、解缠图等结果。

控制点的查看修改方法如下:

1) 在*_reflat_refinement.shp图层右键选择Properties,在打开的矢量属性面板上选择Attribute为AbsRHgtDiff,Color Table为RAINBOW,点击OK,控制点会进行彩色渲染显示,颜色越红的点,说明误差越大。

2) 在SHP图层右键选择View/Edit Attributes,在AbsRHgtDiff字段右键Sort By Selected Column Reverse进行排序。

3) 记录下点号shp_ID(不关闭对话框),然后在轨道精炼工具中,打开加载进来控制点,把误差大的几个点删除,再生成一次点文件,同样的方法再进行一次轨道精炼。

7 相位形变转换——生成形变结果

这一步是将经过绝对校准和解缠的实际相位,转换为形变结果并进行地理编码。生成地理编码的形变结果文件、相干图像,图像精度图像和分辨率图像。

(1)在Toolbox中,选择/SARscape/Interferometry/Phase Processing/5B - Phase to Displacement Conversion and Geocoding,打开Phase to Displacement Conversion and Geocoding面板,在Phase to Displacement Conversion and Geocoding面板

  • 数据输入(Input Files)面板:
    • Coherence File下,选择相干性图_cc输入;
    • Unwrapped Phase File:自动从上一个路径下找到相应的文件输入_reflat_upha;
  • DEM/Cartographic System面板:选择DEM文件的坐标系作为输出DEM产品的坐标系。
  • Parameters面板:
    • 产品相干性阈值(Product Coherence Threshold):25,相干性大于这个值的部分生成形变产品。
    • 垂直形变(Vertical Displacement):对形变结果进行三角解算生成垂直方向的形变。
    • 坡度形变(Slope Displacement):对形变结果进行三角解算生成最大坡度方向的形变。
    • 自定义方向(Displacement Custom Direction):对形变结果进行三角解算生成自定义方向上的形变,由方位向和倾向来确定自定义方向。
    • 方位角(Azimuth Angle):
    • 倾角(Inclination Angle):
    • X Dimension(m):12,X方向上的制图分辨率
    • Y Dimension(m):12,Y方向上的制图分辨率
    • 内插窗口大小(Interpolation Window Size):7
    • 对无效值内插处理(Relax Interpolation):True。对相位缺少区域进行插值处理
    • 对图像外的无效值做掩膜(Dummy Removal):True
  • Output Fiiles面板:默认输出的路径和文件名。

(2)点击Exec执行。

图:相位转形变