ENVISCAPE软件中PS指导手册

1、PS用户手册2014年11月PS处理流程如下图：图1 PS处理流程软件提供了一系列流程化工具来完成PS处理，包括：1、 连接图生成2、 定义感兴趣区域（可选）3、 干涉工作流3.1 配准3.2 干涉生成和去平3.2 平均强度图省工程和振幅离差指数生成（MuSigma）4、第一次反演 4.1相干性、速率和地形残差估算5、第二次反演 5.1 大气相位去除 5.2 相干性、速率和地形残差估算 5.3 形变估算6、地理编码 6.1 速率和精度估算 6.2 形变结果地理编码研究区和数据说明本文档对PS处理进行详细的说明，用户可以参考本文档进行自己的数据处理。文中大多数截图都是来自两组数据，软件版本是SA

2、Rscape5.1+ENVI5.2，有一些特殊的情况，为了提高显示速度，用到了ENVI Classic。研究区是日本莆安市，在离东京不远的千叶市，用到了两组数据集：ENVISAT-ASAR和ALOS-PALSAR，时间空间覆盖相似，便于对比结果，还有GPS测量数据，用来验证PS的结果。ASAR数据集包含34景降轨数据，IS2模式，VV极化，时间从2006年8月26日到2010年10月22日。PALSAR数据集包含23景降轨数据，单极化或双极化（FBS-FBD）模式，HH极化方式，时间从2006年9月4日到2010年8月9日。两组数据是同一个区域的，不同波长、不同传感器、不同入射角，得到的结果很

3、相近，那么说明结果是很接近真实情况的。本文在处理的时候，两个数据集同时处理的，在实际工作中，建议一次性处理一个数据集。图2 ASAR数据集区域（红色）和PALSAR数据集区域（蓝色）两个数据集的时相如下表所示：每个数据都包含三个数据文件：.hdr (ENVI头文件) .sml (SARscape头文件) .kml (Google Earth文件).在本研究区内，有一些地表的沉降或抬升是在预期内的，是油气开采和地下水回填造成的。PS工作流说明开始之前的准备工作在PS工作流开始之前，需要设置ENVI的系统参数和SARscape的系统参数，根据所处理的数据来设置，用一套系统推荐的默认参数。工作流和说

4、明首先，选择既有高的散射又能保持良好干涉的区域（如城市区域）作为研究区是很重要的。出于像素相干性的考虑，对输入数据的个数也是有严格要求的，因为要找到稳定的PS点。如果输入数据不足的话，会对整个区域的相干性估计过高，从而得到估计过高的PS结果，得到错误的形变结果。PS要求至少是20景输入数据。在做干涉的时候，所有数据都匹配到主影像上，利用DEM去平，转换到主影像坐标系。在有精确轨道文件和精确的参考DEM文件，程序是完全自动运行的。如果在卫星轨道或参考DEM的地理位置有误差，需要用几何地面控制点-“Geometry GCP ffile”，对SAR数据以及参考DEM进行校正。在这种情况下，在配准过程

5、中结合GCP计算的偏移，用来对所有输入的从影像的SAR数据校正到主影像上。如果主影像（参考影像）是经过校正的正确的，那么就不需要GCP文件了。干涉步骤完成之后，需要检查这一步的结果，这很重要，尤其是看配准做的如何。PS处理是根据输入的参数：速率和高程区间（velocity and height range），这取决于研究区的实际情况以及我们预期的形变情况。设置参数的关键是要考虑到在两个时相之间，形变敏感度是/2，因此，较小的时间间隔将确定速率的最大范围，以避免速率反演错误。根据处理区域，PS算法采用下面的一种方法： 如果处理的范围大于设置的面积阈值，整个区域会分成若干子区域，每个子区域会选则一

6、个参考点，这样每个区域都会单独估算出参数，大气参数也会不一致，所以最终是以整个镶嵌的结果输出的。 如果处理的区域小于设置的阈值，那么整个区域只用一个参考点。面积阈值是“Area for Single Reference Point (平方公里)”，默认的参数是5*5平方公里，如果处理范围大于5*5平方公里，会采用第一种方法把区域分成若干小块处理，修改这个参数会对所选择的方法有影响，建议就按照默认的5，因为这对大气的估算和去除是比较适合的。PS核心是以下两步：（1）第一次形变反演（inversion: first step） 输出文件：速率、高程、相干性图和差分干涉图，包含大气信息的结果。数据带

