【时序InSAR开源软件】用isce2(stackSentinel.py)+mintpy进行哨兵一号SBAS小基线集InSAR时序处理全流程

---

前言

本博客仅作为个人学习笔记用途，欢迎和各位前辈交流学习

本文记录博主使用isce2+mintpy进行哨兵一号时序InSAR的SBAS形变提取全流程。

isce2 是JPL开发的用于时序InSAR前期干涉操作数据处理的开源软件，支持很多现有的SAR传感器，如ALOS, ALOS2, COSMO_SKYMED, Sentinel1, TERRASARX等等，除单景操作之外，还支持时序多景雷达影像stack处理。

mintpy 是开源时序InSAR分析软件，支持输入ISCE, ARIA, FRInGE, HyP3, GMTSAR, SNAP, GAMMA生成的干涉序列获取时序InSAR形变分析结果，值得一提的是，开发者提供了非常详实的教程、注释和细节功能，很有利于InSAR初学者。

选用数据： 使用20211204-20231124 共59景Sentinel-1

处理思路： 使用isce2软件中的stackSentinel.py获取干涉序列–>利用mintpy处理干涉序列获得形变结果

1.环境配置

使用的环境为WSL/Ubuntu/miniconda

mintpy和isce2可以直接使用conda进行安装（推荐）：参考这篇文章–>哨兵一号开源时序InSAR环境配置：WSL Ubuntu 环境下使用conda安装isce2+mintpy

isce2也可以WLS2+CMAKE+ISCE2完整安装，安装较为繁琐但是内容全面，参考链接–>ISCE2023学习 (0.0) WLS2+CMAKE+ISCE2 安装环境配置

经过上述环境配置，应该实现可以运行stackSentinel.py 和 smallbaselinApp.py两个软件，才能进行后续的工作

2.数据准备

stackSentinel.py需要输入以下文件：

1. SLC：原始单视复数雷达影像SLC文件序列
2. DEM：研究区域的DEM文件
3. OrbitFile：SLC文件对应的哨兵一号轨道文件
4. AuxFile：校准辅助（AUX_CAL）文件用于天线模式校正，（主要用于 2015 年 3 月之前获取的图像）如果日期不在该范围，这一文件夹可以忽略

这里给出我的文件夹结构：

其中Stacks为我的工作文件夹WORK_DIR，也即做好准备工作后，stackSentinel.py输出结果被保存在该文件夹。
手动创建文件夹或mkdir操作后，可以将以上地址添加到环境变量，便于后续操作：

注意将/YOUR_DATA_DIR替换为你的目录

export ORBIT_PATH='/YOUR_DATA_DIR/orbits'
export AUX_PATH='/YOUR_DATA_DIR/AUXFILE'
export DEM_PATH='/YOUR_DATA_DIR/DEM'
export SLC_PATH='/YOUR_DATA_DIR/SLC'
export WORK_DIR='/YOUR_DATA_DIR/Stacks'

如果你的操作平台是ipynb文件（也即jupyter notebook环境），你也可以在python环境中运行：

import os

# 设置环境变量 DEM_PATH
os.environ['DEM_PATH'] = '/YOUR_DATA_DIR/DEM'

# 设置环境变量 ORBIT_PATH
os.environ['ORBIT_PATH'] = '/YOUR_DATA_DIR/orbits'

# 设置环境变量 AUX_PATH
os.environ['AUX_PATH'] = '/YOUR_DATA_DIR/AUXFILE'

# 设置环境变量 SLC_PATH
os.environ['SLC_PATH'] = '/YOUR_DATA_DIR/SLC'

#设置工作地址 WORK_DIR
os.environ['WORK_DIR'] = '/YOUR_DATA_DIR/Stacks'

2.1 SLC数据获取

哨兵一号的SLC数据可以直接在ASF下属的EARTHDATA数据官网进行下载，按照网站上的说明进行数据筛选，类型勾选L1 Single Look Complex(SLC)即可。

也可以使用SSARA工具下载，主要是用于只能使用命令窗口以及python自动化下载的场景，具体可以参考ISCE2023学习（0.7）数据获取
获得的数据压缩包并不需要解压，SLC文件夹范例如下图所示：

2.2 DEM数据获取

DEM数据获取分为两种方法:
1.可以利用ISCE自带的的dem.py功能进行下载
这一方法较为快速，能够直接下载大范围的DEM数据并执行拼接，可控性较强，缺点是参数较为复杂。

