据中国地震台网中心测定：北京时2025年1月7日9时5分，西藏日喀则市定日县（北纬28.5度，东经87.45度）发生6.8级地震，震源深度10千米。
本文以哨兵1A作为数据源，使用DInSAR的方法对本次地震进行干涉测量处理。本文旨在介绍软件处理操作，结果仅供参考，准确结果以官方发布为准。
数据情况如下表所示：
|
卫星 |
Sentinel1A |
|
获取模式 |
IW |
|
数据级别 |
SLC |
|
极化方式 |
VV |
|
轨道方向 |
升轨 |
|
地面分辨率 |
5*20米 |
|
成像时间 |
2025年1月1日（震前）、2025年1月13日（震后） |
|
数据轨道 |
第一景Path：121，Frame：499 第二景Path：121，Frame：494 注：地震区域刚好位于两景影像中间位置 |
|
参考DEM数据 |
ALOS World 3D DEM |
本例子大约需要40g硬盘空间。
本文使用的软件为SARscape6.1，试用申请地址：https://envi.geoscene.cn/sar_license
哨兵1A数据下载网址：https://search.asf.alaska.edu/
1、数据导入
处理数据之前，设置哨兵1数据的系统参数，打开/SARscape/Preferences/Preferences specific面板，选择Load Preferences->Sentinel TOPSAR(IW-EW)。
打开数据导入工具/ SARscape/Import Data/SAR Spaceborne/Single Sensor/Sentinel-1。在打开的面板中，
- 数据输入面板（Input Files）
输入下载好的四个zip文件
- 可选文件面板（Optional Files）
输入所需处理地震范围的地理范围矢量文件
- 参数设置面板（Parameters）：主要参数（Principal Parameters）
- 生成强度数据（Make power QL）：True，数据导入时生成强度数据，以便查看；
- 极化方式（Polarization）：VV Only，只输出同极化数据（DInSAR处理首选同极化数据）；
- 同轨道镶嵌（Make mosaic same track）：True，本实例中需要两景同轨道数据，输入同轨道相邻数据自动进行镶嵌；
- 不进行子区域裁剪（Skip Sample Selection）：False，使用输入的地理范围进行子区域裁剪；
- 对数据重命名（Rename the File Using Parameters）：True。软件会根据输入数据的文件名进行输出文件的自动命名。
其他参数默认。
- 数据输出面板（Output Files）：输出文件（Output file list）：自动读取ENVI默认的数据输出目录以及输入面板中的数据文件名。
注：1.如果要修改输出的路径，在右边单击文件夹图标选择输出文件夹目录。
2、如果要修改输出的路径，在输出文件名右键选择Edit菜单。
图1：数据输入面板参数设置
输出的文件大小约为11g，包括：整景图像的强度图数据（_pwr）、slc索引文件（.slc_list）带有地理坐标的外边框矢量文件（.shp）。
注：本次数据处理时，精密轨道文件（数据获取时间之后21天发布）还未发布，故本次数据导入没有使用精密轨道文件（导入过程中会提示找不到精密轨道文件，忽视即可），轨道误差仅在轨道精炼中使用GCP拟合多项式去除。一般情况下，做InSAR处理，导入哨兵数据建议使用精密轨道文件。
2、基线估算（选做）
打开工具/SARscape/Interferometry/Interferometric Tools/Baseline Estimation，输入主从影像的_slc_list文件，点击Exec按钮执行，输出基线估算的结果。
图2：基线估算面板
基线估算结果：
|
Normal Baseline (m) = -15.211667 Critical Normal Baseline min - max(m) = [-5742.915729] - [5742.915729] Range Shift (pixels) = 1.218750 Azimuth Shift (pixels) = 0.395508 Absolute Time Baseline (Days) = 12 Doppler Centroid diff. (Hz) = 1.091782 Critical min-max (Hz) = [-486.486310] - [486.486310] 2 PI Ambiguity height (InSAR) (m) = 934.040367 2 PI Ambiguity displacement (DInSAR) (m) = 0.027733 1 Pixel Shift Ambiguity height (Stereo Radargrammetry) (m) = 78459.393555 1 Pixel Shift Ambiguity displacement (Amplitude Tracking) (m) = 2.329562 Slant Range Distance reference (m) = 851784.651592 Slant Range Distance seconsary (m) = 851787.489639 Reference Incidence Angle = 36.975620 Absolute Incidence Angle difference = 0.001014 Reference Azimuth Angle = -79.656940 Absolute Azimuth LOS Angle difference = 0.000228 *Potentially suitable for: Standard Interferometry, Persistent Scatterers Stacking, Distributed Scatterers Stacking, Amplitude Tracking. This pair can also be used for: Amplitude Ratio. |
基线估算的结果显示，这两景数据的空间基线为15.2米，远小于临界基线5742.9米，时间基线12天，做DInSAR一个相位变化周期代表的地形变化为0.028米。
3、DInSAR处理
打开/SARscape/Interferometry/DInSAR Displacement Workflow工具，使用DInSAR流程化工具进行差分干涉处理，步骤包括：数据对输入、干涉图生成（包括数据对配准、干涉相位计算、去平）、干涉图滤波、相位解缠、轨道精炼和重去平、相位转形变。
在进行DInSAR工作流之前，需要准备该区域的参考DEM，本文使用的是ALOS World 3D DEM。可使用/SARscape/General Tools/DEM Extraction/ALOS World 3D 30M工具自动下载。
一、 数据输入
DInSAR Displacement Workflow工具左侧第一步Input File面板，输入如下数据：
- Input File面板，输入20250101的VV极化方式的_slc_list作为参考影像，20250113的VV极化方式的_slc_list作为第二影像；
- DEM/Cartographic System面板，输入参考DEM文件；
- Parameters面板，Grid Size：15
图3：输入干涉像对的slc数据
输入的数据文件设置好之后，点击Next。由于之前已经做了数据导入，故Import Generic SAR Date这步略过，直接点击Next到Interferogram Generation界面。
二、 干涉图生成
这一步对干涉像对进行配准，并计算干涉相位，生成干涉图，使用参考DEM对已知地形相位去除，参数说明如下：
- 主要参数（Principal Parameters）：按照默认。
- 全局参数（Global）
生成TIFF数据（Generate Quick Look）：Ture，生成结果图的TIFF格式快视图，方便查看中间结果，其他步骤类似。
处理完成之后，ENVI视窗中自动加载了去平后的干涉图，以及主从影像的强度图，这一步生成的数据结果存放在ENVI默认的临时文件路径下，默认路径为：C:\Users\<用户名>\AppData\Local\Temp，自动生成一个名为SARsTmpDir_DDMMYYYY的文件夹（有些系操作系统版本可能不生成，文件直接在根目录中），这一步生成的结果有：
|
结果命名 |
内容 |
|
INTERF_out_master_pwr/ INTERF_out_slave_pwr |
参考影像/第二影像的多视强度图 |
|
INTERF_out_par |
干涉像对配准时生成的偏移量数据 |
|
INTERF_out_par_orbit_off |
配准时估算轨道偏移时用到的点矢量数据，包括在方位向和距离向 的点的位置坐标、测量到的偏移量、计算出的线性偏移量。 |
|
INTERF_out_par_winCC_off |
配准时从强度数据的相差上估算交叉相关偏移量的点矢量数据。包含每个点上交叉相关的值（CC），范围是0-1。 |
|
INTERF_out_par_winCoh_off |
配准时从相位数据的相差上估算相干性的点矢量数据，包含信噪比（SNR）和相干性，相干性值的范围是0-1 |
|
INTERF_out_srdem |
地距转为斜距的DEM数据 |
|
INTERF_out_sint |
地形合成相位 |
|
INTERF_out_int |
干涉相位 |
|
INTERF_out_dint |
去平的差分干涉相位 |
图4： 去已知地形后（去平）差分干涉图_dint
三、 干涉图滤波和相干性计算
对上一步生成的差分干涉图进行滤波，对干涉图的噪声进行一定程度的抑制，计算干涉像对的相干系数。干涉图滤波有4种方法（Adaptive、Boxcar、Goldstein、Adaptive Non Local InSAR）可选，本次使用了Goldstein滤波方法，这种滤波方法的滤波器是可变的，提高了干涉条纹的清晰度、减少了由空间基线或时间基线引起的失相干的噪声，这种方法是最常用的方法。
- 主要参数（Principal Parameters）：
- 滤波方法（Filtering Method）：Goldstein
如果滤波效果不好，可以修改以下三个参数，本例按照默认不需要修改。
选择滤波参数（Filtering）：
- Goldstein Win Size：默认为32。该值越大，滤波器对小细节(如局部干涉条纹)的灵敏度就越低，该值最好使用2的倍数，该值范围为：32~512。
- Goldstein Min Alpha：默认3，值越大滤波器平滑越强，该值范围为：0.3~3。
- Goldstein Max Alpha：默认5，值越大滤波器平滑越强，该值范围为：0.5~4。
- 全局参数（Global）
- 生成TIFF数据（Generate Quick Look）：Ture，生成结果图的TIFF格式快视图，方便查看中间结果，其他步骤类似。
点击Next按钮，进行干涉图滤波和相干系数计算，处理完成之后，自动加载滤波后的干涉图_fint和相干性系数图_cc。这一步处理之后生成的结果有：
|
结果命名 |
内容 |
|
INTERF_out_fint |
自适应滤波后的差分干涉图 |
|
INTERF_out_cc |
相干系数图 |
图5： 滤波后的差分干涉图
图6：相干系数图
相干系数图是0-1的灰度图像，值越大表示相干性越高，地表发生变化的区域相干性低，从相干系数图上可以明显的看到地震区域地表了发生了变化，该区域相干系数低。
四、 相位解缠
干涉图中的相位，从-π到π呈周期性变化，当真实相位值大于π时，相位会重新从-π开始，以2π为周期循环。相位解缠是将相位由差分相位恢复为真实值的过程。
- 主要参数（Principal Parameters）
- 解缠方法（Unwrapping Method Type）：Minimum Cost Flow
- 解缠分解等级（Unwrapping Decomposition Level）：1
- 解缠最小相干系数阈值（Unwrapping Coherence Threshold）：3。相干系数大于该阈值的像元进行解缠。
点击Next按钮，进行干涉图滤波和相干性生成处理，处理完成之后，自动加载相位解缠结果图_upha。这一步处理之后生成的结果有：
|
结果命名 |
内容 |
|
INTERF_out_upha |
相位解缠结果 |
图8：相位解缠结果
五、 轨道精炼
使用一组代表稳定位置的GPC点，进行多项式拟合，对轨道误差进行去除，优化相位解缠的结果。
（1）GCP选择
在GCP Slection面板，点击 按钮，自动打开流程化的控制点选择工具，软件自动在数据输入位置输入了相应的文件（Input file为_fint，reference file为_upha，DEM file为参考DEM文件）。点击Next按钮，在控制点选择面板，鼠标变为选点的十字丝状态，单击鼠标左键就可以选点。在远离形变条纹、相干性高、解缠结果平滑的位置，选择若干控制点。
- GCPs：GCP列表，每个GCP点的坐标在该面板上显示
- Cartographic System：默认，与参考DEM坐标系一致
- Export：输出GCP文件，注意路径不要输出到中文文件夹。如果默认输出到有中文名文件夹换个路径，否则保存的文件为空。
图9： 选择控制点流程化工具
在控制点生成面板上，点击Finish按钮，生成的控制点文件.xml，自动添加到面板中。点击Next按钮。
（2）轨道精炼
使用多项式模型，结合GCP点的相位，解算多项式系数，拟合多项式，对干涉图和解缠后的相位进行轨道精炼和重去平，得到优化的相位结果。
- 轨道精炼半径（Refinement Radius（m））：默认5。
- 多项式系数（Refinement Res Phase Poly Degree）：3，默认值3表示在距离和方位角方向上的相位残差与恒定的相位偏移将被校正。
- 使用DEM配准（Coregistration With DEM）：True
点击Next按钮，运行得到轨道精炼后的干涉图_reflat_fint和解缠相位图_reflat_upha，这一步处理之后生成的结果有：
|
结果命名 |
内容 |
|
INTERF_out_reflat_fint |
优化后的滤波差分干涉图 |
|
INTERF_out_reflat_sint |
优化后的合成相位图 |
|
INTERF_out_reflat_upha |
优化后的解缠相位 |
|
INTERF_out_reflat_refinement.shp |
轨道精炼用到的有效的控制点文件，斜距坐标系 |
|
INTERF_out_reflat_refinement_geo.shp |
轨道精炼用到的有效控制点文件，地理坐标系 |
|
INTERF_out_reflat.txt |
轨道精炼结果报告 |
同时输出了使用这组GCP点拟合的轨道精炼多项式及误差信息：
|
ESTIMATE A RESIDUAL RAMP Valid points found = 13 Extra constrains = 2 Polynomial Degree choose = 3 Polynomial Type : = k0 + k1*rg + k2*az -4.9548224530 0.0002847385 0.0003529076 Polynomial Coefficients (radians) : k0 = -4.9548224530 k1 = 0.0002847385 k2 = 0.0003529076 Root Mean Square error (m) = 145.3608234122 Mean difference after Remove Residual refinement (rad) = 0.0853880306 Standard Deviation after Remove Residual refinement (rad) = 0.9743980092 |
可以看出，本次使用了13个GCP点，全部位于有效像元，使用这13个GCP点，拟合的轨道精炼的多项式为：-4.9548224530+ 0.0002847385*rg -0.0003529076*az，均方根误差RMSE为145.36米，去除残差后的平均差值为0.08538弧度，去除残差后的标准差为0.974398弧度。
在面板上点击Next按钮，进入下一步。
图10： 轨道精炼和重去平前后的差分干涉图（左：原始结果，右：轨道精炼后）
|
结果命名 |
内容 |
|
INTERF_out_reflat_fint |
优化后的滤波差分干涉图 |
|
INTERF_out_reflat_sint |
优化后的合成相位图 |
|
INTERF_out_reflat_upha |
优化后的解缠相位 |
|
INTERF_out_reflat_refinement.shp |
轨道精炼用到的有效的控制点文件，斜距坐标系 |
|
INTERF_out_reflat_refinement_geo.shp |
轨道精炼用到的有效控制点文件，地理坐标系 |
|
INTERF_out_reflat.txt |
轨道精炼结果报告 |
六、 相位转形变与地理编码
将经过轨道精炼和重去平的解缠相位，转换为形变，并地理编码到地理坐标系。
- 主要参数（Principal Parameters）
- 相干性阈值（Product Coherence Threshold）：0.3。相干性大于该值的相位转为形变值。
- 垂直方向的形变（Vertical Displacement）：False。不计算垂直方向上的形变。
- 斜坡形变（Slope Displacement）：False。不计算斜坡方向上的形变。
- 用户自定义方向的形变（Displacement Custom Direction）：False。
- 方位角（Azimuth Angle）：0。自定义方向的时候设置该角度。
- 入射角（Inclination Angle）：0。自定义方向的时候设置该角度。
- X方向上的水平分辨率（X Dimension（m）：15。X方向制图分辨率。
- Y方向上的水平分辨率（Y Dimension（m）：15。Y方向制图分辨率。
- 地理编码参数（Geocoding）：设置Dummy Removal=ture。可以去除外边框无效像素（如果DEM范围大，可能会出现此类区域）
点击Next按钮，进行相位转形变和地理编码处理。地理编码的坐标系与参考DEM一致。默认得到的是LOS方向的形变结果。这一步得到的结果有：
|
结果命名 |
内容 |
|
_dispwf_disp_precision |
数据质量的估算结果图，代表相位转形变的精度 |
|
_dispwf_disp_ILOS |
视线入射角，视线和水平面垂线的夹角 |
|
_dispwf_disp_dem |
重采样到制图输出分辨率上的参考DEM数据，范围与形变结果一致 |
|
_dispwf_disp_cc_geo |
地理编码的相干性系数图 |
|
_dispwf_disp_ALOS |
视线方位角，正值是正北的顺时针方向，负值是正北的逆时针方向 |
|
_dispwf_disp |
形变结果，单位为m，默认是LOS方向上的形变 |
如果提示以下错误，说明默认输出路径中有中文文件夹。不关闭DinSAR流程化工具面板，修改ENVI中的参数File->Preference->Directories->Output Directory，在DinSAR流程化工具面板中点击Next按钮继续处理。
图11：错误提示面板
七、 结果输出
在Output面板，选择是否删除中间结果，点击Finish，自动对形变结果进行分级设色，输出密度分割图。
形变结果是一个单波段的灰度图像，像元值就是形变量，单位为米，正值代表朝传感器方向运动（抬升），负值代表远离传感器方向移动（沉降）。
图12：相位转形变的结果
可利用/SARscape/Basic/Intensity Processing/Geocoding/Geocoding and Radiometric Calibration工具对INTERF_out_reflat_fint干涉图进行地理编码。
4、结果渲染与表达
形变结果是一个单波段的灰度图像，像元值就是形变量。对形变结果进行彩色渲染有两种方式：颜色分级显示、色带显示。
- 颜色分级显示
使用ENVI的密度分割功能实现，在形变结果图上点击右键，点击New Raster Color Size，在密度分割面板，可根据制图需要设置各级DN值范围与颜色，
图13：密度分割分级显示
- 色带显示
默认形变结果是灰度黑白显示，可在形变结果图上点击右键，选择Change Color Table，点击More，点击
，可选择多种配色方案进行彩色渲染。
按钮可以对色带进行颜色调换。
小技巧：使用彩色渲染形变结果，可通过调整形变结果的直方图，凸显形变信息。点击
打开直方图，默认为2%线性拉伸，可手动调整直方图到最大最小位置，稳定位置为处在色度带中部的相近颜色，沉降和抬升在色带的末端，可凸显形变信息。
- 叠加底图
可叠加雷达数据强度图或光学影像作为底图，形变图设置一定透明度显示。本次处理由于形变区域位于山区，可叠加山体阴影图作为底图，操作方法如下：
（1）打开/Terrain/Topographic Shading Tool，选择DInSAR工作流生成的DEM文件_dispwf_disp_dem（也可选择初始参考DEM）.
（2）在山体阴影生成工具界面，调节拉伸方法、渲染颜色、方位角和高度角等参数得到自己满意的山体阴影效果图。
（3）点击Apply按钮输出结果。
图14：使用DEM数据生成山体阴影图
形变结果图放置在山体阴影图上方，在ENVI工具栏中，使用 
调整形变结果图的透明度，便可看到形变结果图叠加在山体阴影图上的效果。
图15：形变结果渲染叠加山体阴影效果图
也可将地理编码后的差分干涉图叠加到山体阴影图上（干涉条纹的图例来自SARscape的Help，保存为图片，使用ENVI注记加入），效果如下：
图16：差分干涉图叠加山体阴影效果图
ENVI工具栏Annotations下拉菜单，可以根据需要添加制图要素，如公里网、色度带、线段、符号、箭头标识等。
图17：地震形变结果制图表达
标签：相位,结果,干涉,InSAR,INTERF,out,形变,定日县,6.8 From： https://www.cnblogs.com/enviidl/p/18672871