SARscape形变建模工具的使用

本文以云南漾濞地震形变建模为例，介绍该工具的使用操作。

软件版本ENVI5.6+SARscape5.6

1. 生成采样区

生成图像采样区，定义采样区域以及采样点之间的间距，具体操作如下：

（1）打开/SARscape/Interferometry/Displacement Modeling/Sampling Areas，在Sampling Areas面板上，点击Select raster image按钮，选择DInSAR得到的形变结果文件_disp；

（2）在界面中自动显示形变结果文件，在形变结果上，绘制两个多边形，鼠标操作功能如下：

鼠标左键绘制多边形，左键双击结束绘制
在绘制过程中，点击右键取消绘制的最后一个节点
鼠标滚轮可以放大/缩小/平移图像

（3）在右侧的Sampling Areas中，Resolution字段设置对应多边形的采样间隔，外围的区域采样间距2000米，内部的形变区域采样间隔500米。；

（4）点击Save to shapefile，将采样区域保存为shp文件SamplingAreas.shp。

生成采样区

2. 图像采样

图像的子采样可以通过两种方法进行：在特定区域使用规则的点网格或使用四叉树算法。本例中，使用的是第一种方法，即：使用上一步创建的SHP文件。具体操作如下：

（1）打开/SARscape/Interferometry/Displacement Modeling/Image Subsampling工具，打开Image Subsampling面板；

（2）根据默认的扩展名提示，输入DInSAR得到的各个结果数据：

Subsampling image：DInSAR得到的形变结果_disp
Azimuth LOS image：角度数据_ALOS
Incidence LOS image：角度数据_ILOS
DEM：初始参考DEM
Subsampling method：Mesh from vector file

图像的子采样可以通过两种不同的方法进行：在特定区域使用规则的点网格或使用四叉树算法。在本例中，使用第一种方法，Mesh from vector file，即使用SamplingAreas矢量文件。

Sampling Areas：上一步生成的采样区文件
Output Shapefile：设置输出的采样点矢量文件

注：本例中采用的数据坐标系为WGS84。

（3）点击Start，进行图像采样计算，运行完成之后，提示如下对话框，点击“是”，打开采样点可视化界面。

注：对于重叠区域(内多边形完全包含在外多边形中)，采样算法选取分辨率更小的点，在本例中为500 m。

提示：一般情况下，任意数量的数据集可以同时被反转:为每个数据集保持合理数量的点，以避免计算上的额外负载。

提示：总是检查shapefile属性表(使用ENVI或任何GIS软件)以检查内容(参见“InSAR_dataset”联机帮助了解更多细节)

3. Non-Linear Inversion非线性反演

非线性反演用于寻找最佳拟合的源参数，即那些更好地预测观测数据的参数。该建模基于一个XML项目文件进行。

在非线性反演中，利用地球物理震源计算观测位移。我们假设对形变源未知，所有参数都从DInSAR结果中推断出来，由于干涉图中的参考点可能设置在一个变形区域，因此我们通过一个单断层对数据进行建模，也考虑了可能影响数据的偏移量。

（1）点击/SARscape/Interferometry/Displacement Modeling/Non-Linear Inversion工具，打开Non-Linear Inversion面板；

（2）点击New，设置工程文件路径和文件名，新建一个工程文件Yangbi_project.xml。

（3）数据集设置：

点击Add Datasets->add from file，选择第2步生成的采样点文件INTERF_out_reflat_sampled.shp文件，提示该文件的坐标系信息，点击“确定”，输入文件的坐标系定义了该项目的坐标系，如要修改输出坐标系，在源设置面板中修改。

（4）数据集参数设置

对于InSAR数据集，默认选项是评估线性斜坡的参数。在本例中，我们只评估一个持续的变化。双击数据集名称或选中该数据集点击EDIT DATASETS，在数据集参数设置面板，Invert for an orbital surface选项选择Constant offset项。

（5）源设置

设置源初始参数进行反演，有几种方法，本例中使用CMT网站查询到的地震信息作为源初始参数，访问网站：https://www.globalcmt.org/CMTsearch.html，在网站上，设置时间2021年5月21日，搜索得到本次云南漾濞地震的信息为：202105211348A。

在INPUT Sources界面，点击ADD SOURCES->Import->From Global-GMT catalog…，在打开的面板上，输入Insert CMT id：202105211348A，点击Search。

