四维高景雷达星InSAR处理方法

本文主要是介绍四维高景雷达星InSAR处理方法，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

X波段的四维高景卫星系统首批雷达星（简称SV2）于2022年7月16日成功发射，并于2023年3月22日，中国四维测绘技术有限公司发布了国内首个干涉SAR商业卫星星座的数据产品。SV2采用双星编队InSAR模式和分布组网D-InSAR模式。包括2m条带模式（STRIP）、1m滑动聚束模式（SPOT）、0.5m聚束模式（FineSPOT）成像模式。

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

本文是在SARscape5.7下完成，包括数据导入、基线估算、干涉图生成、滤波及相干性计算、相位解缠、轨道精炼与重去平等。

数据导入

SV2搭载了X波段SAR载荷，具备3种成像模式，最高分辨率0.5米。如下为图像产品命令说明。

产品文件名称：

____卫星标识_单/双站_接收地面站代号_成像模式_接收圈号_景中心经度_景中心纬度_采集日期_产品类型_极化方式_产品级别_产品号.格式扩展名

____示例：

SV2-04_MONO_HYKSC1_STRIP_002776_E116.4_N33.6_20230113_SLC_HH_L1A_0000076170 .tiff

____卫星标识= SV2-03，SV2-04

____单/双站= MONO代表单站，BIST代表双站

____接收地面站代号=

北极站：KSBJ71

铜川站：XYTCC1

喀什站：HYKSC1

广州站：SWGZ71

西安站：SHXA71

鸡西站：HYJXC1

精河1.8米S/X站：XJJH11

____成像模式= STRIP代表2m条带模式，SPOT代表1m滑动聚束模式，FineSPOT代表0.5m聚束模式

____采集日期=2023（年）01（月）13（日）

____产品类型=SLC，GEC

____极化方式=HH

____产品级别=L1A，L2A

本文档使用的是SV2 L1A级产品，它是单视复数据（SLC），2m条带模式成像模式。感谢中国四维测绘技术有限公司提供例子数据。

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

注：这里为了提高数据对相干性，设置为6米。

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

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

（2）在打开的面板中，

数据输入面板（Input Files）
输入文件（Input File List）：输入过滤的.meta.xml文件

注：SARscape是支持批量导入文件，如果需要导入文件角度，可以在选择文件的对话框右上角输入.meta.xml进行检索，一次性选择多个搜索的*.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。

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

DEM获取

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

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

本文直接使用工具/SARscape/Import Data/DEM Extraction/ALOS World 3D 30m自动下载SLC范围内的DEM数据。

基线估算

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

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

图：基线估算结果

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

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

该数据对适合InSAR\DinSAR或者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文件，Reference Type：Input DEM，在Input DEM项，输入之前准备好的参考DEM文件AW-3d-30m_dem。
Parameters面板，有制图分辨率大小的设置，这里为Grid Size：6米。

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

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

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

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

图1 干涉处理步骤

点击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将临时路径设置到其他盘的路径下。

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

图2 去平地相位后的干涉图_dint

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

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

提供了三种滤波方法：

Adaptive

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

Boxcar

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

Goldstein

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

全局参数（Global）

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

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

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

图3 干涉图滤波结果fint

图4 滤波后的干涉图_fint和相干系数图_cc

（4）相位解缠（Phase Unwrapping）

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

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

注：可以打开上一步获取的CC相干图浏览直方图，看低于0.2的像素占比。

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

说明：

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：相位解缠结果

图6 相位解缠结果_upha

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

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

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

图7 控制点生成面板

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

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

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

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

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

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

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

图8 生成的控制点文件

（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），否则多项式次数会自动降低。默认使用该方法。

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

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

ESTIMATE A RESIDUAL RAMP

Valid points found = 36

Extra constrains = 2

Polynomial Degree choose = 3

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

2.6152866286

-0.0007015282

-0.0004601463

Polynomial Coefficients (radians) :

k0 = 2.6152866286

k1 = -0.0007015282

k2 = -0.0004601463

Root Mean Square error (m) = 10.2014459812

Mean difference after Remove Residual refinement (rad) = 0.0468077120

Standard Deviation after Remove Residual refinement (rad) = 0.9218314442

这一步得到的结果有：

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

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

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

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

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

主要参数（Principal Parameters）
- 相干性阈值（Product Coherence Threshold）：2。相干性大于该值的进行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））：6。制图分辨率，X方向。
- Y方向上的水平分辨率（Y Dimension（m））：6。制图分辨率，y方向。
地理编码参数（Geocoding）
- 对无效值内插处理（Relax Interpolation）：True。用内插的方法处理无效值。
- 去除图像外的无用值（Dummy Removal）：True。对图像外的区域做掩膜处理