本文参照 ENVI SARscape入门教程 进行实际研究的操作
一、数据准备
1.从ASF网站下载相应地区的SLC数据：https://search.asf.alaska.edu/
首页如上图：其中1为区域选取方式，选择Geographic Search；2为卫星数据，选择Sentinel-1即哨兵1；3为条件过滤器；4为区域绘图工具
设置好区域后打开条件过滤器，设置日期等，其中file type选择SLC数据，direction选择ascending（升轨），subtype选择SA
检索出结果后选择相应的数据，其中红框为当前影像范围，点击购物车按钮将影像加入购物车，可以在购物车中直接下载，也可以批量下载
在进行分析时提示最少应该有16景图像，其中一幅是超级主影像，由于补充下载图像时ASF暂时无法下载了，因此选择欧空局渠道下载：https://scihub.copernicus.eu/dhus/#/forgot-password
欧空局网页无法通过区域范围筛选数据，不太简便，因此我选择先使用更简便的ASF进行筛选，筛选到具体的影像后根据名称到欧空局网站下载
点击需要的SLC影像名称，复制到欧空局网站的搜索栏中，即可搜索到该影像
此时我们可以看到数据处于Offline状态，即未上线，这时候可以将其加入到购物车中，再点击下载，出现下面的提示即说明系统此时开始往服务器上传数据，等待一段时间后可以下载（真的是一段时间，有可能几个小时有可能几天）
在等待欧空局上线数据时，又尝试了一下在ASF中使用其提供的python文件批量下载的方法：
在ASF的购物车界面，点击下方的Data Download——Download Python Script.py，将会下载一个名为download-all+时间的.py文件
将该文件放到要保存的目录下，双击其即开始运行，输入用户名和密码
同样也是由于网络原因无法下载，但在网络可以的情况下应该能够批量下载
2.从ASF网站下载相应的精密轨道文件（非必须）：https://s1qc.asf.alaska.edu/aux_poeorb/
打开后即可看见轨道文件列表，按ctrl+F打开搜索，输入年月日搜索文件
右键点击需要下载的文件，点击链接另存为即可保存，如果无法保存说明没有登录ASF，需要到ASF网站登录账号。
保存好后需要将其作为文件夹等待读取，此处参照SARscape存放sentinel-1精密轨道数据并自动读取_Vigo*GIS_RS的博客-CSDN博客
在文件夹中新建一个二级目录，Sentinel_Orbit \ AUX_POEORB，用于存放常见的轨道数据，也就是POEORB数据，同时在Sentinel_Orbit里面新建其他三个格式的轨道数据（空文件夹）
在AUX_POEORB文件夹中按照按照卫星的型号、年份、月份分别建立文件夹保存轨道文件
数据准备部分到这里就差不多了，剩余的DEM部分在后续操作过程中再获取
二、软件处理
本文使用了ENVI5.6以及相应的SARScape模块（安装到ENVI中中文名是高级雷达处理模块）
1.系统参数设置
首先打开Preferences----Perference specific设置，点击左上角的Load Preferences——Sentinel TOPSAR
设置好默认参数后设置轨道文件的读取路径：点击Preferences----Perference common
2.SLC（SAR）数据预处理：导入，裁剪
在SARScrape中点击Import Data——SAR Spaceborne——Single Sensor——SENTINEL-1，打开导入窗口
将解压好的SLC文件中的manifest.safe文件导入（也可以直接导入zip）
在同一窗口中打开Optional Files，将提前设置好的裁剪区域shp文件导入
将Parameters窗口中的极化方式选择Copolarization only，输出VV极化方式
在output files窗口中我将输出文件根据日期进行了名称修改以更好地分辨，点击右侧按钮修改新路径，点击Exec开始输出
注：由于ENVI没有保存整个工程文件的功能，因此如果裁剪文件处理好后想要中途退出再打开是没有原来裁剪好后的文件的，这时候打开ENVI，点击文件-打开-找到裁剪保存文件夹，将其中的.file文件选中打开即可，这些文件通常为以下格式：sentinel1_11_20200413_103324920_IW_A_VV_slc_list_pwr
3.生成连接图
裁剪过程花了非常久，接下来就是生成连接图的步骤
在SARScape工具目录下点击Interferometric Stacking—SBAS—“1-Connection Graph”
将左侧工具栏中的所有SAR图像拖入Input Files窗口
在参数设置中设置Max normal Baseline为45%~50
MaxTemporalBaseline设置为500（为了让所有数据都参与计算）
在Other Parameters中的Delaunay 3D：是否启动 3D 解缠，设置为 True（需要将Degree of Redundancy设置为High）
Redundacy Criteria设置为Min normal baseline
时间基线阈值主要根据区域条件来设置，干燥地区可设置 500-800 天。在潮湿/植被浓密地区，需要设置小一些，避免出现无相干数据对。
设置好output files路径和名称后点击Exac即开始输出
建立连接图这一步，会生成整个处理的根目录（*_SBAS_processing），包括工作目录以及 SBAS 每一步处理的文件夹，包括下面几项内容：
– 用户可以查看的中间文件，做一些调整和编辑；
– SBAS 处理的工作目录记录所做的处理，如果中途停止的话，可以从停止的地方继续处理；
– 所用的参数列表：work_parameters.sml，该文件记录了 SBAS 处理中所设置的参数，用户可以随时查看某一步的参数设置。
– 连接图文本文件。
– 辅助文件：auxiliary.xml.记录进行的每一个处理，以及数据索引，后面的处理都是用这个文件作为输入文件。
这一步结束后，会自动加载多视后的超级主影像强度数据，通过浏览该数据，可以了解研究区的情况，该强度数据可用于选择研究区的时候选择斜距数据的区域范围。
以下是生成的时间和位置的连接图、以及能做 3D 解缠的子像对。在输出的图中，像对用线连接，输入的数据用点表示，被丢弃的图像标为红色、黄色的是超级主影像，提供了两种形式的绘图：时间和相对于主影像的位置图，时间和空间基线图。
时间和位置的连接图：有助于了解看数据有没很好的配对：每个时相都和其他时相建立了连接，连接的也比较均匀，理想的连接是每个图像至少有 5 个连接，如果有一组连接不能和大多数像对连到一起，这一组的所有像对将丢弃。
4.干涉工作流
对所有的配对的干涉像对进行干涉处理，从相干性生成，去平、滤波和相位解缠，所有的数据对都配准到超级主影像上，为下一步轨道精炼和重去平，以及 SBAS 反演做好数据准备。包括：
·干涉图生成 Interferogram Generation
·干涉图去平 Interferogram Flattening
·自适应滤波和相干系数生成 Adaptive Filter and Coherence Generation
·相位解缠 Phase Unwrapping
生成一系列解缠之后的相位图。所有的干涉图最终都与超级主影像进行了配准，为下一步轨道精炼、重去平以及 SBAS 的反演做准备。
步骤：
（1）打开/SARscape/Interferometric Stacking/SBAS/2 - Interferometric Process
（2）在 Input Files 面板中，选择前面连接图生成文件中的 auxiliary.sml
（3）在 DEM/Cartographic System 面板中，选择 DEM文件，建议选择使用SRTM-1 V3版的DEM数据，参照CSDN论坛中Vigo*GIS_RS的文章(54条消息) SARscape手动下载30mDEM（SRTM1 V3）切片数据-[EC: 40008]_envi中sar的dem下载_Vigo*GIS_RS的博客-CSDN博客
（4）在 Parameters 面板中，选择 Principal Parameters，设置以下参数：
•Range Looks 和 Azimuth Looks：多视视数，与超级主影像一致，在打开auxiliary.sml后会弹出对话框显示相应参数。
•Rebuild All：False
•Coregistration With DEM：True
•Unwrapping Method Type：解缠方法，有三种方法供选择，如果上一步选择了 3D 解
缠默认为 Delaunay MCF。
•Unwrapping 3D：设置为 True，使用 3D 解缠。
•Unwrapping Decomposition Level：解缠分解等级，默认 1。
•Unwrapping Coherence Threshold：解缠相关系数阈值，设置为 0.35.
•Filtering method：滤波方法，有三种方法供选择，默认 Goldstein。
点击Exac开始运行，运行需要两个小时，需要预留好时间
注：这一步是 SBAS 处理中运行时间最长的一步。这一步做完后，会生成一系列结果，包括各个像对去平和滤波后的干涉图(*_fint)，相干系数图 (*_cc) 、解缠结果 (*_upha)和斜距下的强度数据（*_slant_pwr），干涉处理的文件夹在根目录下面( interferogram_stacking)，包含一个索引文件（meta file），能打开所有的干涉处理的结果。生成一个 “interf_tiff”的文件夹，包含所有结果的 TIFF 格式结果，存放在 interferogram_stacking 文件夹下面，便于可视化查看。 查看解缠结果（IS_upha_list_meta）和相应的相干性图(IS_cc_list_meta)，看有无相干性低和解缠结果不好的像对来进行下一步的连接图编辑。
干涉工作流花费的时间极长极长，其实整个SBA-InSAR的过程均需要花较长的时间，在教程中所说的需要2.5小时是有所偏差的，这里给出一份个各流程的时间予以参考：
1.导入数据——每份数据约10-20分钟，如果需要裁剪的话需要更长时间
2.生成连接图——很快，约几分钟
3.干涉工作流——60小时，配置好的电脑快一些
4.第一次反演——50小时
5.第二次反演——25小时
6.地理编码——6小时
5.‘连接图编辑
依次查看上一步生成的各个像对的相干性图（_cc）和解缠结果图（upha），若有相干性低的，就用连接图编辑工具对该像对移除。
可以直接查看…\interferogram_stacking\ interf_tiff 中的 Tiff 格式文件。或者在 ENVI 中打开解缠结果（IS_upha_list_meta）、相干性图(IS_cc_list_meta)、滤波后干涉图（IS_fint_list_meta），在 ENVI 中浏览这些结果
由于每个人的经验不一样，判断需要移除的数据对也会不一样。这里直接贴上来自教程的剔除像对常用标准：
1．如果干涉处理时候用了 3D 解缠的方法，建议保留（得到的结果含有 3D_del_upha 结
果的说明进行的是 3D 解缠）；
2．相干系数图、滤波后的干涉图、相位解缠图，只要其中一个结果较好，建议保留；
3．得到的结果含有明显的大气相位，建议剔除；
步骤：
（1）打开/SARscape/Interferometric Stacking/Stacking Tools/Edit Connection Graph
在 Auxiliary file 选项中选择我们之前得到的auxiliary.sml文件，会自动识别出像对。分别在左侧和右侧面板选择要移除的数据对，左侧为m文件，右侧为s文件，例如结果文件为IS_xxxxxx_m_18_xxxxxx_s_14_*（其中xxx为时间）的相干系数低，解缠效果也不好，则在左侧M面板中选择 ID[18]数据，右侧S面板中中选择 ID[14] 数据，点击 Remove Pair 移除，重复上面的操作，直到移除掉所有需要移除的数据对。也有Automatic选项自动移除。
6.轨道精炼和重去平
这一步是估算和去除残余的恒定相位和解缠后还存在的相位坡道，虽然在SARScape5.6.2版本已经整合在了干涉工作流步骤中，但是那一步没有生成解缠结果图，不方便选择GCP点，因此我还是单独进行。
步骤：
（1）打开SARscape/Interferometric Stacking/SBAS/3 - Refinement and ReFlattening
（2）首先在 Auxiliary file 选项中选择我们之前得到的auxiliary.sml文件，之后点击GCP File选项旁的望远镜图标开始创建GCP文件
在打开的窗口中我们有三个文件需要选择，其中：
Input File文件需要选择一个数据对的解缠结果文件，文件名格式是IS_20160302_m_18_20160202_s_14_upha，这个是必选项
DEM File文件需要选择我们之前得到的DEM文件，这个可选可不选
Reference File需要选择一个参考依据，可以选择相干系数图、干涉图等，如这里选择一个去平和滤波后的干涉图（后缀为fint），这个可选可不选
注意这里有一个大坑：选择的文件夹千万不要错了，要选择work\work_interferogram_stacking中的不带后缀的相关文件才行，要不然选择好GCP点之后保存的GCP点xml文件会瞬间消失.......
选择好文件后点击OK，即进入GCP点选择过程，官方给出的点位选择标准是：
1、没有残余地形条纹；
2、没有形变条纹，远离形变区域，除非已知这个点的形变速率；
3、没有相位跃变，如果 GCP 点位于一个孤立相位上，并且解缠的值非常差，这个位置可能是斜坡相位（phase ramp）的一部分，那么选择的这个 GCP 是不对的；
反映在具体的upha图中是这样的
1处的图像均匀，无噪声无干扰无条纹，选择此处作为GCP点
2处图像有较多噪点，不予选择
（3）建议均匀选择，覆盖全图20-30点位及以上，选择好后设置好坐标系统（通常根据DEM默认为WGS84），在export中设置好输出路径，点击完成即可
（4）回到 Refinement and Re-Flattening 面板中。在 Input Files 面板中会自动读取上一步得到的控制点文件，设置好DEM文件
单击 Parameters面板，选择 Principal Parameters 选项，设置以下参数：
Rebuild All：False
Refinement method：轨道精炼方法，选择Polynomial Refinement，也可以选择Automatic Refinement即自动精炼
Refinement Residual Phase Polynomial Degree：精炼残余相位多项式次数，选择默认3
点击Exac进行轨道精炼和重去平
7.第一步反演（STEP 1 ）
这一步是 SBAS 反演的核心，第一次估算形变速率和残余地形。这一步也会做二次解缠用来对输入的干涉图进行优化，以便进行下一步处理，程序可以在二次解缠之前停止一下，让用户能分析第一次的结果，在 inversion 文件夹中。
步骤：打开SARscape/Interferometric Stacking/SBAS/4 - Inversion: First Step
在Input Files窗口中选择我们之前得到的auxiliary.sml文件
在Parameters窗口中设置参数：
（1）Rebuild All：选择False
（2）Product Coherence Threshold：相关系数阈值，当低于这个阈值的像素将以 NaN 输出。这里设置 0.35
（3）Allow Disconnected Blocks：是否丢弃孤立成像数据，默认为 False
（4）Model Type：Linear
（5）Stop Before Unwrapping：是否在二次解缠前停止，这里如果你不想中途查看一下的话注意要选择False，否则直接进行STEP2的时候会报错
（6）Unwrapping Method Type：解缠方法，有三种方法供选择，如果选择了 3D 解缠默认为 Delaunay MCF
（7）Unwrapping 3D：设置为 True，使用 3D 解缠
（8）Unwrapping Decomposition Level：解缠分解等级，默认 1
（9）Unwrapping Coherence Threshold：解缠相干系数阈值，设置为 0.35.
设置好后点击Exac即开始第一步反演
8.第二步反演（STEP2）
这一步的核心是计算时间序列上的位移，在第一步得到的形变速率(_disp_first)基础上，进行定制的大气滤波，从而估算和去除大气相位，得到更加纯净的时间序列上的最终位移结果。大气高通、大气低通两个选项，对大气影响进行估计，最后每个时间都从测量的位移中减去这些大气部分。在这里“Displacement GCP file”是用来去除还有残余的相位或相位斜坡，在这里，使用轨道精炼的时候所用的轨道控制点。
步骤：
打开SARscape/Interferometric Stacking/SBAS/5 - Inversion: Second Step
(1)在Input Files窗口中选择我们之前得到的auxiliary.sml文件
(2)在 Refinement GCP file 选项中选择我们前面轨道精炼和重配平阶段制作的控制点文件
(3)单击 Parameters 面板，选择 Principal Parameters 选项，设置以下参数：
Product Coherence Threshold：相关系数阈值，当低于这个阈值的像素将以 NaN 输出。这里设置 0.35
Interpol Disconnected Blocks：当时间序列数据不存在的部分，可以通过插值方法估算形变，默认False（如果前面处理中设置了 Allow Disconnected Blocks=True 时，这里可以设置为 True）
Atmosphere Low Pass Size：输入以米为单位的窗口大小。应用于空间分布相关滤波。设置为 1200
Atmosphere High Pass Size：输入以天为单位的窗口大小。应用于时间分布相关滤波。默认 365
Refinement method：轨道精炼方法。这里选择 Polynomial Refinement
Refinement Residual Phase Polynomial Degree：精炼残余相位多项式次数。在重去平处理中用于估算相位解缠结果中的斜坡相位。当 GCP 的数量小于这个次数要求的数量时，该多项式次数会自动降低。默认为 3
设置好后点击Exac开始第二步反演
9.地理编码
此步为最后一步，将输出的结果根据坐标系进行投影，前面的结果均是失真的图像，在这一步将会校正好。
步骤：打开/SARscape/Interferometric Stacking/SBAS/6 - Geocoding
(1)在Input Files窗口中选择我们之前得到的auxiliary.sml文件
(2)在 DEM/Cartographic System 面板中，选择 之前得到的DEM文件
(3)高度和速度阈值（Height/VelocityPrecision Threshold）建议保持默认，10和8
(4)Rebuild All：False
(5)Make Geocoded Raster：是否创建地理编码栅格，选True，将生成输出栅格文件
Make Geocoded Shape：是否创建地理编码矢量，选True，将生成输出矢量文件
Generate Shape Time Series：生成时间序列，选True，将生成输出矢量文件，其中包含位移时间演变
(6)Displacement Custom Direction：False。设置为 True 时，形变和速率将投影到自定方向上，需要设置方位角（Azimuth Angle：以度为单位，从北顺时针方向）和倾斜角（Inclination Angle：以度为单位，从水平面开始计算）注：得到的形变和速率产品，默认情况是视线方向（LOS）
(7)X Dimension (m)：输出像元 X 方向大小，以米为单位
Y Dimension (m)：输出像元 Y 方向大小，以米为单位
注：投影坐标系为经纬度坐标时，当设置像元值大于 0.2，软件内部会粗糙的自动转换为度的单位，当设置像元值小于 0.2，软件自动识别以度为单位
(8)Mean Window Size：对高度图像结果做均值滤波处理
Interpolation Window Size：对结果图像中的无效值用窗口大小内的像素平均值差值，设置 0 表示不进行插值
其他：
Shape Max Number of Points：最大点数，如果输入了-1以外的值，则输出形状文件将拆分为多个部分（每个部分都用渐进编号标记 - _1; _2;形状的每个部分都包含总点数的一部分，该部分对应于输入的值。建议不要超过200000点，以避免可视化问题。
Generate Shape Time Series：生成时间序列，选True，将生成输出矢量文件，其中包含位移时间演变。
Generate Shape Counter Series：生成矢量计数器序列，通过设置此标志，将生成输出矢量文件，其中包含每个采集日期的有效干涉图计数器。当计数器为零时，对应的采集日期已被时间插值。计数器越高，相应的估计位移测量的可靠性越高。
Vertical Displacement：垂直形变，形变和速率将投影到垂直方向，如果实际位移方向未知，则不要使用重新投影。
Slope Displacement：斜面形变，形变和速率将投影到斜距方向，如果实际位移方向未知，则不要使用重新投影。
Displacement Custom Direction：自定义形变方向，设置为 True 时，形变和速率将投影到自定方向上，需要设置方位角（Azimuth Angle：以度为单位，从北顺时针方向）和倾斜角（Inclination Angle：以度为单位，从水平面开始计算）。
Azimuth Angle：方位角，自定义方向 方位角 从北向以度为单位测量 - 顺时针方向。
Inclination Angle：入射角，自定义方向倾角，以与水平平面的度数为单位进行测量
点击Exac即开始地理编码
编码结果：
_term0_geo：反演模型 0 次多项式结果（单位：毫米/年）。没有选择“No DisplacementModel”方法情况下生成。
_term1_geo：反演模型 1 次多项式结果（单位：毫米/年）。没有选择“No DisplacementModel”方法情况下生成。
_term2_geo：反演模型2次多项式结果（单位：毫米/年2）。没有选择“No DisplacementModel”方法情况下生成。
_term3_geo：反演模型3次多项式结果（单位：毫米/年3）。没有选择“No DisplacementModel”方法情况下生成。
_vel_geo：平均形变速率（单位：毫米/年）。没有选择“No Displacement Model”方法情况下生成。
_acc_geo：平均形变加速度(单位：毫米/年 2 )。没有选择“No Displacement Model”方法情况下生成。
_delta_acc_geo：平均形变加速度变化(单位：毫米/年 3 ) 。没有选择“No DisplacementModel”方法情况下生成。
_precision_velocity_geo：估算平均形变速率的精度（单位：毫米/年）。
_disp_geo： 相对一个日期的形变量（单位：毫米）。
_disp_frst_geo：大气校正前，相对一个日期的形变量（单位：毫米）。
_height_geo：用于调整输入 DEM 的高度值（单位：米）。
_precision_height_geo：估算 DEM 调整高度值的平均精度。
_dem：高度调整之后的 DEM。
产品结果主要索引文件包括：
_geo_vel+height_meta： 索引文件，地理坐标下的高度和形变速率测量值。
_geo_otherinfo_meta：索引文件，地理坐标下强度均值、多期相干系数、高度量测精度和高度修正值。
_geo_disp_first_meta：索引文件，大气校正前，在地理坐标下每期形变量。
_geo_disp_meta：索引文件，大气校正后，在地理坐标下每期形变量。
本文参照以下文章进行实操：
Esri公司的《ENVI SARscape入门教程》
CSDN博主「Vigo*GIS_RS」的文章https://blog.csdn.net/qq_41159191/article/details/127473610
