西藏定日县6.8级地震InSAR处理详细操作教程

据中国地震台网中心测定：北京时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：地震形变结果制图表达