CMT source有两个，根据本次地震情况，选Plane2，表格中列出了所有源参数的平均值、最小值和最大值。点击Commit，源被添加到INPUT Sources列表中。

注：虽然CMT解包含点源参数(走向、倾角和倾角、震源位置和深度、地震矩)，但长度、宽度和滑动均可根据Wells & Coppersmith(1994)规则自动推导。

注：可以在这里加入多个数据集，如升降轨的形变结果，GPS数据集，GPS数据集支持离线创建并以矢量格式存储的文件。

（6）经纬度范围设置

经检查，源参数面板默认的经纬度范围，由于地震震源机制的位置相对于实际地震位置，有所偏差，建议在这一步根据形变图来确定反演的经纬度范围。

在INPUT Sources面板，双击源或点击EDIT SOURCES – Edit Parameters…，打开Elastic dislocation（Okada）面板。点击Pick from map按钮，在Select a map面板，选择图像采样文件_sampled.shp，在打开的采样点可视化图像上，鼠标左键绘制形变区范围，LON和LAT的范围实时根据绘制区域调整，确定一个区域即可。点击Close关闭可视化面板。

在Elastic dislocation（Okada）面板，点击Commit提交初始源参数。

（7）进行非线性反演

数据集和源设置好之后，点击Save，保存xml工程文件，点击Start，开始反演。反演完成后，提示对话框运行震源机制解，在此选“否”，后面会进行这一部分操作。

（8）反演结果的查看

处理结束后，四个新的选项卡被添加到面板中：“OUTPUT dataset”，“OUTPUT Sources”，“OUTPUT Statistics”和“Report”。在“OUTPUT dataset”中，选择_nonlinear。然后点击“VIEW DATASETS - Plot 2D”查看DInSAR形变结果、建模形变结果与残差之间的比较。

选择“Report”选项卡查看结果的综合和定量概述。

在结果概述中，检查所有源参数是否都在给定的最小/最大范围内，“***”字符串标识的是该项达到了上限或下限，出现这种情况时，需要扩展该参数的范围（Range）以使该参数达到其最佳拟合值。

（9）源参数调整

本例中，虽然默认源参数得到的反演结果图看起来不错。但是宽度（Width）、深度（Depth）和倾角（Dip）达到了[min,max]范围的边界。显示为***，这种情况下，需要调整范围并重新运行非线性反演。由于倾角不能大于90，而倾角超过90意味着倾角从NE转向SW，因此需要将走向改变180°。将以下参数进行调整：

Strike Range [270, 360] ---> [90, 180] (走向改变180°，此时倾角从NW转向SW)
Dip interval [66, 90] (保持不变，因为走向转到SW之后，断层依然近乎垂直)
Rake Range [123, 213] --> [-213,-123] (当倾角方向改变时，上盘和下盘就互换了)
Depth [1236.4, 10000] ---> [0, 10000] (让源更浅)

在INPUT Sources面板，双击源或点击EDIT SOURCES – Edit Parameters…，打开Elastic dislocation（Okada）面板，修改相应参数。点击Commit。

（10）再次进行非线性反演

在Non-Linear Inversion面板，点击Save，保存当前工程，点击Start，再次进行非线性反演。查看结果数据集，结果查看参考（8）。

（11）结果检查与优化，本次优化之后，得到了更好的结果，发现还有一项深度（Depth）范围需要调整。重复（9）、（10）步骤，继续调整和运行反演，直到得到没有***的指标。

最终本次Okada反演参数和反演结果如下：

注：负的深度值表示高于海平面的顶部断层。当“Compensate for topography”（地形补偿）选项设置为Yes（默认值），可以设置负的深度值。

注：1、2步的采样多边形创建的不同，会得到稍有区别的结果。

在找到最优拟合解后，我们研究了参数不确定性和它们之间可能的权衡的存在性。在“OUTPUT Sources”选项卡，点击“EDIT Sources - Calculate Statistics”，等待运行结束。通过每次添加一个空间相关的随机噪声(运行次数在“选项”选项卡中)，可以执行50次反演。最后可以在“OUTPUT Statistics”选项卡中看到结果，使用“VIEW SOURCES - Plot parameters…”或者使用“VIEW SOURCES - Plot 3D”图形工具

注：在所有图形查看器(2D, 3D，对于数据集和源)中，可以通过鼠标左键功能来调整窗口大小、选择和移动/缩放/旋转单个元素、更改元素属性(双击以访问属性面板)，可以将图保存为TIFF、PDF或JPG格式。

