涪城一号InSAR处理方法详解

涪城一号卫星是由天仪研究院研发并运营的第3颗SAR卫星，也是在海丝一号、巢湖一号基础上研制的首颗面向InSAR业务的商业SAR卫星，于2023年6月7日在酒泉卫星发射中心成功发射。该卫星采用C波段，具备全天时、全天候的干涉成像能力。支持多星组网实现高频次观测；使用GPS+北斗高精度四频GNSS接收机实现精密定轨（精度优于10厘米），并具备优秀的重轨干涉能力。涪城一号支持聚束、条带及扫描等多种成像模式，适应不同的应用场景。

SARscape6.1已经开始支持涪城一号的导入，数据导入后就能进行后续的所有处理，包括一般处理（多视、滤波、地理编码等）、InSAR和DInSAR处理、多孔径干涉测量、振幅偏移量量测、时序SAR处理（PS和SBAS）等。

本文是在SARscape6.1下完成，例子数据为条带成像模式。处理包括数据导入、基线估算、干涉图生成、滤波及相干性计算、相位解缠、轨道精炼与重去平等。

软件下载和试用：https://envi.geoscene.cn/sar_license

涪城一号

涪城一号搭载了C波段SAR载荷，具备5种成像模式，最高分辨率1米。如下表为涪城一号主要卫星指标。

卫星名称

涪城一号

发射时间

2023.6.7

运行轨道

太阳同步轨道

类型

C波段SAR

极化方式

VV

轨道倾角

97.43°

轨道高度

508km

干涉重访周期

11天

左右视能力

均具备

成像重访周期

2-3天

降交点地方时

10:30(±25min)

视角

15~40°

成像模式

聚束、条带、SCANSAR、TOPSAR

定轨精度

实时定轨<10m

事后精密定轨<10cm

分辨率

聚束（SP）：1m

滑聚（SSP）：2m

条带（SM）：3m

TOPS(TP):10m

NS/ES(窄扫/宽扫)：12m/20m　

幅宽

聚束（SP）：7km

滑聚（SSP）：15km

条带（SM）：25km

TOPS(TP):100km

NS/ES(窄扫/宽扫)：100km/175km

在应用领域方面，涪城一号广泛应用于地质灾害监测、城市安全、应急救灾、基础设施监测等。它能够提供高质量、常态化的国产商业InSAR数据，支持对道路桥梁、水库大坝、隧道涵洞等基础设施形变，以及滑坡点（群）监测与隐患点识别等多领域的综合监测应用。此外，涪城一号还成功应用于复杂山区滑坡形变监测，其数据质量与干涉效果已获得业内专家和用户的一致好评。

数据导入

涪城一号产品命标准如下：

本文档使用的一个数据为：

BC3-SP-SLC-1SVV-20240611T054940-007285-000109-001C75.zip

它是单视复数据（SLC），1m聚束成像模式。感谢天仪空间科技研究院提供例子数据。

在处理之前，我们选择一套默认参数，打开Toolbox/SARscape/Preferences/Preferences specific，在打开的界面中，选择Load Preferences->VHR（better than 6m），Cartographic Grid Size（m）：1。

注：如果覆盖区域植被较多等情况，相干性可能较差，可以考虑降分辨率处理。

（1）在Toolbox中，选择/SARscape/Import Data/SAR Spaceborne/Single Sensor/Fucheng-1。

注：打开的时候会弹出帮助文档，以及一个对该数据初步支持的提示框，点击OK即可。

（2）在打开的面板中，

  • 数据输入面板（Input Files）
  • 输入文件（Input File List）：输入过滤的.safe文件
  • 参数设置面板（Parameters）：主要参数（Principal Parameters）
  • 极化方式（Polarization）：ALL，输出所有的极化数据，可以选择只输出同极化或者交叉极化的数据；
  • 对数据重命名（Rename the File Using Parameters）：True。软件会自动在输入文件名的基础上增加几个标识字母，如增加“_HH_slc”。
  • 对无效值填充（Fill dummy during import）：True
  • 数据输出面板（Output Files）

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

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

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

（3）单击Exec按钮开始执行。

图1 数据导入

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

子区裁剪（可选）

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

可使用有地理坐标或无地理坐标的Shapefile、KMZ、kml文件，甚至是个角点坐标对SLC文件进行裁剪。本实例进行了裁剪处理。

DEM获取

DEM一般作用在InSAR处理中的干涉图去平、图像配准、地理编码等步骤中，一般选择SRTM、ALOS World 3D等。SARscape支持自动在线下载、外部DEM导入等方式。详细可参考：

https://www.cnblogs.com/enviidl/p/16469506.html

本文直接使用InSAR流程化处理工具中自动下载SLC范围内的DEM数据。

基线估算

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

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