使用dem.py要求$HOME目录下存在.netrc文件，这一文件储存了下载的源网站、用户和密码。用户账户参考网站。如果没有这一文件，可以运行：

注意不要直接复制运行
将myUsername替换为你的ASF DataSearch用户名，myPassword替换为密码后运行

echo "machine urs.earthdata.nasa.gov login myUsername password myPassword" > $HOME/.netrc

以下是我的.netrc文件内容，aabbccdd为ASF用户名， 123456789为对应的密码

准备工作完成后可以使用dem.py进行下载：

cd $DEM_PATH
dem.py -a stitch -b 22 23 113 114 -r -s 1 -c
cd ..

这里解释以下上边的参数：
-a, --action 执行操作，可选选项为stitch（执行dem拼接）和download（仅下载），默认情况下，若你给出的bbox较大，下载了多幅DEM，会自动执行拼接。
-b, --bbox 研究区域的bounding box, [SNWE]，下载的数据默认都是以度为单位的大小
-r, --report 报告运行情况
-s, --source DEM来源，可选参数为1或者3
-c, --correct 将数据从EGM96转到WGS84，注意我们的数据的基准全部都在WGS84下进行
获得结果如下图所示：

其中.wgs84文件为我们需要的DEM文件

如果dem.py运行存在问题，可以参考文章isce2学习（三）下载哨兵数据、POD精轨数据、DEM数据、AUX_CAL辅助文件下载中更为详尽的解释。

---

2.利用OpenTopography网站进行下载
OpenTopography提供了DEM下载网站，需要进行注册，登录后进入页面，点击Select A REGION按钮，选择对应的研究区域，下方会给出包含该研究区域的对应DEM数据库，选择你想用的数据库点击。通常情况下，我们选择下图中的DEM，博主使用过SRTM GL1，也和上文中dem.py下载的数据一致，但是获取的数据为tif文件。

提交申请后会有等待界面，你也可以在提供的邮箱中收到处理完毕的链接🔗，在里边直接下载即可，获取的压缩包解压到$DEM_PATH下。

2.3 AUX文件获取

官方解释：

The following calibration auxliary (AUX_CAL) file is used for antenna pattern correction to compensate the range phase offset of SAFE products with IPF verison 002.36 (mainly for images acquired before March 2015). If all your SAFE products are from another IPF version, then no AUX files are needed.

也就是说如果你的数据不是在2015年3月前的，一般情况下是不需要AUX_CAL的数据的，当然为了稳妥也可以下载，直接访问链接🔗https://sar-mpc.eu,选择你需要的AUX_CAL文件即可。

2.4 Orbit哨兵一号轨道文件获取

这里使用了开源软件sentineleof进行轨道文件自动下载，很方便
首先安装：

pip install sentineleof

或者用conda安装：

conda install -c conda-forge sentineleof

可以运行eof --help查看是否安装成功
请注意，sentineleof在运行过程中也会读取你的.netrc文件，因此必须按照2.2 DEM数据获取中的教程编辑你的.netrc，如果已经有.netrc文件，则不需要改正。
准备完后，直接在命令行中运行：

eof -p $SLC_PATH --save-dir $ORBIT_PATH

其中-p/--search-path 后跟的是你的SLC文件所在的目录，--save-dir是你要把轨道文件保存在的目录，sentineleof会直接识别你的SLC文件并下载对应的精密轨道文件。

3. 使用stackSentinel.py获取干涉图序列

stackSentinel.py并不是一个一键完成的命令，它的功能是读取你给出的参数，并生成对应的runfile执行文件，以及对应config参数信息，在运行完stackSentinel.py获取runfile之后，我们才正式开始进行干涉操作，这里给出我的运行参数：

stackSentinel.py -s $SLC_PATH -w $WORK_DIR -a $AUX_PATH -d $DEM_PATH/demLat_N22_N23_Lon_E113_E114.dem.wgs84 -o $ORBIT_PATH -n 3 -b '22.28 22.35 113.86 113.97' --useGPU -C geometry -c 3 -z 1 -r 3 -f 0.6 --num_proc 4 

-s\-w\-a\-d\-o分别代表着我们刚刚获得的几种文件（SLC、工作目录、AUX校准、DEM、精密轨道），代码输出的结果会被保存在-w指定的$WORK_DIR工作目录中。

-n指定IW序列号，可选的参数为1、2、3，选定对应的IW编号后将仅提取你指定的IW条带，如果你的研究区域较小，这一参数能够帮你大幅节省空间。