7、有后缀名_first； 用线性模型从所有差分干涉图中估算的形变速率和残余高程信息，残差是大气信息（噪声和非线性目标）对相干性高的象元进行计算。（2）第二次反演（inversion: second step ）： 估算差分干涉图（dint）中的大气信息 目标高程的估算和LOS方向上形变速率的估算，在大气补偿估算之前。 这一步会找到更多的PS点，这对估算大气信息是很重要的； 可以对比相干性文件（_cc_first and _cc ），_cc应该比_cc_first相干性高； 如果研究区有很多植被覆盖（避免选这种区域，选择城市区域），可以减小低通滤波的参数（默认是1200，可以设置为500）； 编辑

8、auxiliary文件（参考后面这部分内容），可以重新运行大气估算和生成形变结果。最后一步是PS地理编码，得到地理坐标的矢量和栅格结果。连接图生成这一步是对多时相的SAR数据进行像对连接，组成若干对干涉像对，和SBAS方法不一样，PS的像对连接，无法进行编辑。这一步确定了在PS处理中每一对干涉像对，每个数据都和主影像（参考影像）组成一对像对，所以N个时相的SAR数据，最多会组成N-1对干涉对，连接图生成工具对数据生成可靠的连接，这些像对会在干涉工作流之后，用于PS反演。程序会自动挑选出主影像，在整个的处理中，主影像都是参考影像，其他的所有影像都会配准到主影像上。用户也可以自己选择主影像，不过还

9、是建议让程序自动挑选，用户自己选择主影像有可能会造成连接较少的情况。启动连接图生成工具，在Input Fie面板中，输入所有的ASAR或PALSAR数据，在Option File面板，Input Master File保持空缺，在Output File面板，设置输出文件路径为（Output Root Name）+_PS_processing。图3 ASAR数据处理，主影像不用输入，程序自动选择图4 PALSAR数据处理，主影像不用输入，程序自动选择生成的连接像对都是以设置的基线阈值进行连接的，设置是在PreferencesPersistent ScatterersBaseline Thresh

10、old (%) setting中，默认的阈值是临界基线的50%，超过这个阈值的图像就不会在PS处理中使用。这一步运行完之后，会得到连接图，每个图像都和选出来的主影像组成了像对，在生成的两个图上，一个点代表一个图像，红色的表示丢弃的图像，绿色的表示可用的图像，黄色的表示主影像。Time-Position plot：该图显示的是每个像对的空间基线，x轴是时相，y轴是空间基线。Time-Baseline plot：该图显示的是每个像对的时间基线，x轴是时相，y轴是该像对的空间基线，像对的连线代表时间基线。该图像从时间和空间基线角度直观的展示了数据集的相干性。连接图可以随时查看，用这个工具/SARsc

11、ape/Interferometric Stacking/Stacking Tools/Plot Viewer。 图5 ASAR数据的连接图 图6 PALSAR数据的连接图多视的视数会根据系统设置中的制图分辨率Cartographic Grid Size自动设置，软件能处理spotlight和stripmap干涉模式，ALOS PALSAR数据，FBD和FBS模式中同样极化方式的数据可以一起使用。图7 ASAR数据集中挑选出的主影像（Asar_slc_pwr）图8 PALSAR数据集中挑选出的主影像（Palsar_FBS_slc_pwr）工程文件Auxiliary file工程文件在整个处理中

12、都是会用到的，就像一个控制中心一样，它记录了一些必要的信息，如所有输入数据的文件路径、SBAS或PS中的每一步骤的处理进度，如果中途停止，再次处理时可以接着上次做的像对接着往下做，而不是重头开始。该文件在处理过程中会被更新，在整个处理中，都会以该工程文件作为输入文件，如果用户想重新做某一步，可以在该文件中改相应步骤的处理标识，如：当浏览完干涉结果之后，想调整一个参数再做一次，那么在PS或者SBAS的工程文件中，把相应步骤设置为“NotOK”，标识改为0，当再运行的时候，可以跳过前面已经运行好的步骤，从这一步重新开始做，然后这一步之后的所有步骤也会重新做。如SBAS处理中，想重新做重去平。重去平