图2 基线估算结果

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

  • 空间基线为104米（远小于临界基线11842米）；
  • 时间基线为11天；
  • 高程的变化是013710米（2 PI Ambiguity height (InSAR) (m)）；
  • 相位变化一个2π周期，代表着地表发生了027759米的形变；
  • 立体像对1个像素移动模糊高程为951635米；
  • 振幅跟踪1个像素移动模糊形变为454231米；

该数据对适合标准InSAR处理、Persistent Scatterers Stacking, Distributed Scatterers Stacking、或者Amplitude Tracking处理。

InSAR流程化处理工具

启动流程化的InSAR处理工具：/SARscape/Interferometry/InSAR DEM Workflow，打开InSAR DEM面板。流程化的处理工具的左侧是工作流的各个步骤，右侧是相应步骤的输入输出及参数设置。右下方是步骤控制按钮。

（1）数据输入（Input）

  • Input File面板，Input Reference file项，输入前时相SLC作为参考影像，Input Secondary file项，输入前时相SLC作为第二影像。
  • DEM/Cartographic System面板，输入参考高程，提供了两种类型，分别是输入已有的参考DEM文件、自动下载相应区域的DEM数据。这里自动下载相应区域的DEM数据：DEM Download。
  • Parameters面板，有制图分辨率大小的设置，这里为Grid Size：1米。

点击Next按钮，到Inport Generic SAR Data步骤，这一步是进行数据导入的，数据已经是导入后的SLC，故点击Next，到下一步。

（2）干涉图生成（Interferogram Generation）

  • 主要参数（Principal Parameters）
    • 距离向视数（Range Looks）：1。自动添加。
    • 方位向视数（Azimuth Looks）：3。自动添加。
    • 制图分辨率（Grid Size for Suggested Looks）：1。根据之前的设置自动添加。
    • 配准时使用DEM（Coregistration With DEM）：True。激活该项，配准时光谱偏移计算就会考虑到局部的地形信息，若轨道参数不是非常精确的话建议不激活。以下几种情况激活：1）长条带的数据；2）高纬度区域的数据；3）带有非0多普勒注释的数据（尤其是波段较长的数据如ALOS PALSAR）。激活该选项所需的数据处理时间会比较长。
  • 全局参数（Global）
  • 生成TIFF数据（Generate Quick Look）：Ture，生成TIFF格式的中间结果，经过彩色渲染的干涉快视图，可作为插图使用，其他步骤类似。

注：1.建议在开始的时候设置好制图分辨率，在此处使用自动添加的视数，不再做手动修改。2.绝大多数情况使用默认能得到较好的结果。

图3 干涉处理步骤

点击Next按钮，进行干涉图生成处理。处理完成后，自动加载了去平后的干涉图、以及主从影像的强度图。生成的结果存放在ENVI默认的输出路径下，默认为：C:\Users\<用户名>\AppData\Local\Temp，自动生成了一个名为SARsTmpDirXXXXXXXXX的文件夹，这一步生成的结果有：

  • _dint：去平后的干涉图
  • _int：干涉图
  • _ reference _pwr：参考影像强度图
  • _ secondary _pwr：第二影像强度图
  • _par：文本文件保存的配准时的偏移量数据，
  • _sint：合成的相位
  • _srdem：斜距几何下的参考DEM
  • shp：估算轨道偏移时用到的点矢量数据，包括在方位向和距离向 的点的位置坐标、测量到的偏移量、计算出的线性偏移量。
  • shp：从强度数据的相差上估算交叉相关偏移量的点矢量数据。包含每个点上交叉相关的值（CC），范围是0-1。
  • shp：从相位数据的相差上估算相干性的点矢量数据，包含信噪比（SNR）和相干性，相干性值的范围是0-1。

注：C:\Users\<用户名>\AppData文件夹默认是隐藏，需要显示隐藏文件夹；可以事先在ENVI中File->preference将临时路径设置到其他盘的路径下。

如下图为去平后干涉图。

图4 去平地相位后的干涉图_dint（快视图）

（3）滤波和相干性计算（Adaptive Filter and Coherenace Generation）

  • 选择主要参数（Principal Parameters）：
  • 滤波方法（Filtering Method）：Goldstein

提供了三种滤波方法：

  • Adaptive

这种方法适用于高分辨率的数据（如TerraSAR-X或COSMO-SkyMed）

  • Boxcar

使用局部干涉条纹的频率来优化滤波器，该方法尽可能的保留了微小的干涉条纹。

  • Goldstein

最常用的方法。这种滤波方法的滤波器是可变的，提高了干涉条纹的清晰度、减少了由空间基线或时间基线引起的失相干的噪声。

如果滤波效果不好，可以修改滤波参数，如下为Goldstein滤波参数的修改：