非线性反演的输出文件有：

_nonlinear.shp——反演得到的源矢量文件

_sampled_nonlinear.shp——反演得到的形变结果网格点矢量文件

_statistics.shp——50次反演的源矢量统计结果

非线性反演结果

4. Linear Inversion线性反演

对断层平面上的滑动分布进行线性反演。在本教程中，我们将检索一个带有固定耙子的分布式正滑，但也提供了一套完整的替代方案(在给定的时间间隔内可变耙子，完全自由耙子，等等)。

（1）点击/SARscape/Interferometry/Displacement Modeling/Linear Inversion工具，打开Linear Inversion面板。

（2）点击Open，选择在上一步非线性反演中生成的项目文件Yangbi_project.xml，软件检查是否存在完整的非线性反演，点击“是”，提示允许自动设置输入数据集和源，点击“确定”。

（3）数据集设置：

在Input Datasets面板，已经自动输入INTERF_out_reflat_sampled.shp采样数据，双击数据集打开设置面板，Invert for an orbital surface选项选择Constant offset项，评估一个恒定的轨道面。

（4）源设置：

在INPUT Sources面板，双击源或点击EDIT SOURCES – Edit Parameters…，打开Elastic dislocation（Okada）面板。该面板的参数是根据非线性反演的输出参数进行设置的。需要注意的是，这个源代表的是具有均匀滑移的平均解，必须在宽度和长度上加以扩展，以完全包含滑移分布。如下图所示：

然后，在走向和倾向方向上划分为若干个小格网，参数面板如下所示：

默认使用了“Distributed and positive slip with fixed rake”约束，下拉菜单中提供了更多的方案。

（5）进行非线性反演

（6）反演结果的查看

处理结束后，三个新的选项卡被添加到面板中：“OUTPUT dataset”，“OUTPUT Sources”，和“Report”。点击“Report”选项卡查看结果的综合和定量概述。

点击OUTPUT Sources面板的VIEW SOURCES->Plot 2D – Frontal View（fault only），查看断层滑动分布的二维图，发现过于平滑，如下图所示：

线性反演断层滑动分布的效果由“Damping Factor”（阻尼因子）控制，在“Options”面板，将默认的Damping Factor 0.1降低到0.02，并再次开始反演，以获得更可靠的结果。如下图所示：

注：1、阻尼因子是一个试错参数，可不断尝试，直到得到一个满意的结果。

2、阻尼因子的设置，文献中提出了几种自动定义的方法，但它是一个经验参数，其设置主要与用户经验和判断有关。

在“OUTPUT dataset”中，选择_sampled_linear。然后点击“VIEW DATASETS - Plot 2D”查看DInSAR形变结果、建模形变结果与残差之间的比较。

线性反演的输出文件有:

_linear.shp——具有模型位移点的矢量文件

_fault_linear——基于模型滑移分布的多边形网格矢量文件

5. 查看震源机制解

在包含地震源的非线性和线性反演之后，会自动调用震源机制解面板。如果要手动打开查看震源机制解，可以打开工具/SARscape/Interferometry/Displacement Modeling/Modeling Tools/Calculate and Draw Focal Mechanism，在面板的INPUT Sources选项卡中，点击ADD SOURCES->Import->From XML project file…，选择之前反演时生成的工程文件Yangbi_project.xml。

在下图的XML文件管理器中，Section下拉菜单中选择LINEAR INVERSION OUTPUT，然后选择2021_05_21-YUNNAN,CHINA项，点击Add。

在列表中，选择打开的源，点击Start，可以显示震源机制解。

注：可以点击Save to PDF…按钮，将结果保存为pdf文件。

6. CFF应力转移

CFF应力转移（CFF Stress Transfer）工具可用于计算从一列输入源到一列接收源所引起的应力变化。在漾濞地震的情况下，只涉及一个震源，因此我们在此计算其断层本身引起的应力变化。自感应力变化可以用来验证断层面上的余震分布是否与预期一致，即是否位于应力变化最大的地方。

（1）点击/SARscape/Interferometry/Displacement Modeling/CFF Stress Transfer工具，打开CFF Stress Transfer面板；

（2）点击Open，选择在之前反演时生成的项目文件Yangbi_project.xml（在此也可以新建项目文件）；