13、之后，这一部分的工程文件记录如下：OK 20 0把这部分内容改为下面内容，程序就会重新做冲重去平以及之后的步骤。NotOK 0 0注：图像ID是从0开始依次编号的，可以在连接图生成之后的_SBAS_processingconnection_graph CG_report.txt文件里面查看到编号，像对的命名格式为：_MasterDate_m_ID_SlaveDate_s_ID，这样比较容易查看。如PS处理中，想重新做配准。在配准之后，工程文件中有如下内容：OK 20修改内容如下，程序就会从配准这一步重新开始做。NotOK 0注：PS没有主从影像的标识，因为主影像只有一个，图像的ID编号也是在工

14、程文件的上面记录的。还有一个很重要的文件是work_parameters.sml，在work路径下，这个文件记录了SBAS或者PS处理过程中用到的所有参数。incremental模式，只有SBAS处理中有，目的是节省处理时间，当一个新的影像添加进来时，不会对所有的像对都重新进行处理，而是只对新加进来的这个图像组成的像对进行处理。这种情况需要设置和之前的处理相同的处理参数。如果使用了incremental模式，就不能再执行Connection Graph generation了。开始进行SBAS或者PS处理，在干涉工作流这一步做完之前，都不要改变slc数据集存放路径，这一步完成之后，用户则可以移

15、动slc数据集到其他位置，进行后续的处理。感兴趣区域的裁剪如果已知形变区域，可以先裁剪出一块研究区来，这样可以节省处理时间和空间。用该工具裁剪的话，范围会作用于所有的图像，裁剪的子区域数据集结果会输出到一个新的文件夹下。这一步是可选项，如果不做这一步，则会对整个影像进行PS处理。图9 感兴趣区域裁剪应该定义一个比实际感兴趣区域（至少200像素）要大的区域，因为处理中过滤器可能会导致一些边界效应。在这些过滤器中，大气过滤器（第二次反演中）。该过滤器的低通滤波是对应于地面上的距离（通常为默认的1.2公里）：所以，感兴趣区域应该包含一个缓冲区2-3公里左右的真实感兴趣的区域。在任何情况下，研究区不应

16、小于2300*2300像素。如果用户在SBAS中用了incremental 模式，在之前的PS处理中又是对感兴趣区域进行的处理，那么需要重新对新加入的数据裁剪一下。在本文中，样例数据是经过裁剪的。干涉工作流这一步会自动进行以下处理：配准、干涉图生成、去平、振幅离差指数计算。l 在配准的一步，所有的图像都会配准到主影像上。在距离向上会过采样4倍（PreferencesPersistent ScatterersRange Looks setting 设置），这样可以避免由于基线太长而产生的快速变化的干涉条纹，这和DInSAR是不同的，因为PS是针对点目标进行的分析，不会执行频谱偏移和常见的多普勒带

17、宽滤波。每个图像都和主影像生成了干涉图。l 使用输入的DEM或者平均高程进行干涉相位的去平处理，如果输入“Geometry GCP file”，可以对主影像和DEM进行几何校正。参考DEM的精度越高，分辨率越高，地形去除效果越好。值得注意的是，如果主图像已经过手动或自动的轨道校正（general toolsOrbital Correction）就不需要输入GCP信息了。l 振幅离差指数（D）的计算公式为：，是时序上的振幅的标准差，是时序上的某个点的振幅平均值，如果一个象元，在所有的时相上，都能有相似的、比较大的振幅信息，那么这个象元的相位离差是比较小的。这样可以识别潜在的相干点，而不需要分析相

18、位，因为在这个阶段的相位仍然包含很多未知的信息，很难分析。振幅离差指数D越小，相位的标准差越大（如自然地表的郊区、植被区域的等），振幅离差指数越大，相位标准差越小（如城市区域）。因此，用离差指数的阈值进行选点是最实用的方法，相应的相位离差比较小。SAR数据会自动进行辐射定标，从而得到精确的和。如上所述，PS干涉过程会有一个过采样的处理，生成干涉结果。不过由于几何形状，这些数据不能直接进行可视化显示。如果要看快视图的话，可以输出多视后的去平干涉图dint，能看到在距离向和方位向一样分辨率的快视图。多视的视数是和主影像相关的。视数最好超过1：1，原因是：l 可以提高干涉的信噪比（SNR），以及得到