-b 指定研究区域，注意stackSentinel并不会帮你按照bbox裁剪文件，而是根据给出的bbox来提取burst，并且至少会提取3个burst

--useGPU在运行过程中使用GPU进行运算加速

-C配准方案，有两个选项:geometry和NESD, NESD是在geometry几何配准的基础上，根据burst重叠区域进一步进行精化配准，选定NESD将大幅增加运算时间和储存空间，理论上讲相应的配准效果也会有所提升（然而在我的个人案例中并没有感受到多大的提升）

-c（注意大小写c是不一样的） 连接数量，主影像和相邻该间隔数的副影像形成干涉对，如我这里选定的数为3，则会生成 ab\ac\ad\bc\bd\be\cd\ce\cf······

-z\-r多视处理参数，分别为方位向\距离向

-f滤波强度，默认0.5

--num_proc 指定同时运行的任务数量，可以理解为多线程运行，后边的每个runfile文件其实就是一堆python文件运行命令，我这里设置为4，就是同时运行4个py命令，能够大幅减少运行时间，可能会干爆CPU内存，我的是非常小的城市区域，4线程全程吃内存30G左右，诸位量力而行

更多的参数，可以通过stackSentinel.py --help获取。

运行stackSentinel.py后将会在$WORK_PATH找到对应的文件夹，run_file和config，运行run_file中的所有文件，才算正式进行干涉图生成，很吃内存也很耗时间

建议是一步一步去做，防止中间出错或忽略掉一些提示信息，例如：

cd $WORK_PATH/run_files
bash run_01_unpack_topo_reference

你也可以在python中完成：

import os
os.system(directory+'/run_01_unpack_topo_reference')
print('run_01_unpack_topo_reference complete.')

directory为$WORK_PATH/run_files

所有的结果都被保存在$WORK_PATH中，可以实时地去查看，运行完所有的run_files后，我们所需的结果都在merged文件夹中，其中geom_reference里的文件记录了雷达坐标系和地心坐标系之间的转换关系，用于后续地理编码的操作，interferograms文件夹里即为各个日期对应的干涉结果：
.unw 解缠相位结果
filt_fine.int 滤波后干涉图
fine.int 原始干涉结果
.cor 相干性结果

注意ISCE生成的所有的vrt+xml+meta源文件形式的文件，都可以用mdx.py检查结果，
例如我查看某一日期对的滤波后干涉影像：

cd $WORK_PATH/merged/interferograms/20220226_20220310
mdx.py filt_fine.int

4. 使用MintPy计算时序InSAR形变序列

在获取了干涉序列后，最后一步就是用MintPy去获得最终的时序InSAR形变产品,首先在$WORK_PATH中创建一个mintpy文件夹，作为MintPy的工作文件夹：

cd $WORK_PATH
mkdir mintpy; cd mintpy
smallbaselineApp.py -g

smallbaselineApp.py是mintpy的核心，这里我们创建了mintpy文件夹，并使用smallbaselineApp.py -g生成了一个cfg参数配置文件，名字叫smallbaselineApp.cfg。

该文件中对所有参数都给出了超级详细的解释，按照需求进行修改即可，其他参数可以默认auto，不写在cfg文件中也可以的。这里给出我的cfg文件：

# vim: set filetype=cfg:
##------------------------ smallbaselineApp.cfg ------------------------##
########## computing resource configuration
mintpy.compute.maxMemory = 24 #分配的运算内存[float > 0.0], auto for 4, max memory to allocate in GB
mintpy.load.processor       = isce  #指定你的干涉处理平台[isce, aria, hyp3, gmtsar, snap, gamma, roipac, nisar], auto for isce
mintpy.load.updateMode      = no  #是否每次都是忽略已有的结果，重新生成结果[yes / no], auto for yes, skip re-loading if HDF5 files are complete