（3）在INPUT Sources面板，点击ADD SOURCES->Import->From XML project file…，选择线性反演时保存的工程文件Yangbi_project.xml；

（4）在打开的XML文件管理器中，Section下拉菜单中选择LINEAR INVERSION OUTPUT，然后选择2021_05_21-YUNNAN,CHINA项，点击Add；

（5）切换到Receiver Sources面板，点击ADD SOURCES->Import->From XML project file…，选择线性反演时保存的工程文件Yangbi_project.xml；

（6）在打开的XML文件管理器中，Section下拉菜单中选择LINEAR INVERSION OUTPUT，然后选择2021_05_21-YUNNAN,CHINA项，点击Add；

（7）点击Save，保存工程文件；

（8）点击Start，进行CFF应力计算。

计算完成之后，加载了OUTPUT Sources面板，在该面板上，VIEW SOURCES->Plot 2D Frontal View（fault only），可以查看二维应力的结果。

Options界面下提供了标准模式（Standard Formulation）和不排水模式（Undrained Formulation）的选项，Mask singularities at the input source proximity选项能将输入源和接受源距离太小的位置进行掩膜。如下两图是该参数勾选与不勾选的结果。

CFF应力转移计算的输出结果有：

输出CFF值存储在XML工程文件中，在“ModelingRoot”-“CFFModeling”-“CFFOutput”-“CFFTarget”部分

2021_05_21-YUNNAN, CHINA_cff.shp——矢量文件

7. Forward Modeling正演模型

该功能可以由一个或多个源预测地表位移。输出可以是栅格或矢量点文件。在本教程中，我们使用线性反演后定义的源来预测表面位移，输出栅格文件。

（1）点击/SARscape/Interferometry/Displacement Modeling/Forward Modeling工具，打开Forward Modeling面板；

（2）点击Open，选择在之前反演时生成的项目文件Yangbi_project.xml（在此也可以新建项目文件）；

（3）在INPUT Sources面板，点击ADD SOURCES->Import->From XML project file…，选择线性反演时保存的工程文件Yangbi_project.xml；

（4）在打开的XML文件管理器中，Section下拉菜单中选择LINEAR INVERSION OUTPUT，然后选择2021_05_21-YUNNAN,CHINA项，点击Add；

（5）输出定义：在Options界面，Forward Model Output默认为Raster，点击Set raster info按钮，在Raster Parameters界面上，设置COORDINATE SYSTEM为Geo WGS84，点击Get from map按钮，在底图上绘制预测形变的区域，绘制好之后点击Close，区域的范围就自动填写到面板上了，设置CellSize为0.002，点击Commit。

（6）选择DEM文件，Output File Root Name设置为2021_05_21-YUNNAN_disppl

（7）点击Save，点击Start，进行正演计算。计算完成之后，面板中不会创建额外的选项卡，输出的结果包含三个栅格数据，分别是三个方向的位移分量：垂直、北、东方向：

2021_05_21-YUNNAN_disppl_up

2021_05_21-YUNNAN_disppl_north

2021_05_21-YUNNAN_disppl_east

如下图是位移的北分量，在ENVI中彩色渲染的效果，从断层的横断面图来看，断层量测发生了向北和向南的滑移。

8. 投影到LOS方向

将位移的东、北和垂直分量投影到卫星视线(LOS)中。投影过程需要传感器方位角和观察角度，这些在_ALOS和_ILOS文件中，是DInSAR的输出结果。这个面板还允许添加一个轨道表面。这有助于更好地模拟观测到的信号。在本教程中，在非线性和线性反演中评估了常量偏移量。要检索轨道表面方程(在本例中是一个简单的偏移量)，使用XML Project File Explorer(带有“Get from XML…”按钮)，并导航到下拉菜单中的“LINEAR INVERSION OUTPUT”部分，检索评估的轨道贡献:

（1）打开工具/SARscape/Interferometry/Displacement Modeling/Modeling Tools/Project raster to LOS，在面板上依次输入_east、_north、_up、_ALOS、_ILOS数据，其中，前三个来自于正演计算，后两个来自于DInSAR。

（2）在ADVANCED OPTIONS选项中，点击Set from XML，输入前面生成的工程文件Yangbi_project.xml，在XML Project File Explorer面板，下拉菜单选择“LINEAR INVERSION OUTPUT”，选中线性反演得到的轨道分量，点击Add。