19、更可行的相干性估算结果；l 可以提高处理效率。如果增加多视的视数的话，象元的空间分辨率会降低，要估算最佳的视数，可以打开主影像的.sml文件，找到象元采样间隔和入射角的字段.ASAR 处理中，打开主影像Asar_slc_pwr.sml：7.41699 4.18004 22.3999GroundRangePixelSpacing = PixelSpacingRg / sin(IncidenceAngle) = 20.114要得到20米的多视后的分辨率： Multi-look range = 1 Multi-look azimuth = groundRangePixelSpacing/PixelS

20、pacingAz = 4.048 4图10 ASAR的多视视数PALSAR处理中，打开主影像Palsar_FBS_slc_pwr.sml，查看字段：4.00002 3.68399 38.9998groundRangePixelSpacing = PixelSpacingRg / sin(IncidenceAngle) = 7.422要得到多视后15 米的象元分辨率，视数为: Multi-look range = 1 Multi-look azimuth = groundRangePixelSpacing/PixelSpacingAz = 4.878 5图11 PALSAR的多视视数注：在PS干

21、涉处理开始时，应该输入工程文件、参考DEM，和正确的视数（便于生成快视图），在本处理中，ASAR视数是方位向4，距离向1，PALSAR视数是方位向5距离向1。因为输出的结果非常多，程序生成了三个meta文件，可以一次性打开所有相关的结果：l slant_dint_meta，包含所有的去已知地形后的干涉图 l slant_pwr_meta，包含所有的强度数据 注：在处理过程中，不要移动文件夹中的任何文件。第一次反演这一步是第一次模型反演估算残余地形和形变速率。在干涉图生成之后，就从所有的干涉图中去除了相位偏移，程序基于振幅离差指数自动选出一个或者多个像素（参考点），计算出相位偏移，参考点的数量取

22、决于感兴趣区的大小，默认的是25平方公里一个参考点（PreferencesPersistent ScatterersArea for Single Reference Point）。该方法是基于一定数量的“永久散射体”（Persistent Scatterers）的识别。然后集中对这些可靠的独立目标（像素）的历史相位进行分析。符合PS点的条件是：它必须是稳定的（波动小于1毫米/年），它必须是能被SAR天线探测到的方向。有两点必须说明：l 在研究区内的所有PS点的恒定形变是监测不出来的。l 适用于线性趋势的形变，形变速率的变化是监测不出来的好的PS点，如屋顶、电杆、桥梁，通常是城市中的点，或其他

23、人造地物，如温室、大坝、金属和混凝土的地物（如井盖周边，管道等）。除了这些之外，自然目标如裸露的岩石也可以是PS点。SAR数据的时相分布应该和预期的形变力学特点相一致。对于这些点，会用以下两种方法进行计算：l 一定范围内（Area for Single Reference Point，默认是25平方公里），只使用一个参考点。l 对大范围分析时，整个区域分为多个子区域，设置个重叠的百分比（Overlap for SubAreas%），每一个子区域都用的同一组参数，每个子区域独立处理。最后，会对所有子区的结果进行镶嵌，输出一个完整的结果。注：“area for Single Reference P

24、oint”定义的是一个子区域的范围。和SBAS不同的是，PS只有一个模型：线性模型，估算残余高程和形变速率。模型是：，Disp是时间t的形变量，V是形变速率。在第二次反演之前，检查第一次反演的结果，看镶嵌的结果是否有错误。实际上，在对velocity_first”和“heigh_first”镶嵌时，在重叠区会估算出一个偏移量，会用于每一块子区域的相位重计算，这个偏移量是用Coherence Threshold for Merging这个参数估算的。在ENVI中加载“velocity_first”和“heigh_first”这两个数据结果，再打开subareas.shp文件，首先，看速率和高程残

25、差有没有明显的像素值的跳跃，有时候会出现跳跃的情况，一些小的跃变会在大气滤波时去除。注：用彩色显示这两个结果，可以看出像素值跃变的情况。 图12 ASAR处理得到的Velocity_first (左) 和 Heigh_first (右)图13 PALSAR处理得到的Velocity_first (左) 和Heigh_first (右) 从上图可以看出，PALSAR数据得到的速率图中有一块值是有跳跃的，不过，因为这个区域比较小，这个结果可以用于第二次反演。注：在“cc_first”图上，可以看到不同子区域的重叠区，选择的是相干性更好的，这会导致不同子区域之间有值的跃变，这是正常的，如图14和15

