陆探一号（LT-1）InSAR处理方法

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://pan.baidu.com/s/1-K6kmi20To-JAaJYJw0Zew?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格式的图像轮廓线。

子区裁剪（可选）

这一步为可选项，使用的工具为：SARscape/General Tools/Sample Selections/Sample Selection SAR Geometry Data。

2基线估算

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

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

图：基线估算结果

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

  • 空间基线为114米（远小于临界基线23165米）；
  • 时间基线为4天；
  • 高程的变化是464.53米（2 PI Ambiguity height (InSAR) (m)）；
  • 相位变化一个2π周期，代表着地表发生了0.119米的形变；
  • 立体像对1个像素移动模糊高程为6503.422米；
  • 振幅跟踪1个像素移动模糊形变为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）：0.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执行。

图：相位转形变

posted @ 2023-06-12 20:44  ENVI-IDL技术殿堂  阅读(5368)  评论(1)  编辑  收藏  举报