选择滤波参数（Filtering）：调整以下三个参数：

  • Goldstein Win Size：默认为32，设置大一些，如128。该值越大，滤波器对小细节(如局部干涉条纹)的灵敏度就越低，该值最好使用2的倍数，该值范围为：32~512。
  • Goldstein Min Alpha：默认3，设置大一些，如1.5。值越大滤波器平滑越强，该值范围为：0.3~3。
  • Goldstein Max Alpha：默认5，设置大一些，如3，值越大滤波器平滑越强，该值范围为：0.5~4。

本例中不修改这三个参数。

  • 全局参数（Global）

生成TIFF数据（Generate Quick Look）：Ture，生成TIFF格式的中间结果，经过彩色渲染的干涉快视图，可作为插图使用。

点击Next按钮，进行干涉图滤波和相干性生成处理，处理完成之后，自动加载滤波后的干涉图_fint和相干性系数图_cc。这一步处理之后生成的结果有：

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

如下图为滤波之后的干涉图，如下图为去平后干涉图，图上可以看到有个明显的干涉条纹，有可能是11天内有形变现象发生，也有可能是使用的ALOS World 3D生产后地表发生了比如挖掘活动造成的。

图5 干涉图滤波结果fint（快视图）

图6 相干系数图_cc

（4）相位解缠（Phase Unwrapping）

相位的变化是以2π为周期的，所以只要相位变化超过了2π，相位就会重新开始和循环。相位解缠是对去平和滤波后的相位进行解缠处理，使之与线性变化的地形信息对应，解决2π模糊的问题。

  • 主要参数（Principal Parameters）
    • 解缠方法（Unwrapping Method Type）：Minimum Cost Flow
    • 解缠分解等级（Unwrapping Decomposition Level）：1
    • 解缠最小相干性阈值（Unwrapping Coherence Threshold）：25。相干系数小于该阈值的像元不进行解缠。

注：可以打开上一步获取的CC相干图浏览直方图，看低于0.25的像素占比判断这个0.25是否合理。

图7 相位解缠参数设置面板

说明：

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。

点击Next按钮，进行干涉图滤波和相干性生成处理，处理完成之后，自动加载滤波后的相位解缠结果图_upha。这一步处理之后生成的结果有：

  • _upha：相位解缠结果

图8 相位解缠结果_upha

（5）控制点选择（GCP Selection）

输入用于轨道精炼的控制点文件，可以用已有的GCP点文件也可在此选择控制点。

  • 在Refinement GCP File（Mandatory）项中，点击按钮，自动打开流程化的控制点选择工具，并输入了相应的参考文件。

图9 控制点生成面板

在控制点生成面板上，点击Next，打开控制点选择工具，鼠标变为选点状态，在图像上单击鼠标左键，选择控制点。

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

  • 选择平地区域
  • 选择相干性高的区域；
  • 控制点尽量分布于整个范围内；
  • 将去平地之后的干涉图上（_dint或_fint）和解缠图，切换显示，选择控制点，避免干涉条纹密集区域（说明存在地形变化）以及解缠跃变区域（解缠不平滑）；
  • 选择相干性高的区域；
  • 控制点尽量分布于整个范围内；
  • 避免解缠错误的区域，不能位于解缠错误的相位跃变上（phase jump），如相位孤岛等；

通俗来说，GCP要远离陡峭的地形区域和有残余地形相位区域，当地形起伏大的山区，最好是选择山谷底部的平地区域。

  • 单击Cartographic System按钮，查看控制点的参考坐标系统，该坐标系是从参考DEM上自动读取的。
  • 单击Export按钮，查看控制点的存放路径和文件名。生成的控制点文件为xml。

注：如果默认输出的路径有中文字符的文件夹或者文件名，需要改到只有英文字符的文件路径中。

在控制点生成面板上点击Finish，生成了控制点文件，并自动添加到InSAR流程化处理面板的Refinement GCP File（Mandatory）项。在面板上点击Next按钮，进入下一步。

图10 生成的控制点文件

（6）轨道精炼和重去平（Refinement and Reflattening）

进行轨道精炼和相位偏移的计算，消除可能的斜坡相位，对卫星轨道和相位偏移进行纠正。这一步对解缠后的相位是否能正确转化为高程或形变值很关键。

  • 主要参数（Principal Paremeters）
    • 轨道精炼方法（Refinement Method）：Polynomial Refinement
    • 轨道精炼的多项式次数（Refinement Res Phase Poly Degree）：3。在重去平的过程中用到的估算相位斜坡的多项式次数，若输入的控制点个数较少，次数会自动降低。默认的3表示在距离向和方位向加上一个恒定的相位偏移的相位斜坡会被修正，如果仅需要相位偏移校正，这个次数可以设置为1。
    • 配准时是否考虑到地形因素（Coregistration With DEM）：True。

说明：

轨道精炼的三种方法（在Refinement参数面板中可设置）：

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

图11 轨道精炼和重去平面板