26、。图14 ASAR处理得到的cc_first结果，有一些垂直方向上的相干性跃变图15 PALSAR处理得到的cc_first结果，中间有一块相干性跃变 如果在速度和残余高程图上发现不同子区域的值有跃变的现象，用户应该调整Area for single reference point、Overlap for SubAreas、和/或Coherence Threshold for Merging这几个参数，直到得到连续的结果。针对速率或高程有跃变的情况，首先可以调整Coherence Threshold for Merging这个参数（在other参数中），如果处理区域是城市区域，则可以把这个值调

27、小，这样可以增加用于偏移估计的点的密度（数量）。如果处理区域大多是植被覆盖的自然表面，可以调大这个参数，避免太多的无效信息。如果调整这个参数并没有得到好一些的结果，可以调整Area for single reference point和overlap for SubAreas这两个参数（减小分块大小或者是增加重叠区百分比），从而在重叠区能得到更高的相干性，不过需要注意的是：分块数量越多，处理时间越长。注：在进行第一次反演时，输入工程文件Auxiliary.sml，在ASAR处理时，Area for Single Reference Point (sqkm) 保持默认的5km，Coherence

28、 Threshold for Merging 保持默认的0.66。在PALSAR处理中，也用的是默认参数。图16 ASAR处理参数图17 PALSAR处理参数第二次反演第二次反演是对第一次的线性模型反演得到的结果估算大气的相位。进行线性模型估算，是在大气校正之前，从所有的差分干涉相位中估算出了速率和残余高程。在大气估算过程中，PS技术利用散射体的稠密性分布，以去除大部分信号的传播延迟的波动，这主要是由于在对流层中的变化。大气滤波是用时间上的高通滤波和空间上的低通滤波来实现的：大气低通滤波：和大气的空间特性有关，大气在空间上是高相干性的。是用一个窗口进行滤波，大规模的空间上的变化适合用大的窗口，

29、小窗口是更好地纠正孤立的工件，由于局部变化窗口越小，滤波效果越强。大气高通滤波：和大气的时间特性有关。大气在时间上是低相干的。是通过使用一个时间窗口，大窗口更适合纠正与低时间变化的影响，而小窗口更好地纠正频繁的大气变化。窗口尺寸越大，滤波效果越强。第二次反演的模型，得到的是去除大气之后的相位组分，得到最终的形变速率。这一步是必须做的。注：运行这一步，输入工程文件Auxiliary.sml，设置大气滤波的窗口。图18 ASAR处理的参数图19 PALSAR处理的参数这一步得到的结果有：Height：大气校正后的地形残差（米）precision_height：高程残差的平均精度（更多信息可参考相位

30、转高程）Velocity：大气校正之后的平均形变速率（毫米/年）precision_vel：平均形变速率的精度（更多信息可参考相位转形变）cc：多时相的相干系数。显示了有多少形变趋势和模型是相符合的此外，还生成了几个meta结果文件：slant_atm_meta：这是每个时相的大气组分，斜距坐标系的结果，在work文件夹下有这个文件。slant_dint_reflat_meta：这是大气校正之后的每个时相的去地形之后的干涉图，斜距坐标系的结果。slant_disp_meta：这是大气校正之后的每个时相的形变量，斜距坐标系的结果。地理编码PS结果产品进行地理编码，用户可以选择栅格和/或矢量两种格

31、式输出结果。为了获得更精确的形变结果，可以输入一个或多个地面控制点（例如从GPS或其他测量手段得到的）“Displacement GCP file”。这些点可用来优化形变估算的结果。如果只输入一个控制点，会计算出一个恒定的平均形变速率的偏移量，没有任何空间变化；如果输入更多的控制点，会用这些GCP拟合出一个最佳的偏移量。这些地面控制点必须是地理坐标。不过这种优化不是必须做的。Product Coherence Threshold参数的意思是，小于该阈值的象元将会去除，不会输出这些点。本文进行地理编码处理，输入工程文件Auxiliary.sml，输入参考DEM文件，其他参数保持默认，make G