##---------for ISCE only:
mintpy.load.metaFile        = ../reference/IW*.xml  #指定IW对应的xml数据[path of common metadata file for the stack], i.e.: ./reference/IW1.xml, ./referenceShelve/data.dat
mintpy.load.baselineDir     = ../baselines  #基线文件夹[path of the baseline dir], i.e.: ./baselines
##---------interferogram stack:
mintpy.load.unwFile         = ../merged/interferograms/*/filt_fine.unw  #所有的解缠后相位文件unw[path pattern of unwrapped interferogram files]
mintpy.load.corFile         = ../merged/interferograms/*/filt_fine.cor  #[path pattern of spatial coherence       files]
mintpy.load.connCompFile    = ../merged/interferograms/*/filt_fine.unw.conncomp  #[path pattern of connected components    files], optional but recommended
##---------geometry用于地理编码的文件:
mintpy.load.demFile         = ../merged/geom_reference/hgt.rdr  #[path of DEM file]
mintpy.load.lookupYFile     = ../merged/geom_reference/lat.rdr  #[path of latitude /row   /y coordinate file], not required for geocoded data
mintpy.load.lookupXFile     = ../merged/geom_reference/lon.rdr  #[path of longitude/column/x coordinate file], not required for geocoded data
mintpy.load.incAngleFile    = ../merged/geom_reference/los.rdr  #[path of incidence angle file], optional but recommended
mintpy.load.azAngleFile     = ../merged/geom_reference/los.rdr  #[path of azimuth   angle file], optional
mintpy.load.shadowMaskFile  = ../merged/geom_reference/shadowMask.rdr  #[path of shadow mask file], optional but recommended
########## 4. correct_unwrap_error (optional)
## supported methods:
## a. phase_closure          - suitable for highly redundant network
## b. bridging               - suitable for regions separated by narrow decorrelated features, e.g. rivers, narrow water bodies
## c. bridging+phase_closure - recommended when there is a small percentage of errors left after bridging
mintpy.unwrapError.method          = bridging  #选择解缠误差改正方法[bridging / phase_closure / bridging+phase_closure / no], auto for no
mintpy.troposphericDelay.method = no  #是否进行大气校正，这里不进行大气校正[pyaps / height_correlation / gacos / no], auto for pyaps
########## 9.2 reference_date
## Reference all time-series to one date in time
## reference: Yunjun et al. (2019, section 4.9)
## no     - do not change the default reference date (1st date)
mintpy.reference.date = no   #指定参考日期，no为第一个日期[reference_date.txt / 20090214 / no], auto for reference_date.txt
mintpy.plot.maxMemory = 16  #分配给绘图的内存[float], auto for 4, max memory used by one call of view.py for plotting.

如果你是按照上述流程进行的stackSentinel生成，那么直接用这个参数是可以获得一个结果的。mintpy系统自带使用pyAPS进行大气校正，需要额外下载ERA5数据，这里比较麻烦就不赘述，后续可以再开一篇博客。

设置好smallbaselineApp.cfg配置文件后，直接运行：

smallbaselineApp.py smallbaselineApp.cfg

将会自动运行完所有步骤并获取结果，

mintpy获取的结果均用hdf5文件（.h5）进行储存，可以看到mintpy文件夹下生成了很多.h5文件，其中比较重要的列举一下：
timeseries.h5累计形变序列文件
timeseries_demErr.h5经过地形DEM误差改正后的累计形变序列文件
velocity.h5形变速率文件

\geo文件夹下储存了所有系统自动进行地理编码后的结果，如geo_timeseries_demErr.h5包含了各个日期对应的累积形变序列，储存于hdf5数据格式中

\pic文件夹包括了系统给出的一些可视化结果

值得一提的是，mintpy给出了大量的程序进行DIY操作，便于我们进行出图、提取数据等工作。

1. 首先，对于所有mintpy生成的h5文件，均可使用info.py查看其文件结构以及相关信息。如info.py timeseries_demErr.h5，这一功能类似于GDAL中的gdalinfo。
2. 所有的timeseries文件，均可以用tsview.py进行交互式可视化查看，如tsview.py timeseries_demErr.h5。
3. 非时间序列的 .h5文件，可以用view.py进行查看，view.py velocity.h5
4. save_gdal.py从.h5文件中提取栅格文件转为ENVI或Tiff格式
5. geocode.py对文件进行地理编码处理，可以实现rdr雷达坐标系下的栅格文件和WGS84坐标栅格文件的互相转化

以上所有的代码均可通过-h命令查看对应的参数，从而进行更加精细的DIY操作。

---

以上就是一个简单的全流程操作，可能还有没有说清楚的地方，或者有纰漏之处，请大家多多指正。附上刚刚学习的时候参考过的一些写的比较好的博客或有用的网站，欢迎交流。

参考

【INSAR形变监测】SENTINEL-1 TOPS INSAR全流程 (ISCE2+MINTPY)

isce2学习（三）下载哨兵数据、POD精轨数据、DEM数据、AUX_CAL辅助文件下载

mintpy另外给出了许多使用jupyter notebook写的step-by-step教程，对新手很友好，可以参见：

https://github.com/insarlab/MintPy-tutorial