点击Next按钮，进行轨道精炼和重去平处理，处理完成之后，将优化的结果显示在Refinement Results的面板，内容如下：

ESTIMATE A RESIDUAL RAMP

Valid points found = 18

Extra constrains = 2

Polynomial Degree choose = 3

Polynomial Type : = k0 + k1*rg + k2*az

 -1.1408671713

 0.0000328031

 0.0001629409

Polynomial Coefficients (radians) :

              k0 = -1.1408671713

              k1 = 0.0000328031

              k2 = 0.0001629409

Root Mean Square error (m) = 5.3055124701

Mean difference after Remove Residual refinement (rad)  = 0.0054496954

Standard Deviation after Remove Residual refinement (rad)  = 0.5666274461

 

这一步得到的结果有：

  • _reflat_fint：重去平后的干涉图
  • _reflat_sint：重去平后的合成相位图
  • _reflat_upha：重去平后的解缠结果
  • _reflat_srdem：重去平后的斜距DEM
  • txt：轨道精炼所用的轨道修正参数
  • shp：轨道精炼用到的有效的控制点文件，斜距坐标系。
  • shp：轨道精炼用到的有效控制点文件，地理坐标系。

浏览轨道精炼后的干涉图、解缠图等结果。也可以分析控制点的质量，方法如下：

  1. 在ENVI中打开SARsTmpDir_DDMMYYY_HHMMSS路径下的shp。临时路径默认为C:\Users\<用户名>\AppData\Local\Temp\SARsTmpDir_ DDMMYYYY _HHMMS）。
  2. 在SHP图层右键选择Properties，在打开的矢量属性面板上选择Attribute为AbsRHgtDiff，Color Table为RAINBOW，点击OK，控制点会进行彩色渲染显示，颜色越红的点，说明误差越大。
  3. 在SHP图层右键选择View/Edit Attributes，在AbsRHgtDiff字段右键Sort By Selected Column Reverse进行排序。
  4. 记录下点号shp_ID(不关闭对话框)，然后在InSAR工作流中，点击Back，回到选择控制点的一步，加载进来控制点，把误差大的几个点删除，再生成一次点文件，同样的方法再进行一次轨道精炼。

（7）相位转高程以及地理编码（Phase to Height Conversion and Geocoding）

将经过绝对校准和解缠的相位，结合合成相位，转换为高程数据以及地理编码到制图坐标系统，

  • 主要参数（Principal Parameters）
    • 相干性阈值（Product Coherence Threshold）：25。相干性大于该值的进行DEM产品的输出。
    • 空间滤波大小（Spatial Wavelet Size（m））：该值结合雷达数据的分辨率以及参考DEM的分辨率，决定了在相位转换高程时所保留的相位。一般设置为参考DEM分辨率的4倍，此处设置为120m。
    • 生成栅格文件（Genetate Raster Format）：True
    • 生成矢量文件（Genetate Shape Format）：False。将栅格DEM转为矢量数据，数据量很大，耗费时间长，按照默认参数False，不进行栅矢转换。
    • 生成LAS格式（Generate Las Format）：False。不生成点云格式的DEM。
    • 生成的DEM类型（Output Type）：默认是Ellipsoidal椭球高
    • 大地水准面（Geoid Type）：EGM96
    • X方向上的水平分辨率（X Dimension（m））：1。制图分辨率，X方向。
    • Y方向上的水平分辨率（Y Dimension（m））：1。制图分辨率，y方向。
  • 地理编码参数（Geocoding）
    • 对无效值内插处理（Relax Interpolation）：True。用内插的方法处理无效值。
    • 去除图像外的无用值（Dummy Removal）：True。对图像外的区域做掩膜处理

图12 相位转高程

点击Next按钮，进行相位转高程和地理编码处理。地理编码的坐标系是以参考DEM的坐标系为准。这一步得到的结果有：

  • _demwf_dem：DEM产品
  • _demwf_cc_geo：地理编码的相干性系数图
  • _demwf_precision：数据质量的估算结果图，代表DEM的精度。
  • _demwf_resolution：基于局部入射角得到的空间分辨率

（8）结果输出（Output）

   结果默认输出在ENVI的默认输出路径下，文件名中包含 _demwf。若想保留中间结果便于查看，将Delete Temporary Files不勾选。

图13 结果输出

点击Finish，输出结果。结束InSAR DEM处理的工作流。得到的DEM数据自动进行密度分割配色展示。

图14 提取结果

注：Back和Next按钮，可切换至中间某一步查看参数或调整参数重新处理。

与30米的ALOS World 3D对比，很多小地形能表现出来，如下图所示。

图15 结果对比（右-80%缩放）

图16 结果对比（右-100%缩放）

posted @ 2025-04-17 17:39  ENVI-IDL技术殿堂  阅读(1516)  评论(0)    收藏  举报