32、eocoded Shape改为Ture：会得到矢量文件的PS点，矢量文件可以用时序分析工具进行分析，本处理中没有使用GPS数据。图20 ASAR的处理面板图21 PALSAR的处理面板地理编码得到的产品：l Ref_GCP_geo： 矢量文件，在反演step 1中自动生成的参考点。l SubAreas_geo：转成地理坐标的子区域分布。l mean_geo：SAR平均强度图像和附带文件(.sml, .hdr). l Meta files (_meta)：索引文件，包括形变速率、残余高度、相关系数图、每期形变量（需要勾选Geocoded Raster Products选项）。 l Meta fi

33、les (_meta)：索引文件，包括形变沿着最大坡度方向投影结果(_SD)和垂直投影面投影结果(_VD)。l work 子文件夹：存储中间结果。l PS渲染图和相关信息（.shp和.kml）。l Maximum slope direction values（_ADF）：最大坡度方向值（）及附带文件 (.sml, .hdr)。l Maximum slope inclination values（_IDF） ：最大坡度倾斜值及附带文件 (.sml, .hdr)。l Azimuth Line of Sight （_ALOS)：方位向视线及附带文件 (.sml, .hdr)。正角度表示 从北顺时针方

34、向，负角度表示从北逆时针方向。l Incidence angle of the Line of Sight（_ILOS）：入射角视线及附带文件 (.sml, .hdr)。表示视线和垂直方向的夹角。一般的产品索引文件包括：l _geo_vel+height_meta： 索引文件，地理坐标下的高度和形变速率测量值。 l _geo_otherinfo_meta：索引文件，地理坐标下强度均值、多期相干系数、高度量测精度和高度修正值。l _geo_disp_first_meta：索引文件，大气校正前，在地理坐标下每期形变量。l _geo_disp_meta：索引文件，大气校正后，在地理坐标下每期形变量。

35、在地理编码之后，PS处理就结束了，这时候如果想移动数据的话，就可以移动PS处理的整个文件夹_ps_processing了。注：如果选择make geocoding raster选项的话，会输出最大坡度方向_SD和垂直方向_VD的形变结果。时序分析可以用时序分析工具显示PS点在各时相上的形变量，如果选择了Make Geocoded Shape 为Ture的话，可以用General ToolsVector analyzer工具进行时序分析，如果选择了Make Geocoded Raster为Ture的话，可以用General ToolsRaster Analyzer来进行时序分析。本文只生成了矢量

36、结果。要对PS结果进行时序分析的话，首先在ENVI中加载mean_geo文件和PS点的.shp文件，然后加载显示，打开工具General toolsTime Series AnalyzerVector tool，选择数据范围，点击Color apply，在图中选择一个PS点，点击Plot Time Series，就可以打开历史形变图。如果有好几个.shp图层的话，在要分析的图层上点击鼠标右键，选择”Set as active layer”，图层在激活状态，再进行时序分析。注：（1）对不同的.shp图层选择同样的数据范围（2）要先加载显示mean_geo图像，然后再使用时序分析工具。图22 彩色

37、显示PS_75_0.shp中的速率字段，自定义范围为-15到15mm/y图23 可以在图层上面右键，把该图层设置为活动状态，时序分析工具就可以应用于该图层了结果说明可以在图24和图26看出，浦安市新填造的区域（A区）显示了明显的沉降，而在内陆地区（B区），表面稳定。而在时间序列上，ASAR得到的时序图比PALSAR更稳定；事实上，ASAR输入的时相比PALSAR的多。此外，两个数据的波长不同，对地面的探测也会有差异。可以看到，两组数据得到的结果是一致的。特别是两个沉降区。建议在可能的情况下，可以用不同传感器的多组数据，来对比结果。注：一个shp文件中最多能存储99999个点，如果PS点特别多，超过了这个数，会生成多个shp文件，如果不想输出这么多的PS点，可以把Product Coherence Threshold这个阈值设置小，PS点就不会特别多，可以生成在一个shp文件中，不过，PS点的空间分辨率会降低。 有时候因为电脑显卡原因，不能对远程桌面连接的数据进行颜色渲染。这种情况下，直接连接到远程位置并重新启动计算机或将结果复制