基于PS-InSAR技术的形变监测分析流程
文章目录
1. 背景知识
1.1 PS-InSAR技术
PS-InSAR技术是“永久散射体合成孔径雷达干涉测量”的多重嵌套缩写，即Persistent Scatterer Interferometric Synthetic Aperture Radar。其中PS（永久散射体）指对雷达波的后向散射较强，并且在时序上较稳定的各种地物目标，如建筑物与构筑物的顶角、桥梁、栏杆、裸露的岩石等目标；InSAR(合成孔径雷达干涉测量)是指利用同一地区不同期次SAR数据中的相位信息进行干涉测量的技术。
1.1.1 雷达干涉测量
雷达影像反映了雷达所发射的电磁波和目标物相互作用的结果。雷达干涉测量技术，综合了合成孔径雷达（SAR）成像原理和干涉测量技术，利用传感器的系统参数、姿态参数和轨道之间的几何关系等精确测量地表某一点的三维空间位置及其微小变化。SAR本身是一种主动式微波传感器，由于其全天候、全天时获取数据，并能穿透云雾、烟尘和大面积获取地表信息的特点而成为对地观测领域不可或缺的传感器，尤其适用于传统光学传感器成像困难的地区。
1.1.2 InSAR技术
InSAR技术以合成孔径雷达复数据提取的相位信息为信息源获取地表的三维信息和变化信息。InSAR通过两副天线同时观测（单轨模式），或两次近平行的观测（重复轨道模式），获取地表同一景观的复图像对。由于目标与两天线位置的几何关系，在复图像上产生了相位差，形成干涉条纹图。干涉条纹图中包含了斜距向上的点与两天线位置之差的精确信息。因此，利用传感器高度、雷达波长、波束视向及天线基线距之间的几何关系，可以精确地测量出图像上每一点的三维位置和变化信息。
1.1.3 技术原理
PS-InSAR技术中的PS（永久散射体）指对雷达波的后向散射较强，并且在时序上较稳定的各种地物目标，如建筑物与构筑物的顶角、桥梁、栏杆、裸露的岩石等目标；InSAR(合成孔径雷达干涉测量)是指利用同一地区不同期次SAR数据中的相位信息进行干涉测量的技术。
基于InSAR技术，对K+1幅SAR单视复数影像，经配准、辐射定标、PS探测和干涉处理，并借助已知DEM进行差分干涉处理，得到K幅干涉和差分干涉图、H个PS点以及各PS点在各差分干涉图中的差分干涉相位集。在考虑地表形变、高程误差、大气影响及失相关的情况下，得到每个PS点在每幅差分干涉图上的差分干涉相位组成，其中，对形变速率增量和高程误差增量积分，可以得到每个PS点相对于主参考点的形变速率和高程误差。同时，根据求解结果在PS离散点上进行相位解缠，经过积分，还可以获得解缠的线性相位残差（相对于主参考点）。
1.1.4 技术特征
利用雷达卫星进行PS-InSAR干涉测量，具有以下特征：
无需地面测站
由于雷达卫星干涉测量监测无需地面测站，因而可使监测时空范围的设计更为自由、方便。同时，可以避免地面控制点的限制，尤其是许多中间过渡点（采用常规大地测量方法进行变形监测时，为传递坐标经常要设立许多中间过渡点），且不必建标，从而可节省大量的人力物力，大大提高监测效率。
主动发射微波
雷达卫星干涉测量由地面控制站根据监测任务安排，制定卫星数据获取计划，卫星根据编程指令，绕行地球通过制定区域时，向地面主动发射微波并接收回波完成测量。
观测点密度高
常规监测条件下，1平方公里内的监测点数量一般为1-100个，离散孤立的监测点，仅能近似地反映区域形变的情况。雷达干涉测量监测点数平均密度可达20000个/平方公里，高密度分布的观测点，为观测区域内不同目标的形变分析提供客观数据支持，进而实现区域内连续形变特征分析。
全天候观测
雷达卫星干涉测量不受气候条件的限制，在夜晚或是风雪雨雾条件下仍能进行有效观测。这一点对于汛期、多云多雨地区的崩塌、滑坡、泥石流等地质灾害监测是非常有利的。
全自动化观测
由于雷达卫星干涉测量的数据采集工作是自动进行的，同时卫星与接收站、接收站与用户之间通过数据链路进行联系，故用户可以较为方便地把雷达卫星干涉测量监测系统建成全自动化的监测系统。这种系统涉及不但可保证长期连续运行，而且可大幅度降低变形监测成本，提高监测资料的可靠性。
mm级精度
mm级的精度已可满足一般崩滑体变形监测的精度要求，因而可在滑坡、崩塌、泥石流等地质灾害的监测中得到了广泛的应用，成为一种新的有效的监测手段。
1.1.5 技术优化
利用雷达卫星进行干涉测量进行监测时也存在一些不足之处，主要表现在，大气参数的变化（对流层水汽含量和电离层）、地形变化剧烈或植被覆盖茂密区域的去相关引起的相位噪声及失相干，复杂地形条件下的相位解缠，轨道参数（基线）等的精确校准和地形快速纠正等问题。
1.1.6 应用
鉴于Ｄ-InSAR算法的缺陷，意大利雷达遥感专家Alessandro Ferretti于1999年提出了PS-InSAR方法。其基本思想是：第一，利用覆盖同一研究区的多景单视ＳＡＲ影像，选取其中一景ＳＡＲ影像作为主影像，其余ＳＡＲ影像分别与主影像配准，依据时间序列上的幅度和（或）相位信息的稳定性选取永久散射体（Persisttent Scatterer, PS）目标；第二，经过干涉和去地形处理，得到基于永久散射体目标的差分干涉相位，并对相邻的永久散射体目标的差分干涉相位进行再次差分；第三，根据两次差分后的干涉相位中各个相位成分的不同特性，采用构建形变相位模型和时空滤波或方式估计形变和地形残余信息。
PS－InSAR技术不是针对ＳＡＲ 影像中的所有像元进行数据处理，而是选取在时间上散射特性相对稳定、回波信号较强的ＰＳ点作为观测对象。这些ＰＳ点通常包括人工建筑物、灯塔、裸露的岩石以及人工布设的角反射器等。ＰＳ点的准确选取可以确保即便在干涉对的时间或空间基线很长的条件下（甚至达到临界基线），ＰＳ点依然呈现出较好的相干性和稳定性。PS－InSAR技术已在多个城市的高分辨率地面沉降监测中得到广泛应用，特别是城市重点基础设施的高分辨率形变监测。通过对比同期的水准和GPS测量数据，证实了PS－InSAR技术具有较高的可靠性,其精度可以到mm级。
然而PS－InSAR方法也存在自身的缺陷，主要表现在两个方面：第一，其通常要求覆盖同一区域的ＳＡＲ影像数目较多（通常要求大于２５景），便于保证模型解算的可靠性。其次，PS－InSAR技术由于是基于大量ＰＳ点的迭代回归或网平差计算，运算效率不高，因此不适合大范围地区。
1.2 Sentinel-1数据
1.2.1 Sentinel-1简介
哨兵1号（sentinel-1）包括哨兵-1A和哨兵-1B两颗卫星。这两颗卫星是处于同一轨道平面的极轨卫星，分别于2014年4月3日和2016年4月25日成功发射。这两颗卫星搭载C波段合成孔径雷达，具有4种成像模式，可为陆地和海洋服务提供全天时、全天候的雷达图像，提供一系列运营服务，包括北极海冰，日常海冰测绘，海洋环境监视监测科研，监测地面运动风险，森林制图，水和土壤管理和测绘，以支持人道主义援助和危机情况。
1.2.2 Sentinel-1扫描模式
SM（Stripmap）：一种标准的SAR条形图成像模式，其中地面区域被连续的脉冲序列照亮，而天线波束指向一个固定的方位角和仰角。SM模式仅用于小岛屿，在紧急情况管理等特殊事件时使用。
IW（Interferometric Wide swath）：IW模式是陆地上的主要采集模式，满足了大部分业务需求。它以5米 x 20米的空间分辨率(单视)获取250公里长的数据。IW模式使用渐进扫描SAR (TOPSAR)地形观测捕获三个子区域。在TOPSAR技术中，除了像扫描雷达一样控制波束的范围外，波束还可以在每个爆发的方位角方向上由后向前进行电子控制，避免了扇形现象，并导致整个区域的图像质量均匀。
EW（Extra Wide swath）：使用TOPSAR成像技术在五个区域获取数据。EW模式以牺牲空间分辨率为代价提供了非常大的区域覆盖。(言外之意是空间分辨率低)。EW模式主要用于沿海监测，包括海运监测、溢油监测和海冰监测。
WV（Wave）：数据是在被称为“小片段”的小型条形地图场景中获取的，这些场景在轨道沿线每隔100公里定期设置一次。通过交替获得小点，以近距离入射角获得一个小点，而以远距离入射角获得下一个小点。WV是哨兵1号在海上的操作模式。（言外之意用于海洋）。
对于每一种模式，都可以生产SAR的0级、1级SLC、1级GRD和2级OCN的产品。
1.2.3 Sentinel-1数据产品
- Raw Level-0 data （特定情况下使用）：0级产品；
- SLC（ Single Look Complex）：已被处理后的一级产品，能获得相位和振幅信息。相位信息是时间的函数，根据相位信息和速度可实现距离的测量。（可用于测距和形变观测）。
- GRD（Ground Range Detected）：一级产品，有多视强度数据，该强度数据与后向散射系数有关。（可用于土壤水分反演）。
- OCN（Ocean）：二级产品，用于检索海洋地球物理参数。（即应用于海洋）。
所有的产品都是从0级产品直接加工的。每种模式都可以潜在地生成一级SLC、一级GRD和二级Ocean产品。
1.2.4 Sentinel-1应用
1.2.5 Sentinel-1极化方式
- 单极化方式：HH或VV；
- 双极化方式：HH+HV 或 VV+VH；
- WV扫描模式只有单极化方式HH或VV。其他扫描模式单极化方式（HH或VV）和双极化方式（HH+HV或VV+VH）都有。
- IW：用（VV+VH）极化方式–>观测陆地；
- WV: 用VV极化方式–>观测海洋
1.2.6 Sentinel-1产品分辨率
-
SLC一级产品分辨率
-
GRD一级产品分辨率
1.2.7 Sentinel-1数据文件命名格式
- MMM：表示数据来源于A星或B星，有S1A和S1B两个选择。
- BB：条带扫描模式，有IW、EW、WV 3种选择。
- TTT：表示产品的类型，有SLC、GRD、OCN 3种产品选择。
- R： 为分辨率类别。F表示(Full resolution),H表示High resolution，M表示Medium resolution(分辨率类别仅用于GRD)。
- L：数据处理等级，为1级，或2级。
- F：产品类可以是Standard (S)或Annotation (A)。Annotation产品只在PDGS内部使用，不分发。
- PP：极化方式，具体见图。
- 中间一串是起始时间和数据截止时间。
- OOOOOO：产品开始时的绝对轨道号，轨道号范围：000001-999999。
- DDDDDD：任务数据获取标识符,范围:000001-FFFFFF。
- CCCC：产品唯一标识符,是使用CRC-CCITT在清单文件上计算CRC-16生成的十六进制字符串。
在产品文件夹中，测量数据集和注释数据集遵循类似的命名约定，用破折号(-)分隔小写字母和数字字符。
1.2.8 Sentinel-1数据下载
官网下载地址：https://sentinel.esa.int/web/sentinel/toolboxes/sentinel-1 (数据检索不方便，无法批量下载）。
其他下载地址：https://search.asf.alaska.edu/（地球数据网站，推荐）
处理哨兵数据的工具：
-
GEE（Google Earth Engine）：https://earthengine.google.com
-
SNAP（Sentinel Application Platform）：https://step.esa.int/main/download/snap-download/
1.3 SNAP软件
1.3.1 SNAP介绍
A common architecture for all Sentinel Toolboxes is being jointly developed by Brockmann Consult, SkyWatch and C-S called the Sentinel Application Platform (SNAP).
SNAP全称Sentinel Application Platform，是欧空局开发的卫星数据科学探索通用工具箱。SNAP 为处理、建模和可视化卫星图像提供了一个直观的平台，特别是对于哨兵任务。该软件经过优化，可以处理大量卫星数据。它还有助于处理 SAR 数据，例如来自 Sentinel-1 的数据。
1.3.2 SNAP的特点
- 所有工具箱的通用架构；
- 快速的图像显示和导航，即使是千兆像素图像；
- 图形处理框架（GPF）：用于创建用户定义的处理链；
- 高级图层管理：允许添加和操作新的叠加层，例如其他波段的图像、来自 WMS 服务器的图像或 ESRI shapefile；
- 用于自定义感兴趣区统计和各种丰富的绘图；
- 简单的位掩码定义和覆盖；
- 使用任意数学表达式的灵活频带计算
- 对常见地图投影进行准确的重投影和正射校正
- 使用地面控制点进行地理编码和校正
- 自动下载SRTM DEM 和支持切片选择
- 用于高效扫描和编目大型档案的产品库；
- 支持多线程和多核处理器；
- 集成可视化的WorldWind；
1.3.3 SNAP 所用到的技术
- NetBeans平台 —桌面应用程序框架
- Install4J —跨平台安装
- GeoTools —地理空间工具库
- GDAL —读/写栅格和矢量地理空间数据
- Jira —问题跟踪器
- Git —版本控制系统
1.4 名词对照
| 名词 | 说明 | |
|---|---|---|
| coregistration | 配准 | |
| TOPS | 渐进扫描地形观测工作模式，Terrain Observation by Progressive Scans | |
| SLC | 单视复数影像，Single Look Complex |
2. 数据准备
至少15幅同一区域不同时段的Sentinel-1 IW采集模式下的SLC产品数据，并且极化方式为DV，示例数据如下：
S1A_IW_SLC__1SDV_20201015T004834_20201015T004901_034800_040E3E_59DE.zip
S1A_IW_SLC__1SDV_20201108T004834_20201108T004901_035150_041A4C_59CC.zip
S1A_IW_SLC__1SDV_20201226T004832_20201226T004859_035850_043283_3255.zip
S1A_IW_SLC__1SDV_20210119T004831_20210119T004858_036200_043EBF_7BD0.zip
S1A_IW_SLC__1SDV_20210212T004831_20210212T004858_036550_044AE4_13B8.zip
S1A_IW_SLC__1SDV_20210308T004830_20210308T004857_036900_045718_FB8A.zip
S1A_IW_SLC__1SDV_20210401T004831_20210401T004858_037250_04633D_EE19.zip
S1A_IW_SLC__1SDV_20210425T004832_20210425T004859_037600_046F51_41DC.zip
S1A_IW_SLC__1SDV_20210519T004833_20210519T004900_037950_047A9D_E8F1.zip
S1A_IW_SLC__1SDV_20210612T004834_20210612T004901_038300_04850E_91FF.zip
S1A_IW_SLC__1SDV_20210706T004836_20210706T004903_038650_048F8D_4242.zip
S1A_IW_SLC__1SDV_20210730T004837_20210730T004904_039000_0499FE_DFFC.zip
S1A_IW_SLC__1SDV_20210823T004838_20210823T004905_039350_04A5B8_7F06.zip
S1A_IW_SLC__1SDV_20210928T004840_20210928T004907_039875_04B7B7_EC4F.zip
S1A_IW_SLC__1SDV_20211010T004840_20211010T004907_040050_04BDBF_F498.zip
3. 软件准备
3.1 SNAP
需要利用SNAP软件对SLC数据进行影像分割（TOPS Split）、轨道校正（Apply Orbit File）、主影像获取（InSAR Stack Overview）、影像配准（Sentinel-1 Back Geocoding）、条带Deburst（Sentinel-1 TOPSAR Deburst and Merge）、干涉处理（Interferogram formation）、StaMPS Export等数据预处理（Pre-processing）操作。
SNAP安装步骤：
系统环境：Ubuntu 20.04
安装包下载：https://step.esa.int/main/download/snap-download/
参考文档：http://step.esa.int/docs/tutorials/SNAP_CommandLine_Tutorial.pdf
注：安装前注意关闭终端的X11转发，否则可能会在调用本地的图形界面时出现问题（在使用云服务器的堡垒机连接时出现此问题，直连时没有）。
安装命令：
chmod +x esa-snap_sentinel_unix_9_0_0.sh
sh esa-snap_sentinel_unix_9_0_0.sh
# 安装完成后配置全局变量
vi ~/.bashrc
# 在文件的最后一行添加一行:export PATH=$PATH:bin目录路径,例如
export PATH=$PATH:/usr/local/snap/bin
# 使配置生效
source ~/.bashrc
# 命令测试
gpt -h
执行gpt -h后，如果出现如下内容，说明已经安装成功：
INFO: org.esa.snap.core.gpf.operators.tooladapter.ToolAdapterIO: Initializing external tool adapters
INFO: org.esa.s2tbx.dataio.gdal.GDALVersion: GDAL not found on system. Internal GDAL 3.2.1 from distribution will be used. (f0)
INFO: org.esa.s2tbx.dataio.gdal.GDALVersion: Internal GDAL 3.2.1 set to be used by SNAP.
INFO: org.esa.snap.core.util.EngineVersionCheckActivator: Please check regularly for new updates for the best SNAP experience.
INFO: org.esa.s2tbx.dataio.gdal.GDALVersion: Internal GDAL 3.2.1 set to be used by SNAP.
Usage:
gpt <op>|<graph-file> [options] [<source-file-1> <source-file-2> ...]
Description:
This tool is used to execute SNAP raster data operators in batch-mode. The
operators can be used stand-alone or combined as a directed acyclic graph
(DAG). Processing graphs are represented using XML. More info about
processing graphs, the operator API, and the graph XML format can be found
in the SNAP documentation.
Arguments:
<op> Name of an operator. See below for the list of <op>s.
<graph-file> Operator graph file (XML format).
<source-file-i> The <i>th source product file. The actual number of source
file arguments is specified by <op>. May be optional for
operators which use the -S option.
Options:
-h Displays command usage. If <op> is given, the specific
operator usage is displayed.
-e Displays more detailed error messages. Displays a stack
trace, if an exception occurs.
-t <file> The target file. Default value is 'target.dim'.
-f <format> Output file format, e.g. 'GeoTIFF', 'HDF5',
'BEAM-DIMAP'. If not specified, format will be derived
from the target filename extension, if any, otherwise the
default format is 'BEAM-DIMAP'. Ony used, if the graph
in <graph-file> does not specify its own 'Write'
operator.
-p <file> A (Java Properties) file containing processing parameters
in the form <name>=<value> or a XML file containing a
parameter DOM for the operator. Entries in this file are
overwritten by the -P<name>=<value> command-line option
(see below). The following variables are substituted in
the parameters file:
${system.<java-sys-property>}
${operatorName} (given by the <op> argument)
${graphFile} (given by the <graph-file> argument)
${targetFile} (pull path given by the -t option)
${targetDir} (derived from -t option)
${targetName} (derived from -t option)
${targetBaseName} (derived from -t option)
${targetFormat} (given by the -f option)
-c <cache-size> Sets the tile cache size in bytes. Value can be suffixed
with 'K', 'M' and 'G'. Must be less than maximum
available heap space. If equal to or less than zero, tile
caching will be completely disabled. The default tile
cache size is '1,073,741,824M'.
-q <parallelism> Sets the maximum parallelism used for the computation,
i.e. the maximum number of parallel (native) threads.
The default parallelism is '8'.
-x Clears the internal tile cache after writing a complete
row of tiles to the target product file. This option may
be useful if you run into memory problems.
-S<source>=<file> Defines a source product. <source> is specified by the
operator or the graph. In an XML graph, all occurrences of
${<source>} will be replaced with references to a source
product located at <file>.
-P<name>=<value> Defines a processing parameter, <name> is specific for the
used operator or graph. In an XML graph, all occurrences
of ${<name>} will be replaced with <value>. Overwrites
parameter values specified by the '-p' option.
-D<name>=<value> Defines a system property for this invocation.
-v <dir> A directory containing any number of Velocity templates.
Each template generates a text output file along with the
target product. This feature has been added to support a
flexible generation of metadata files.
See http://velocity.apache.org/ and option -m.
-m <file> A (Java Properties) file containing (constant) metadata
in the form <name>=<value> or any XML file. Its primary
usage is to provide an additional context to be used
from within the Velocity templates. See option -v.
--diag Displays version and diagnostic information.
Operators:
Aatsr.SST Computes sea surface temperature (SST) from (A)ATSR products.
AATSR.Ungrid Ungrids (A)ATSR L1B products and extracts geolocation and pixel field of view data.
AdaptiveThresholding Detect ships using Constant False Alarm Rate detector.
AddElevation Creates a DEM band
AddLandCover Creates a land cover band
ALOS-Deskewing Deskewing ALOS product
Apply-Orbit-File Apply orbit file
Arc.SST Computes sea surface temperature (SST) from (A)ATSR and SLSTR products.
ArviOp Atmospherically Resistant Vegetation Index belongs to a family of indices with built-in atmospheric corrections.
Azimuth-Shift-Estimation-ESD Estimate azimuth offset for the whole image
AzimuthFilter Azimuth Filter
Back-Geocoding Bursts co-registration using orbit and DEM
BandMaths Create a product with one or more bands using mathematical expressions.
BandMerge Allows copying raster data from any number of source products to a specified 'master' product.
BandPassFilter Creates a basebanded SLC based on a subband of 1/3 the original bandwidth
BandsDifferenceOp No description available.
BandSelect Creates a new product with only selected bands
BandsExtractorOp Creates a new product out of the source product containing only the indexes bands given
Bi2Op The Brightness index represents the average of the brightness of a satellite image.
This index is sensitive to the brightness of soils which is highly correlated with the humidity and the presence of salts in surface
Binning Performs spatial and temporal aggregation of pixel values into cells ('bins') of a planetary grid
BiOp The Brightness index represents the average of the brightness of a satellite image.
Biophysical10mOp The 'Biophysical Processor' operator retrieves LAI from atmospherically corrected Sentinel-2 products
BiophysicalLandsat8Op The 'Biophysical Processor' operator retrieves LAI from atmospherically corrected Landsat8 products
BiophysicalOp The 'Biophysical Processor' operator retrieves LAI from atmospherically corrected Sentinel-2 products
c2rcc.landsat8 Performs atmospheric correction and IOP retrieval with uncertainties on Landsat-8 L1 data products.
c2rcc.meris Performs atmospheric correction and IOP retrieval with uncertainties on MERIS L1b data products.
c2rcc.meris4 Performs atmospheric correction and IOP retrieval with uncertainties on MERIS L1b data products from the 4th reprocessing.
c2rcc.modis Performs atmospheric correction and IOP retrieval on MODIS L1C_LAC data products.
c2rcc.msi Performs atmospheric correction and IOP retrieval with uncertainties on Sentinel-2 MSI L1C data products.
c2rcc.olci Performs atmospheric correction and IOP retrieval with uncertainties on SENTINEL-3 OLCI L1B data products.
c2rcc.seawifs Performs atmospheric correction and IOP retrieval on SeaWifs L1C data products.
c2rcc.viirs Performs atmospheric correction and IOP retrieval on Viirs L1C data products.
Calibration Calibration of products
Change-Detection Change Detection.
ChangeVectorAnalysisOp The 'Change Vector Analysis' between two dual bands at two differents dates.
CiOp Colour Index was developed to differentiate soils in the field.
In most cases the CI gives complementary information with the BI and the NDVI.
Used for diachronic analyses, they help for a better understanding of the evolution of soil surfaces.
CloudProb Applies a clear sky conservative cloud detection algorithm.
Coherence Estimate coherence from stack of coregistered images
Collocate Collocates two products based on their geo-codings.
Compactpol-Radar-Vegetation-Index Compact-pol Radar Vegetation Indices generation
Compute-Slope-Aspect Compute Slope and Aspect from DEM
Convert-Datatype Convert product data type
CoregistrationOp Coregisters two rasters, not considering their location
CP-Decomposition Perform Compact Polarimetric decomposition of a given product
CP-Simulation Simulation of Compact Pol data from Quad Pol data
CP-Stokes-Parameters Generates compact polarimetric Stokes child parameters
CreateStack Collocates two or more products based on their geo-codings.
Cross-Channel-SNR-Correction Compute general polarimetric parameters
Cross-Correlation Automatic Selection of Ground Control Points
CrossResampling Estimate Resampling Polynomial using SAR Image Geometry, and Resample Input Images
DarkObjectSubtraction Performs dark object subtraction for spectral bands in source product.
DeburstWSS Debursts an ASAR WSS product
DecisionTree Perform decision tree classification
DEM-Assisted-Coregistration Orbit and DEM based co-registration
Demodulate Demodulation and deramping of SLC data
Double-Difference-Interferogram Compute double difference interferogram
DviOp Difference Vegetation Index retrieves the Isovegetation lines parallel to soil line
EAP-Phase-Correction EAP Phase Correction
Ellipsoid-Correction-GG GG method for orthorectification
Ellipsoid-Correction-RD Ellipsoid correction with RD method and average scene height
EMClusterAnalysis Performs an expectation-maximization (EM) cluster analysis.
Enhanced-Spectral-Diversity Estimate constant range and azimuth offsets for a stack of images
Faraday-Rotation-Correction Perform Faraday-rotation correction for quad-pol product
Fill-DEM-Hole Fill holes in given DEM product file.
FlhMci Computes fluorescence line height (FLH) or maximum chlorophyll index (MCI).
Flip flips a product horizontal/vertical
ForestCoverChangeOp Creates forest change masks out of two source products
FUB.Water MERIS FUB-CSIRO Coastal Water Processor to retrieve case II water properties and atmospheric properties
FuClassification Colour classification based on the discrete Forel-Ule scale
GemiOp This retrieves the Global Environmental Monitoring Index (GEMI).
Generalized-Radar-Vegetation-Index Generalized Radar Vegetation Indices generation
GenericRegionMergingOp The 'Generic Region Merging' operator computes the distinct regions from a product
GLCM Extract Texture Features
GndviOp Green Normalized Difference Vegetation Index
GoldsteinPhaseFiltering Phase Filtering
GRD-Post Applies GRD post-processing
HorizontalVerticalMotion Computation of Horizontal/Vertical Motion Components
IEM-Hybrid-Inversion Performs IEM inversion using Hybrid approach
IEM-Multi-Angle-Inversion Performs IEM inversion using Multi-angle approach
IEM-Multi-Pol-Inversion Performs IEM inversion using Multi-polarization approach
Image-Filter Common Image Processing Filters
Import-Vector Imports a shape file into a product
IntegerInterferogram Create integer interferogram
Interferogram Compute interferograms from stack of coregistered S-1 images
IonosphericCorrection Estimation of Ionospheric Phase Screens
IpviOp Infrared Percentage Vegetation Index retrieves the Isovegetation lines converge at origin
IreciOp Inverted red-edge chlorophyll index
KDTree-KNN-Classifier KDTree KNN classifier
KMeansClusterAnalysis Performs a K-Means cluster analysis.
KNN-Classifier K-Nearest Neighbour classifier
Land-Cover-Mask Perform decision tree classification
Land-Sea-Mask Creates a bitmask defining land vs ocean.
LandWaterMask Operator creating a target product with a single band containing a land/water-mask.
LinearToFromdB Converts bands to/from dB
Maximum-Likelihood-Classifier Maximum Likelihood classifier
McariOp Modified Chlorophyll Absorption Ratio Index, developed to be responsive to chlorophyll variation
Mci.s2 Computes maximum chlorophyll index (MCI) for Sentinel-2 MSI.
Merge Allows merging of several source products by using specified 'master' as reference product.
Meris.Adapt.4To3 Provides the adaptation of MERIS L1b products from 4th to 3rd reprocessing.
Meris.CorrectRadiometry Performs radiometric corrections on MERIS L1b data products.
Meris.N1Patcher Copies an existing N1 file and replaces the data for the radiance bands
Minimum-Distance-Classifier Minimum Distance classifier
MndwiOp Modified Normalized Difference Water Index, allowing for the measurement of surface water extent
Mosaic Creates a mosaic out of a set of source products.
MphChl This operator computes maximum peak height of chlorophyll (MPH/CHL).
Msavi2Op This retrieves the second Modified Soil Adjusted Vegetation Index (MSAVI2).
MsaviOp This retrieves the Modified Soil Adjusted Vegetation Index (MSAVI).
MtciOp The Meris Terrestrial Chlorophyll Index, aims at estimating the Red Edge Position (REP).
This is the maximum slant point in the red and near-infrared region of the vegetal spectral reflectance.
It is useful for observing the chlorophyll contents, vegetation senescence, and stress for water and nutritional deficiencies, but it is less suitable for land classification
Multi-size Mosaic Creates a multi-size mosaic out of a set of source products.
Multi-Temporal-Speckle-Filter Speckle Reduction using Multitemporal Filtering
Multilook Averages the power across a number of lines in both the azimuth and range directions
MultiMasterInSAR Multi-master InSAR processing
MultiMasterStackGenerator Generates a set of master-slave pairs from a coregistered stack for use in SBAS processing
Multitemporal-Compositing Compute composite image from multi-temporal RTCs
Ndi45Op Normalized Difference Index using bands 4 and 5
NdpiOp The normalized differential pond index, combines the short-wave infrared band-I and the green band
NdtiOp Normalized difference turbidity index, allowing for the measurement of water turbidity
NdviOp The retrieves the Normalized Difference Vegetation Index (NDVI).
Ndwi2Op The Normalized Difference Water Index, allowing for the measurement of surface water extent
NdwiOp The Normalized Difference Water Index was developed for the extraction of water features
Object-Discrimination Remove false alarms from the detected objects.
Offset-Tracking Create velocity vectors from offset tracking
Oil-Spill-Clustering Remove small clusters from detected area.
Oil-Spill-Detection Detect oil spill.
OlciAnomalyFlagging Adds a flagging band indicating saturated pixels and altitude data overflows
OlciO2aHarmonisation Performs O2A band harmonisation on OLCI L1b product. Implements update v4 of R.Preusker, June 2020.
OlciSensorHarmonisation Performs sensor harmonisation on OLCI L1b product. Implements algorithm described in 'OLCI A/B Tandem Phase Analysis'
Orientation-Angle-Correction Perform polarization orientation angle correction for given coherency matrix
Oversample Oversample the datset
OWTClassification Performs an optical water type classification based on atmospherically corrected reflectances.
PCA Performs a Principal Component Analysis.
PduStitching Stitches multiple SLSTR L1B product dissemination units (PDUs) of the same orbit to a single product.
PhaseToDisplacement Phase To Displacement Conversion along LOS
PhaseToElevation DEM Generation
PhaseToHeight Phase to Height conversion
PixEx Extracts pixels from given locations and source products.
Polarimetric-Classification Perform Polarimetric classification of a given product
Polarimetric-Decomposition Perform Polarimetric decomposition of a given product
Polarimetric-Matrices Generates covariance or coherency matrix for given product
Polarimetric-Parameters Compute general polarimetric parameters
Polarimetric-Speckle-Filter Polarimetric Speckle Reduction
PpeFiltering Performs Prompt Particle Event (PPE) filtering on OLCI L1B
Principle-Components Principle Component Analysis
ProductSet-Reader Adds a list of sources
PssraOp Pigment Specific Simple Ratio, chlorophyll index
PviOp Perpendicular Vegetation Index retrieves the Isovegetation lines parallel to soil line. Soil line has an arbitrary slope and passes through origin
Rad2Refl Provides conversion from radiances to reflectances or backwards.
Radar-Vegetation-Index Dual-pol Radar Vegetation Indices generation
Random-Forest-Classifier Random Forest based classifier
RangeFilter Range Filter
RayleighCorrection Performs radiometric corrections on OLCI, MERIS L1B and S2 MSI L1C data products.
Read Reads a data product from a given file location.
ReflectanceToRadianceOp The 'Reflectance To Radiance Processor' operator retrieves the radiance from reflectance using Sentinel-2 products
ReipOp The red edge inflection point index
Remodulate Remodulation and reramping of SLC data
RemoteExecutionOp The Remote Execution Processor executes on the remote machines a slave graph and then on the host machine it executes a master graph using the products created by the remote machines.
Remove-GRD-Border-Noise Mask no-value pixels for GRD product
RemoveAntennaPattern Remove Antenna Pattern
ReplaceMetadata Replace the metadata of the first product with that of the second
Reproject Reprojection of a source product to a target Coordinate Reference System.
Resample Resampling of a multi-size source product to a single-size target product.
RiOp The Redness Index was developed to identify soil colour variations.
RviOp Ratio Vegetation Index retrieves the Isovegetation lines converge at origin
S2repOp Sentinel-2 red-edge position index
S2Resampling Specific S2 resample algorithm
SAR-Mosaic Mosaics two or more products based on their geo-codings.
SAR-Simulation Rigorous SAR Simulation
SARSim-Terrain-Correction Orthorectification with SAR simulation
SaviOp This retrieves the Soil-Adjusted Vegetation Index (SAVI).
SetNoDataValue Set NoDataValueUsed flag and NoDataValue for all bands
SliceAssembly Merges Sentinel-1 slice products
SM-Dielectric-Modeling Performs SM inversion using dielectric model
SmacOp Applies the Simplified Method for Atmospheric Corrections of Envisat MERIS/(A)ATSR measurements.
SnaphuExport Export data and prepare conf file for SNAPHU processing
SnaphuImport Ingest SNAPHU results into InSAR product.
Speckle-Divergence Detect urban area.
Speckle-Filter Speckle Reduction
SpectralAngleMapperOp Classifies a product using the spectral angle mapper algorithm
SRGR Converts Slant Range to Ground Range
Stack-Averaging Averaging multi-temporal images
Stack-Split Writes all bands to files.
StampsExport Export data for StaMPS processing
StatisticsOp Computes statistics for an arbitrary number of source products.
SubGraph Encapsulates a graph within a graph.
Subset Create a spatial and/or spectral subset of a data product.
Supervised-Wishart-Classification Perform supervised Wishart classification
TemporalPercentile Computes percentiles over a given time period.
Terrain-Correction RD method for orthorectification
Terrain-Flattening Terrain Flattening
Terrain-Mask Terrain Mask Generation
ThermalNoiseRemoval Removes thermal noise from products
Three-passDInSAR Differential Interferometry
TileCache Experimental Operator which provides a dedicated cache for its source product.
A guide on how this operator is used is provided at https://senbox.atlassian.net/wiki/x/VQCTLw.
TileWriter Writes a data product to a tiles.
TndviOp Transformed Normalized Difference Vegetation Index retrieves the Isovegetation lines parallel to soil line
ToolAdapterOp Tool Adapter Operator
TopoPhaseRemoval Compute and subtract TOPO phase
TOPSAR-Deburst Debursts a Sentinel-1 TOPSAR product
TOPSAR-DerampDemod Bursts co-registration using orbit and DEM
TOPSAR-Merge Merge subswaths of a Sentinel-1 TOPSAR product
TOPSAR-Split Creates a new product with only the selected subswath
TsaviOp This retrieves the Transformed Soil Adjusted Vegetation Index (TSAVI).
Undersample Undersample the datset
Unmix Performs a linear spectral unmixing.
Update-Geo-Reference Update Geo Reference
Warp Create Warp Function And Get Co-registrated Images
WdviOp Weighted Difference Vegetation Index retrieves the Isovegetation lines parallel to soil line. Soil line has an arbitrary slope and passes through origin
Wind-Field-Estimation Estimate wind speed and direction
Write Writes a data product to a file.
3.2 MATLAB
需要利用MATLAB作为运行环境调用StaMPS脚本来提取地面位移变化，并根据结果矩阵绘制变化曲线。
MATLAB安装步骤：
系统环境：Ubuntu 20.04
操作终端：MobaXterm
参考文档：https://zhuanlan.zhihu.com/p/590500384
准备好的安装文件包括：
Matlab910R2021a_Lin64.iso和Crack文件夹，Crack文件夹里有4个文件：
libmwlmgrimpl.so、license.lic、license_server.lic和license_standalone.lic，文件上传到服务器/opt/matlab路径下
安装前注意是否已经有java环境，可以通过java -version测试，如果没有则可以通过apt install default-jdk -y安装最新版jdk，或者通过apt install openjdk-11-jdk安装指定版本。
3.2.1. iso文件挂载
- 创建一个文件夹，用于挂载iso文件：
mkdir /media/matlab - 创建一个文件夹，作为matlab的安装位置：
mkdir /usr/local/matlab/2021a - 进行挂载：
mount -o loop /opt/matlab/Matlab910R2021a_Lin64.iso /media/matlab
3.2.2. 创建激活配置文件
- 切换到配置文件路径：
cd /usr/local/matlab/2021a - 创建配置文件：
touch activate.ini - 添加文件内容：
vim activate.ini
isSilent=true
activateCommand=activateOffline
licenseFile=/opt/matlab/Crack/license_standalone.lic
- 使文件生效：
source activate.ini
3.2.3. 开始安装
- 切换到安装文件路径：
cd /media/matlab - 执行命令并等待：
./install -mode silent -fileInstallationKey 09806-07443-53955-64350-21751-41297 -agreeToLicense yes -licensePath /opt/matlab/Crack/license_standalone.lic -destinationFolder /usr/local/matlab/2021a -activationPropertiesFile /usr/local/matlab/2021a/activate.ini
3.2.4. 破解并激活
- 复制破解文件：
cp /opt/matlab/Crack/libmwlmgrimpl.so /usr/local/matlab/2021a/bin/glnxa64/matlab_startup_plugins/lmgrimpl、cp /opt/matlab/Crack/license.lic /usr/local/matlab/2021a/licenses - 运行激活命令：
/usr/local/matlab/2021a/bin/activate_matlab.sh -propertiesFile /usr/local/matlab/2021a/activate.ini
3.2.5. 取消iso挂载
umount -l /media/matlab
3.2.6. 将matlab添加到环境变量
- vim ~/.bashrc，文件最后添加下面内容：
MATLAB_HOME=/usr/local/matlab/2021a
export PATH=$PATH:$MATLAB_HOME/bin
- 使配置生效：
source ~/.bashrc
3.2.7. 命令测试
matlab -h
Usage: matlab [-h|-help] | [-n | -e]
[v=variant]
[-c licensefile] [-display Xdisplay | -nodisplay]
[--noFigureWindows]
[-nosplash] [-debug]
[-softwareopengl | -nosoftwareopengl]
[-desktop | -nodesktop | -nojvm]
[-batch MATLAB_command | -r MATLAB_command]
[-sd folder | -useStartupFolderPref]
[-logfile log]
[-singleCompThread]
[-jdb [port]]
[-Ddebugger [options]]
[-nouserjavapath]
-h|-help - Display arguments.
-n - Display final environment variables,
arguments, and other diagnostic
information. MATLAB is not run.
-e - Display ALL the environment variables and
their values to standard output. MATLAB
is not run. If the exit status is not
0 on return then the variables and values
may not be correct.
v=variant - Start the version of MATLAB found
in bin/glnxa64/variant instead of bin/glnxa64.
-c licensefile - Set location of the license file that MATLAB
should use. It can have the form port@host or
be a colon separated list of license files.
The LM_LICENSE_FILE and MLM_LICENSE_FILE
environment variables will be ignored.
-display Xdisplay - Send X commands to X server display, Xdisplay.
Linux only.
-nodisplay - Do not display any X commands. The MATLAB
desktop will not be started. However, unless
-nojvm is also provided the Java virtual machine
will be started.
-noFigureWindows - Disables the display of figure windows in MATLAB.
-nosplash - Do not display the splash screen during startup.
-softwareopengl - Force MATLAB to start with software OpenGL
libraries. Not available on macOS.
-nosoftwareopengl - Disable auto-selection of software OpenGL
when a graphics driver with known issues is detected.
Not available on macOS.
-debug - Provide debugging information especially for X
based problems. Linux only.
-desktop - Allow the MATLAB desktop to be started by a
process without a controlling terminal. This is
usually a required command line argument when
attempting to start MATLAB from a window manager
menu or desktop icon.
-nodesktop - Do not start the MATLAB desktop. Use the current
terminal for commands. The Java virtual machine
will be started.
-singleCompThread - Limit MATLAB to a single computational thread.
By default, MATLAB makes use of the multithreading
capabilities of the computer on which it is running.
-nojvm - Shut off all Java support by not starting the
Java virtual machine. In particular the MATLAB
desktop will not be started.
-jdb [port] - Enable remote Java debugging on port (default 4444)
-batch MATLAB_command - Start MATLAB and execute the MATLAB command(s) with no desktop
and certain interactive capabilities disabled. Terminates
upon successful completion of the command and returns exit
code 0. Upon failure, MATLAB terminates with a non-zero exit.
Cannot be combined with -r.
-r MATLAB_command - Start MATLAB and execute the MATLAB_command.
Cannot be combined with -batch.
-sd folder - Set the MATLAB startup folder to folder, specified as a string.
Cannot be combined with -useStartupFolderPref.
-useStartupFolderPref - Set the MATLAB startup folder to the value
specified by the Initial working folder option
in the General Preferences panel.
Cannot be combined with -sd.
-logfile log - Make a copy of any output to the command window
in file log. This includes all crash reports.
-Ddebugger [options] - Start debugger to debug MATLAB.
-nouserjavapath - Ignore custom javaclasspath.txt and javalibrarypath.txt files.
matlab -display Xdisplay
MATLAB is selecting SOFTWARE OPENGL rendering.
< M A T L A B (R) >
Copyright 1984-2021 The MathWorks, Inc.
R2021a (9.10.0.1602886) 64-bit (glnxa64)
February 17, 2021
To get started, type doc.
For product information, visit www.mathworks.com.
>>
3.3 StaMPS
StaMPS is a software package that allows to extract ground displacements from time series of synthetic aperture radar (SAR) acquisitions. The package incorporates persistent scatterer and small baseline methods plus an option to combine both approaches. It is compatible with the TRAIN software and therefore allows to incorporate various tropospheric correction methods in the processing workflow.
StaMPS（Stanford Method for Persistent Scatterers）软件可以从SAR时序数据（ISCE, SNAP, GAMMA, and ROI PAC and DORIS等软件通过对SLC数据预处理可得到）中提取地面位移。
安装步骤：
系统环境：Ubuntu 20.04
安装包下载：https://github.com/dbekaert/StaMPS/releases/tag/v4.1-beta
参考文档：https://homepages.see.leeds.ac.uk/~earahoo/stamps/StaMPS_Manual_v4.1b1.pdf
# 解压后配置环境变量（可以下载已经编译后的压缩文件）
vi ~/.bashrc
# 在文件的最后一行添加一行:export PATH=$PATH:bin目录路径,例如
export PATH=$PATH:/data/apps/StaMPS-4.1-beta/bin
# 使配置生效
source ~/.bashrc
# 将StaMPS下的matlab文件夹添加到MATLAB搜索路径，首先进入MATLAB命令行
matlab -display Xdisplay
% 将文件夹及其子文件夹添加到搜索路径
>> addpath(genpath('/data/apps/StaMPS-4.1-beta'))
% 测试是否配置成功
>> stamps --version
另外需要安装两个依赖软件，包括Snaphu（用来进行解缠）和TRAIN（用来进行大气噪声校正），安装流程与上面一样。
Snaphu：https://web.stanford.edu/group/radar/softwareandlinks/sw/snaphu/
TRAIN：http://github.com/dbekaert/TRAIN
4. 数据处理
4.1 预处理（Pre-processing）
预处理过程基于SNAP软件，具体流程如下：
4.1.1. Sentinel-1 TOPSAR Split
The TOPSAR Split operator provides a convenient way to split each subswath with selected bursts into a separate product. This operator is the first processing step in the TOPS InSAR processing chain.
The user may select the desired subswath with desired bursts and polarisations.
4.1.2. Apply Orbit File
The orbit state vectors provided in the metadata of a SAR product are generally not accurate and can be refined with the precise orbit files which are available days-to-weeks after the generation of the product.
The orbit file provides accurate satellite position and velocity information. Based on this information, the orbit state vectors in the abstract metadata of the product are updated.
4.1.3. InSAR Stack Overview
This function gives a general information about the interferometric stack. The information about the acquisition date, sensor, mode, as well as information about perpendicular and temporal baselines are being listed. Also an estimate for the modeled (expected) coherence is being computed, and used in selection of the optimal reference image for the InSAR stack.
The reference image is selected such that the dispersion of the perpendicular baseline is as low as possible. The reference image is selected maximizing the (expected) stack coherence of the interferometric stack. The “optimal” reference implies improved visual interpretation of the interferograms and aids quality assessment.
4.1.4. Sentinel-1 Back Geocoding
This operator co-registers two S-1 SLC split products (reference and secondary) of the same sub-swath using the orbits of the two products and a Digital Elevation Model (DEM).
In resampling the secondary images into reference frame, deramp and demodulation are performed first to the secondary image, then the truncated-sinc interpolation is performed. Finally, the reramp and remodulation are applied to the interpolated secondary image.
4.1.5. Sentinel-1 TOPSAR Deburst and Merge
For the TOPSAR IW and EW SLC products, each product consists of one image per swath per polarization. IW products have 3 swaths and EW have 5 swaths. Each sub-swath image consists of a series of bursts, where each burst was processed as a separate SLC image. The individually focused complex burst images are included, in azimuth-time order, into a single subswath image, with black-fill demarcation in between, similar to the ENVISAT ASAR Wide ScanSAR SLC products.
For IW, a focused burst has a duration of 2.75 sec and a burst overlap of ~50-100 samples. For EW, a focused burst has a duration of 3.19 sec. Overlap increases in range within a sub- swath.
Images for all bursts in all sub-swaths of an IW SLC product are re-sampled to a common pixel spacing grid in range and azimuth. Burst synchronisation is ensured for both IW and EW products.
Unlike ASAR WSS which contains large overlap between beams, for S-1 TOPSAR, the imaged ground area of adjacent bursts will only marginally overlap in azimuth just enough to provide contiguous coverage of the ground. This is due to the one natural azimuth look inherent in the data.
For GRD products, the bursts are concatenated and sub-swaths are merged to form one image. Bursts overlap minimally in azimuth and sub-swaths overlap minimally in range. Bursts for all beams have been resampled to a common grid during azimuth post-processing.
In the range direction, for each line in all sub-swaths with the same time tag, merge adjacent sub-swaths. For the overlapping region in range, merging is done midway between subswaths.
In the azimuth direction, bursts are merged according to their zero Doppler time. Note that the black-fill demarcation is not distinctly zero at the end or start of the burst. Due to resampling, the data fades into zero and out. The merge time is determined by the average of the last line of the first burst and the first line of the next burst. For each range cell, the merging time is quantised to the nearest output azimuth cell to eliminate any fading to zero data.
4.1.6. Interferogram formation (InSAR operator)
This operator computes (complex) interferogram, with or without subtraction of the flat-earth (reference) phase. The reference phase is subtracted using a 2d-polynomial that is also estimated in this operator.
If the orbits for interferometric pair are known, the flat-earth phase is estimated using the orbital and metadata information and subtracted from the complex interferogram. The flat-earth phase is the phase present in the interferometric signal due to the curvature of the reference surface. The geometric reference system of the reference surface is defined by the reference system of satellite orbits (for now only WGS84 supported, which the reference system used by all space-borne SAR systems).
The flat-earth phase is computed in a number of points distributed over the total image, after which a 2d-polynomial is estimated (using least squares) fitting these ‘observations’, (e.g. plane can be fitted by setting the degree to 1.)
A polynomial of degree 5 normally is sufficient to model the reference phase for a full SAR scene (approx 100x100km). While, a lower degree might be selected for smaller images, and higher degree for ‘long-swath’ scenes. Note that the higher order terms of the flat-earth polynomial are usually small, because the polynomial describes a smooth, long wave body (ellipsoid). To recommended polynomial degree, that should ensure the smooth surface for most image sizes and areas of the world is 5th degree.
In order to reduce the noise, as the post-processing step, you can perform multilooking (with Multilook Operator). Multilooking has to be performed separately on ‘virtual’ bands phase or intensity. In future releases complex Multilook operator will be released. Note that in case of ESA’s ERS and Envisat sensors, the factor 5:1 (azimuth:range) or similar ratio between the factors is chosen to obtain approximately square pixels (20x20 m^2 for factors 5 and 1). Of course the resolution decreases if multilooking is applied
4.1.7. StaMPS Export
Use the StaMPS Export to produce data that can be used within the StaMPS applicaiton for Persistent Scattering Interferometry (PSI).
另外可参考snap2stamps自动进行链式处理。
4.2 永久散射体处理（PS Processing）
Linux命令行执行命令：mt_prep_snap 20211022 /home/data/mexico/mexico_ps 0.4 3 2 50 200
0.4 = amplitude dispersion (0.4-0.42 are reasonable values)
3 = number of patches in range (default 1)
2 = number of patches in azimuth, (default 1)
50 = overlapping pixels between patches in range (default 50)
200 = overlapping pixels between patches in azimuth (default 200)
The number of patches you choose will depend on the size of your area and the memory on your
computer. Generally, patches containing < 5 million SLC pixels are OK.
Linux命令行调用matlab脚本并设置为后台运行：
nohup matlab -nodisplay -batch runsteps >running.log &
扫一扫
4318
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
