写在前面：
本文工作旨在搭建开源sbas-InSAR环境，着重分享个人之前的开发经验。由于需要对ISCE源码做一些变动，所以源码采用编译安装的方式，但是很多依赖库还是基于conda安装。为了方便起见，以docker镜像的方式构建了代码仓库，覆盖处理过程中所有使用到的工具和脚本。镜像构建和程序运行指令总结见README。
另外，本文工作针对S1A数据处理即ISCE的topsStack模块，源于自己之前的一些数据处理经验和任务并行需求，改动之处多在细微之处，并未修改主要的处理算法。由于基础框架的处理步骤在官网和许多博客都有解读，这里不再赘述，着重讲述本文代码仓库的命令运行和方法变动。
Tips：
1. 在克隆代码仓库后构建镜像前，需要先修改config_files中的netrc和cdsapirc文件，在文件中写入自己的账户信息，分别用于DEM数据和大气校正数据的下载，详细配置方法见ISCE和PyAPS；
2. docker镜像已成功测试，熟悉docker环境的同学可通过docker build一键搭建，放心食用~；
3. 对于不想使用docker环境且自行配置环境编译代码的同学，关注代码仓库相较于ISCE-2.6.2源码仓库的改动即可；
4. 对于想继续使用conda环境的同学，如若想应用本文所述的优化方法，需要将代码仓库中修改过的py脚本以及编译后的文件替换原始conda环境中的对应文件。这里以我的conda环境为例，给出编译后的源码程序和需要替换的conda环境程序具体路径：
一、ISCE篇
本文代码仓库中的ISCE基础版本为2.6.2。该部分主要通过ISCE2中的topsStack模块，将下载好的S1A数据提取为SLC格式，并进行配准、干涉、滤波、相位解缠等处理，主要输出产品为用于时序分析的干涉图数据及其相关信息。处理流程图如下所示：
熟悉ISCE topsStack的同学应该很眼熟了。需要说明的是，流程图中的黑色文本框表示原ISCE中的处理方法，红色文本框意味着在原ISCE处理方法的基础上做了某些添加或修改。 首先，给出流程中需要用到的命令参数。
| 参数名称 | 参数说明 | 示例 |
| nasadem_directory | DEM存档数据路径 | /opt/data/nasaDEM |
| dem_directory | DEM输出数据存储路径 | /opt/data/descending/DEM |
| dem_bbox | DEM数据覆盖范围 | 17 20 -107 -98 |
| runfiles_path | 执行文件路径 | /opt/data/descending/process/run_files |
| working_directory | 工作路径 | /opt/data/descending/process |
| slc_directory | SAR影像存储路径 | /opt/data/descending/SLC |
| aux_directory | 辅助文件存储路径 | /opt/data/descending/AuxDir |
| orbit_directory | 轨道数据存档路径 | /opt/data/Orbits |
| slc_bbox | 研究区范围 | '19.3 19.5 -99.15 -98.86' |
| azimuth_looks | 方位向多视数 | 2 |
| range_looks | 距离向多视数 | 10 |
| reference_date | 主影像日期 | 20201206 |
| num_connections | 干涉图连接数 | 3 |
| filter_strength | 干涉图滤波权重因子 | 0.4 |
| swath_num | 要处理的条带数 | '3' |
| cc_thres | 干涉相干性阈值 | 0.35 |
| pwr_thres | 相对强度阈值 | 0.075 |
| num_process | run_files并行任务数 | 2 |
| num_process4topo | 地理编码并行任务数 | 1 |
| snr_misreg_threshold | 信噪比阈值 | 10 |
| num_overlap_connections | 叠掩干涉连接数 | 5 |
| esd_coherence_threshold | 相干性阈值 | 0.85 |
| nfft | 滤波FFT窗口尺寸 | 32 |
| subset | 干涉图是否裁剪 | 1 |
| parallel_num | run_files分段处理 | 4 |
大部分参数都来自于stackSentinel.py程序，其中新增的一些参数如“cc_thres”，“pwr_thres”会在后续每步程序执行时进行说明。接下来，我将以墨西哥城2021年33景降轨S1A数据和上表中的参数为例，给出流程图中每个步骤需要执行的脚本命令，并解释其相对原始方法的改动。
1. 构建镜像并启动容器
在系统终端中运行下面的命令构建镜像并启动容器。
docker build --rm --force-rm -t insar_process:v1.0 -f Dockerfile .
xhost +
docker run -itd --privileged=true --net=host --shm-size 4g --gpus=all --name insar -e DISPLAY=$DISPLAY -v /media/comma/InSAR2/Mexico/:/opt/data insar_process:v1.0
docker exec -it insar /bin/bash2. DEM获取及处理
在容器内运行下面的命令以获取干涉处理过程中需要的DEM数据。
mkdir -vp $dem_directory && cd $dem_directory && \
dem.py -a stitch -b $dem_bbox -r -t nasadem -c -l -d $nasadem_directory && \
fixImageXml.py -i *.dem.wgs84 -f && \
rm -f $nasadem_directory/demLat* $nasadem_directory/*.hgt这里使用的是ISCE中的“dem.py”方法。由于dem下载链接有时候不稳定，包括后面要使用的轨道数据也是一样，所以自己爬了所有的nasaDEM数据和精密轨道数据，分别存放在本地数据存档路径“nasadem_directory”和“orbit_directory”中，在该命令中启用“-l”和“-d”选项跳过原始DEM的下载阶段，读取本地的DEM数据并处理生成ISCE格式的研究区dem。
3. 可执行文件生成
在容器内运行下面的命令以生成InSAR数据处理需要运行的可执行文件。
mkdir -vp $working_directory $slc_directory $aux_directory $orbit_directory && \
dem=`find $dem_directory -name "*.dem.wgs84"` && \
stackSentinel.py --working_directory="$working_directory" \
--slc_directory="$slc_directory" \
--dem="$dem" --aux_directory="$aux_directory" --orbit_directory="$orbit_directory" \
--bbox="$slc_bbox" --azimuth_looks="$azimuth_looks" --range_looks="$range_looks" \
--reference_date="$reference_date" --num_connections="$num_connections" \
--filter_strength="$filter_strength" --nfft="$nfft" --swath_num="$swath_num" \
--cc_thres="$cc_thres" --pwr_thres="$pwr_thres" --num_process="$num_process" \
--num_process4topo="$num_process4topo" \
--snr_misreg_threshold="$snr_misreg_threshold" \
--num_overlap_connections="$num_overlap_connections" \
--esd_coherence_threshold="$esd_coherence_threshold" \
--useGPU这部分内容基于ISCE2中的“stackSentinel.py”方法，对应本文的第一张表，具体有以下一些修改：
（1）流程和并行参数设置
参数：--num_process
这个参数其实是“stackSentinel.py”本身就有的参数，用于设置可执行文件中同时运行的任务数量。以干涉图滤波步骤的可执行文件run_15_filter_coherence为例，如果不设置该参数即保持默认值为1，程序默认输出的可执行文件内容格式如下：
SentinelWrapper.py -c /opt/data/descending/process/configs/config_igram_filt_coh_20201206_20201218
SentinelWrapper.py -c /opt/data/descending/process/configs/config_igram_filt_coh_20201206_20201230
SentinelWrapper.py -c /opt/data/descending/process/configs/config_igram_filt_coh_20201206_20210111
SentinelWrapper.py -c /opt/data/descending/process/configs/config_igram_filt_coh_20201218_20201230即在默认情况下运行可执行文件时，终端以串行的方式依次执行每一行命令。如果启用该并行参数，如设置“num_process”的值为2，那么程序输出的可执行文件内容会变成下面这样：
SentinelWrapper.py -c /opt/data/descending/process/configs/config_igram_filt_coh_20201206_20201218 &
SentinelWrapper.py -c /opt/data/descending/process/configs/config_igram_filt_coh_20201206_20201230 &
wait
SentinelWrapper.py -c /opt/data/descending/process/configs/config_igram_filt_coh_20201206_20210111 &
SentinelWrapper.py -c /opt/data/descending/process/configs/config_igram_filt_coh_20201218_20201230 &
wait在这种情况下运行可执行文件，终端中每次会同时运行两行命令。显然，在计算资源足够的情况下，“num_process”的值越大，可执行文件运行的时间就越短。
在ISCE源码中，启用该参数时所有步骤的可执行文件都会被修改。但是之前我在测试的时候遇到一些问题，比如有些步骤在同时运行多个任务时，出现“挂起”的异常现象，即任务完成后，终端命令不会自行结束，而一直停留在运行状态。还有一种异常现象是，某些步骤通过这种任务并行的方式，输出文件莫名丢失了一部分，进而导致后续任务找不到需要的文件而出错。
当时是在云计算环境中处理时遇到上面的问题，无法判断本地主机运行时是否会出现同样的问题，最终对“stackSentinel.py”、“Stack.py”中的这部分代码进行了修改，排除了相关问题。修改后启用该参数只会对部分步骤设置任务并行。
（2）干涉图滤波
参数：--nfft
这个参数是新增参数，用来调整滤波中的窗口大小。这里并没有修改ISCE中的干涉图滤波算法，而是将原本算法中的该参数暴露出来。另外，原方法中只进行了一次滤波，本文仓库中默认对干涉图进行了两次滤波以提升相干性。这里涉及到修改的代码有“isce2-2.6.2/contrib/stack/topsStack/FilterAndCoherence.py”以及“isce2-2.6.2/components/mroipac/filter”中的部分源码。
（3）相位解缠
参数：--cc_thres，--pwr_thres
这两个参数都是新增参数，用于在相位解缠时进行相干性和强度值掩膜。ISCE源码中本身没有为Snaphu解缠提供掩膜的接口，通过相干性掩膜噪声区域、强度值掩膜水体，可以提高相位解缠的精度，并且为后续的时序分析提供一个掩膜文件，相比单一的相干性掩膜可以更直接有效地掩膜研究区中的水域。这里涉及到修改的代码有“isce2-2.6.2/contrib/stack/topsStack/unwrap.py”以及“isce2-2.6.2/contrib/Snaphu”中的部分源码。另外，原ISCE仓库中应用的Snaphu-1.4.2，本文代码仓库中的Snaphu模块是在Snaphu-2.0.6的基础上进行了一些修改。
（4）GPU处理模块
参数：--useGPU
这个参数是“stackSentinel.py”本身具有的参数，用来启用相关步骤的GPU处理，默认不设置该选项完全依赖CPU计算。
我不太清楚当前的ISCE版本是否对这方面内容做了更新，但是两年前，当时的ISCE官方仓库中有GPUgeo2rdr，GPUresampslc，GPUtopozero这三个GPU模块，但是通过编译后好像只有GPUgeo2rdr模块可用，至于其他两个模块是编译失败还是运行时抛出异常，详细情况我也记不清了。我当时也是在官方实现的基础上做了些修改，使这三个模块都编译可用，但是印象只有GPUtopozero模块可以为run_01的处理带来比较明显的提速。原谅当时的我对于cuda也是一知半解，现在更是记不清修改逻辑，最后总归是测试通过了，GPUtopozero的加速还算有些意义。感兴趣的同学可以查看“isce2-2.6.2/components/zerodop”中的这部分代码并和官方实现进行对照。
4. 可执行文件分割
在容器内运行下面的命令以分割可执行文件，并为可执行文件添加权限。
cd $runfiles_path && split_file.sh $parallel_num
chmod 777 -R $runfiles_path/*“parallel_num”也是一个并行处理参数，这里需要和前面的“num_process”参数区分下。 “num_process”作为“stackSentinel.py”程序的输入参数，用于控制每个可执行文件内以串行还是并行的方式处理任务。而“parallel_num”作为“split_file.sh”脚本的输入参数，是将某些步骤中，原本ISCE生成的单个可执行文件分割为多个。
同样以可执行文件run_15_filter_coherence为例。假设文件内有93行命令，以“split_file.sh”按“parallel_num”等于4进行分割，可能会将其分割为run_15_filter_coherence_1、run_15_filter_coherence_2、run_15_filter_coherence_3、run_15_filter_coherence_4、run_15_filter_coherence_5等多个文件，这几个文件中的命令行总数为93。后续在运行可执行文件时就不再运行原本的run_15_filter_coherence，而是具有编号后缀的可执行文件。
这个脚本设计的初衷是针对云上计算时，可以多节点同时处理一个步骤的可执行文件。在本地环境这个参数选项则看上去有些多余，但是如果你的本地工作站性能足够强大，你依然可以设置“$parallel_num”参数，在后续运行可执行文件时，可以打开多个终端并同时运行不同编号后缀的可执行文件，实现手动并行加速。如果不需要的话，将“$parallel_num”设为1并顺序运行原本的可执行文件就好了。
5. 可执行文件运行
类似ISCE官方的处理过程，这部分是在容器内依次运行前面生成的可执行文件。带有for循环的是流程中对可执行文件做了分割处理的步骤，默认通过for循环串行处理。你也可以通过上面提到的多终端同时运行方法手动执行。
${runfiles_path}/run_01*
for ((idx=1; idx<=$(($(ls ${runfiles_path}/run_02* 2>/dev/null | wc -l)-1)); idx++)); do ${runfiles_path}/run_02*_${idx}; done
${runfiles_path}/run_03*
${runfiles_path}/run_04*
for ((idx=1; idx<=$(($(ls ${runfiles_path}/run_05* 2>/dev/null | wc -l)-1)); idx++)); do ${runfiles_path}/run_05*_${idx}; done
for ((idx=1; idx<=$(($(ls ${runfiles_path}/run_06* 2>/dev/null | wc -l)-1)); idx++)); do ${runfiles_path}/run_06*_${idx}; done
for ((idx=1; idx<=$(($(ls ${runfiles_path}/run_07* 2>/dev/null | wc -l)-1)); idx++)); do ${runfiles_path}/run_07*_${idx}; done
${runfiles_path}/run_08*
for ((idx=1; idx<=$(($(ls ${runfiles_path}/run_09* 2>/dev/null | wc -l)-1)); idx++)); do ${runfiles_path}/run_09*_${idx}; done
for ((idx=1; idx<=$(($(ls ${runfiles_path}/run_10* 2>/dev/null | wc -l)-1)); idx++)); do ${runfiles_path}/run_10*_${idx}; done
${runfiles_path}/run_11*
${runfiles_path}/run_12*
for ((idx=1; idx<=$(($(ls ${runfiles_path}/run_13* 2>/dev/null | wc -l)-1)); idx++)); do ${runfiles_path}/run_13*_${idx}; done
for ((idx=1; idx<=$(($(ls ${runfiles_path}/run_14* 2>/dev/null | wc -l)-1)); idx++)); do ${runfiles_path}/run_14*_${idx}; done
if [ "${subset}" -eq 1 ]; then subset_isce.py -d ${working_directory} -b "$slc_bbox" -r ${range_looks} -a ${azimuth_looks}; fi
for ((idx=1; idx<=$(($(ls ${runfiles_path}/run_15* 2>/dev/null | wc -l)-1)); idx++)); do ${runfiles_path}/run_15*_${idx}; done
pre_unwrap.sh ${working_directory} ${range_looks} ${azimuth_looks}
for ((idx=1; idx<=$(($(ls ${runfiles_path}/run_16* 2>/dev/null | wc -l)-1)); idx++)); do ${runfiles_path}/run_16*_${idx}; done
unw=`find ${working_directory} -name "filt_fine.unw" | sed -n "1p"` && \
gdal_translate ${unw} -b 1 -of ISCE ${working_directory}/magnitude && \
generate_mask.py ${working_directory}/magnitude --min 1 -o ${working_directory}/waterMask.h5 && \
rm ${working_directory}/magnitude*相较于原始的ISCE topsStack方法只需要依次运行16个可执行文件，上面的命令中穿插了额外的命令，如在run_14和run_15之间、run_15和run_16之间、run_16之后，这些额外的命令对应于流程图中的红色文本框。
在run_14和run_15之间，即干涉图滤波之前，我们通过“subset_isce.py”对前面生成的干涉图进行裁剪，裁剪范围和“stackSentinel.py”中输入的“slc_bbox”保持一致。这一步是可选项，因为默认的ISCE会保留输入“slc_bbox”范围内的整个Burst。这个步骤设计初衷是针对研究小范围区域，通过提前裁剪干涉图使得相位解缠的速度大幅增加。当时我最开始的想法是裁剪SLC，但是后来没能找到具体方法，索性就在这一步裁剪了干涉图。
在run_15和run_16之间，即相位解缠之前，通过“pre_unwrap.sh”生成相位解缠需要读取的相关掩膜文件，平均强度图ave.mli和平均相干图ave.cor。
在run_16后，通过MintPy模块中的“generate_mask.py”等命令生成掩膜文件waterMask.h5，用于后续Mintpy时序分析。
6. 一些经验
（1）“stackSentinel.py”运行时研究区范围“slc_bbox”过小可能导致报错IndexError: list index out of range；
（2）run_08运行完之后，可以使用下面命令查看misreg文件夹中的配准精度；
grep std ./misreg/range/pairs/*/* # 距离向干涉图配准精度
grep std ./misreg/azimuth/pairs/*/* # 方位向干涉图配准精度
grep 0. ./misreg/range/dates/* # 距离向时序配准精度
grep 0. ./misreg/azimuth/* # 方位向时序配准精度按我自己的理解，方位向干涉图配准精度可能会出现某行配准精度为0意味着配准失败，过多干涉图配准失败的话可以返回调整配准参数（esd_coherence_threshold、snr_misreg_threshold和num_overlap_connections），并返回运行run_07和run_08重新配准。上述理解有误的话也请大家批评指正哈。
（3）需要修改流程中相关的处理参数并重新运行相关步骤时，可以使用下面的命令直接修改配置文件的参数，不需要返回执行stackSentinel命令并从头开始处理。
cd ./configs # 进入configs文件夹，里面有全部执行步骤的配置文件
sed -i "s/strength : 0.5/strength : 0.8/g" config_igram_filt_coh_* # 以修改滤波器系数为例，运行该命令将所有滤波配置文件中的strength : 0.5修改为strength : 0.8
cd run_files && ./run_15_filter_coherence # 修改之后重新从滤波步骤开始运行（4）确认ISCE处理完无误之后，需要baselines、merged、reference、secondarys四个文件夹做时序分析，剩余过程文件可以删掉腾出空间。
二、MintPy篇
本文代码仓库克隆了MintPy-1.6.1。关于MintPy的配置和运行方法，官方手册已经相当完备了，这里简单介绍下处理过程。
当具备了上述ISCE输出的baselines、merged、reference、secondarys四个文件夹后，我们在同级路径下（其他路径的话需要在配置文件编辑）创建一个新的mintpy文件夹并进入，然后将编辑好的配置文件（如smallbaselineSentinel.txt）也一同放进去并运行smallbaselineApp.py就OK了。这里有一点需要注意，配置文件名中需要包含Sentinel字段。
mkdir -vp ${working_directory}/mintpy && cd ${working_directory}/mintpy
smallbaselineApp.py smallbaselineSentinel.txtmintpy的处理参数均在配置文件中设置，官方给出的配置文件中包含每个参数的详细说明和引用算法， 这里列一下常用的处理参数。
mintpy.compute.maxMemory = 16 # 程序运行的最大内存
mintpy.compute.cluster = auto #[local / slurm / pbs / lsf / none], auto for none, 设为local可开启多线程
mintpy.compute.numWorker = auto #[int > 1 / all / num%], auto for 4 (local) or 40 (slurm / pbs / lsf), num of workers
mintpy.subset.yx = auto #[y0:y1,x0:x1 / no], auto for no 根据xy坐标裁剪范围
mintpy.subset.lalo = auto #[S:N,W:E / no], auto for no 根据地理坐标裁剪范围
mintpy.network.coherenceBased = yes # 开启根据相干性进行网络校正
mintpy.network.minCoherence = 0.4 # 去除相干性低于0.4的干涉图
mintpy.reference.yx = auto #[257,151 / auto] 输入参考点的xy坐标
mintpy.reference.lalo = 41.20373,121.99064 #[31.8,130.8 / auto] 输入参考点的地理坐标
mintpy.unwrapError.method = phase_closure # 选择相位解缠误差校正的方法
mintpy.networkInversion.minTempCoh = 0.6 # 掩膜时间相干性低于0.6的像元
mintpy.deramp = quadratic # 去除相位趋势面的模型参数
mintpy.topographicResidual = yes # 开启去除DEM误差
mintpy.topographicResidual.polyOrder = 3 # DEM误差模型参数
mintpy.timeFunc.polynomial = auto #[int >= 0]，auto for 1 形变模型参数，默认1为线性形变速率
mintpy运行结束后，过程文件主要以图片格式保存在pic文件夹中，数据文件以h5格式保存在geo文件夹中。
MintPy仓库包含很多实用的方法，希望大家不要局限于它的主程序运行，多多分析它的源码方法。这里列出几个查看和导出数据文件的方法，比如“info.py”查看数据信息，“view.py”图形显示数据，“save_*.py”命令导出数据，下面是一些常用命令和查看的形变速率结果：
view.py ./geo/geo_velocity.h5 velocity -v -10 10 -u mm # 查看形变速率，颜色栏阈值为[-10,10]，单位为mm
tsview.py ./geo/geo_timeseries_ERA5_ramp_demErr.h5 -v -20 20 # 查看时序形变结果,可以通过设置--ref-date参数将参考日期设为首个日期
save_gdal.py ./geo/geo_velocity.h5 # 形变速率导出为GTiff格式
geocode.py ifgramStack.h5 -d unwrapPhase -l geometryRadar.h5
save_gdal.py geo_unwrapPhase.h5 -d unwrapPhase-20220608_20220831 --of ISCE
三、MSBAS篇
MSBAS同样是一款出色的开源时序InSAR数据处理软件，相较于MintPy致力于去除各种相位残差和干涉图优化，它更侧重于实现一维至多维时序形变求解，其主页讲述了三个版本MSBASv3、MSBASv6和MSBASv10的迭代过程。本文代码仓库克隆了MSBASv6，对MSBAS不熟悉的同学可以阅读论文并下载软件学习。为了防止大家走丢，我在代码仓库中一并放入了MSBASv3，其中包含源码、样例数据、论文和使用说明。
不同版本MSBAS主程序运行的头文件参数可能有所不同，但是整体上处理过程是类似的。详细的程序运行方法大家参考下我的运行脚本“run_msbas.sh”。
当你通过上面的ISCE+MintPy处理得到同一区域升降轨的数据结果，那么就可以通过MSBAS将LOS向形变分解为二维形变了（三维形变解算我还没了解过~）。MSBAS程序运行的核心思路大概分为三步，数据准备、写入参数文件、运行。“run_msbas.sh”中大部分内容都是在从MintPy处理的时序结果中导出MSBAS所需的数据，最后通过“msbas header.txt”运行核心的时序形变解算，最后还可以使用“msbas_extract”命令以文本的方式提取某一形变点的时间序列。下面是一个头文件header.txt样例。
FORMAT=2
FILE_SIZE=1033,1408
C_FLAG=0
R_FLAG=0
I_FLAG=0
SET=0,000000,-12.262834697288602,44.21425,asc.txt
SET=0,000000,-167.7257087622666,44.197063,dsc.txt下图是结合墨西哥城升降轨数据计算的UD向和EW向形变速率图。
四、GMT篇
这里是GMT作图专区，误认GMTSAR的同学绕行哈~
如果觉得MintPy自动保存的图片不满足需求，或是在寻求比Arcgis省时的批量出图方法，这里就要力荐一波GMT了。需要注意一下，镜像中通过apt安装的GMT版本为6.0.0，而代码仓库中的制图脚本为GMT5格式。通常情况下GMT6可以运行GMT5语法的脚本，如果遇到问题的话大家可以自行在系统上安装GMT5或是将脚本语法升级为GMT6。
由于GMT具有自己复杂的语法，又有很强的定制化属性，本文给出两个制图脚本“plot_deformation.sh”、“plotp_deformation.sh”供大家参考，用于从MintPy的输出结果中绘制形变速率图和时序形变图，建议初学者对照GMT官方手册理解脚本命令。这两个绘图脚本的不同之处在于，“plotp_deformation.sh”需要一个由“select_pslike.py”脚本提取的“geo_tcp”文件，只绘制提取后的形变点，感兴趣的同学可以自行查看脚本。这里放一张“plotp_deformation.sh”绘制的形变速率图。也可在脚本中添加“gmt psimage”命令添加DEM或光学影像底图，可以说GMT制图效果取决于大伙的需求和对GMT语法的掌握程度。
写在最后：
以上内容算是对之前某段时期工作的部分总结了，就先写到这儿吧。当时感觉差强人意的内容，如今看来也乏善可陈。时隔两年，整理不易，希望大家能以批判的眼光看待，也希望能为新人开发者提供一些解决问题的思路。
扫一扫
6181
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
