| SENTINEL-1 TOOLBOX软件在InSAR数据处理中的应用 |
2. 重庆市智能感知大数据产业技术协同创新中心,重庆,401121
2. 2 Chongqing Collaborative Innovation Center of Smart Sensing & Big Data, Chongqing 401121, China
合成孔径雷达干涉测量(interferometic synthetic aperture radar, InSAR)是一种快速发展的大地测量技术, 能够全天候获取高精度、连续覆盖的地面高程和地表形变信息[1], 已成为对地观测的重要前沿技术之一, 相关的数据处理工作研究与软件开发业得到了各国的重视。SENTINEL-1 TOOLBOX(S1TBX)是由欧空局开发的一款界面友好的开源SAR图像处理软件[2], 它能够处理1级以及更高级的SAR(synthetic aperture radar)数据。S1TBX软件的主要功能是读取由ESA或者其他机构提供的SAR数据, 并实现校准、正射纠正、图像融合、干涉处理、影像重投影、滤波以及其他数据处理[3], 还能将数据转换成通用格式供第三方SAR处理软件使用。S1TBX支持的数据来源包括ALOS PALSAR、ENVISAT ASAR、ERS-1/2、Radarsat-1/2、JERS SAR、TerraSAR-X、Cosmo skymed等SAR卫星数据。本文对采用S1TBX软件进行InSAR数据处理的具体流程及相关原理进行了详细说明。
1 S1TBX处理InSAR数据的具体流程在利用SAR数据进行配准及生成干涉图之前, 需要对影像进行预处理, 以消除或减弱噪声, 改善影像质量, 如应用精密轨道文件替代原始数据中的卫星轨道参数, 以提高配准和地理编码的精度, 并为去除平地效应提供精密轨道数据; 对影像进行辐射校正, 减少辐射偏差等。预处理完成后, 就可以进行影像配准, 生成干涉图以及相干系数图等操作[4], 具体流程如图 1所示。
![]() |
| 图 1 S1TBX数据处理流程 Figure 1 Flow Chart of S1TBX Data Processing |
1.1 更新轨道信息
由于SAR影像自带的卫星轨道信息文件一般不够精确, 可以用精密轨道文件进行优化, 以提高影像配准及地理编码的精度[5]。精轨文件提供了卫星在获取影像时的精确位置和速度信息, 一般在影像获取几天之后就能得到。S1TBX软件能够自动下载的精轨文件有:CTDP(Centre de Traitement Doris Poseidon)提供的ENVISAT精轨文件, DEOS(Delft Institute for Earth-Oriented Space Research)提供的ENVISAT及ERS-1/2精轨文件, 以及Delft大学提供的ERS-1/2精轨文件。应用精轨文件之后, 影像元数据中的轨道状态向量信息将更新。
1.2 辐射定标在对同一地区进行多时相、多角度观测时, 由于SAR轨道和观测角度的变化导致成像几何关系不同, 导致同一地面观测单元的后向散射观测值发生变化, 这种后向散射变化需要采用严格的方法加以校正或补偿, 以产生独立于地形变化的只包含地物表面散射特征信息的SAR数据产品, 即辐射定标[6]。针对不同成像模式的影像数据, 定标采用的参数及公式各不相同:对于地距影像, 需要入射角和定标常数, 而对于斜距复数影像, 需要入射角、定标常数、天线增益、距离向传播损失系数等参数。本文采用的ERS数据的绝对定标公式为:
| $ \sigma _{ij}^0 = \frac{{{\rm{DN}}_{ij}^2}}{K} \times \frac{{{\rm{sin}}{\alpha _{ij}}}}{{{\rm{sin}}{\alpha _{{\rm{ref}}}}}} $ | (1) |
式中, σij0为第i行、第j列像元的后向散射系数; DNij为第i行、第j列像元的原始散射强度; K为绝对定标系数; αij是第i行、第j列像元的雷达波入射角; αref为参考入射角。
1.3 影像配准通过影像配准, 可以将几幅不同的影像转换到同一个坐标系下, S1TBX通过空间配准进行影像的平移、旋转和尺寸改正以及辅影像重采样实现两幅或多幅影像的精确配准。选择一幅影像作为主影像, 其他作为辅影像, 设置好地面控制点的个数、搜索窗口的大小以及插值方法等参数, 即可实现主辅影像的自动配准[7], 对实数影像, 其配准精度可优于0.2个像元; 对于复数影像, 可优于0.05个像元。
1.4 生成干涉纹图和相干系数图将影像配准之后, 就可以将复数影像进行共轭相乘生成干涉相位图以供相位解缠, 同时生成干涉强度图和相干系数图, 以便对相位数据的质量进行评价。
设两幅复数影像分别表示为:
| $ {S_1} = {A_1}{{\rm{e}}^{{\rm{j}}{\varphi _1}}};{S_2} = {A_2}{{\rm{e}}^{{\rm{j}}{\varphi _2}}} $ | (2) |
式中, A1、A2分别表示两幅影像的幅度; φ1、φ2分别表示两幅影像的相位。则干涉图像可表示为:
| $ \begin{array}{l} S = {S_1} \times S_2^* = {A_1}{{\rm{e}}^{{\rm{j}}{\varphi _1}}} \times {({A_2}{{\rm{e}}^{{\rm{j}}{\varphi _2}}})^*} = \\ \;\;\;\;\;\;\;\;\;{\rm{ }}{A_1}{A_2}{{\rm{e}}^{{\rm{j}}\left[ {{\varphi _1} - {\varphi _2}} \right]}} = A{{\rm{e}}^{{\rm{j}}\varphi }} \end{array} $ | (3) |
式中, *表示两幅复数影像进行共轭相乘, 即幅度相乘, 相位相减。
由式(3)可知, 干涉图像的幅度A等于两幅影像的幅度之积; 相位φ等于两幅影像的相位之差。
S1TBX在生成干涉纹图的过程中, 采用基于精密轨道数据去平地效应的方法自动计算并去除了参考相位, 其基本原理为:选择椭球基准面, 在同一参照坐标系(如WGS84坐标系)下, 确定基准面中任一点的坐标和该点所对应的两次成像卫星的位置, 就可以直接计算场景内基准面上的点到卫星的距离差, 对多个点的计算结果进行最小二乘拟合, 就可以得到基准面的参考相位, 从干涉相位中减去参考相位就可以得到由地形高差以及形变引起的相位图。
1.5 从干涉纹图中分离出地形相位干涉相位可以分为参考相位、地形相位、形变相位和噪声等4部分。S1TBX软件在生成干涉纹图的过程中已经去除了参考相位, 为了进行形变分析或其他研究, 可以将地形相位从干涉纹图中分离出来。其基本原理是:将观测区域的外部DEM(digital elevation model)转换到干涉影像对应的雷达坐标系中, 然后采用最邻近插值法对转换后的DEM进行赋值并转换成相位值, 最后从干涉相位中减去由DEM计算出的相位值就能得到去除地形相位的干涉相位。S1TBX软件能够自动下载的外部DEM有SRTM3、ASTER GDEM、GETASSE30等。
1.6 干涉纹图去噪由于系统误差、影像配准误差以及基线去相干等因素的影响, 干涉纹图中往往存在着斑点噪声, 会造成相位数据的不连续和不一致, 严重影响相位解缠的效果。因此尽可能地降低噪声, 提高信噪比, 对提高相位解缠的精度和效率有着十分重要的意义[8]。S1TBX软件提供了两种去除斑点噪声的方式:①通过多视平滑处理去除斑点噪声; ②采用空间域滤波算法去除斑点噪声。S1TBX提供了均值滤波、中值滤波[9]、Frost滤波[10]、Gamma滤波[11]、Lee滤波等滤波算法, 可根据影像的特点进行选择。
1.7 地形改正由于合成孔径雷达是侧视成像, 而且一景影像中地表形态不断变化, 因此实际地物在雷达影像中会发生一系列畸变, 如透视收缩、叠掩、底顶位移等, 因此需要进行相应的改正, 使影像能真实地反映地表情况, 并将其从距离/多普勒坐标系中投影到通用的坐标系中[12]。S1TBX根据卫星轨道状态向量信息、雷达时间参数、斜距改正参数以及参考DEM来获取精确的几何定位信息, 从而对影像进行几何校正, 并投影到WGS84坐标系中。
2 实例分析以意大利西西里岛的埃特纳(Etna)火山地区的两景ERS-1/2影像为例, 对S1TBX软件处理InSAR数据的流程和功能进行实证研究。
2.1 数据预处理首先, 对该区域的两景ERS-1/2影像应用Delft大学提供的ERS-1/2精轨文件更新轨道信息; 然后根据影像元数据中提供的入射角、定标常数、天线增益、距离向传播损失系数对影像进行辐射校正。经过数据预处理后的两幅影像强度图如图 2(a)所示。
![]() |
| 图 2 经过数据预处理后和配准后的主、辅影像强度图 Figure 2 Intensity Graphs of Preprocessed and Coregistered Master, Slave Images |
2.2 影像配准
在影像对中选取200个地面控制点, 以64像素×64像素的搜索窗进行粗配准, 然后缩小窗口, 以32像素×32像素的搜索窗进行精配准。配准后的主、辅影像强度图如图 2(b)所示。
2.3 生成干涉纹图和相干系数图利用S1TBX生成干涉纹图的过程中, 软件自动计算并去除了由平地效应引起的参考相位, 在计算参考相位时, 选择了501个参考相位计算点, 采用5次多项式进行最小二乘拟合。去除参考相位后的干涉纹图如图 3(a)所示。在生成干涉相位图后, 为了对相位数据进行评价分析, 需要生成相干系数图, 如图 3(b)所示。通过比较图 3(a)和3(b)可知, 干涉条纹较为密集的区域对应的相干系数图亮度较暗, 即相干系数较小, 信噪比低, 说明此区域地形复杂, 受噪声干扰较严重; 反之, 在干涉条纹比较稀疏的区域, 对应的相干系数图亮度较亮, 说明此区域地形较为平坦, 相位数据质量好。
![]() |
| 图 3 去除参考相位后的干涉纹图和相干系数图以及转换为分贝值的干涉强度图 Figure 3 Interferogram Subtracted Reference Phase, Coherence Coefficient Diagram and Interference Intensity Diagram Converted to dB |
由于影像配准之后强度信息减弱, 在图像显示时较暗, 因此一般将其转换为分贝值, 以提高亮度, 便于判别分析, 如图 3(c)所示。
S1TBX软件中, 强度图可以分为常规值(linear)和分贝值(dB)两种表示方式, 其转换关系为dB=10×lg(linear)。
2.4 从干涉纹图中分离出地形相位为了从干涉纹图中分离出地形相位, 需要提供外部DEM计算地形相位, 本文采用SRTM(shuttle radar topography mission)提供的DEM计算地形相位, 如图 4(a)所示。从干涉相位中减去地形相位就能得到去除地形相位的干涉相位, 如图 4(b)所示。
![]() |
| 图 4 从干涉相位图中分离出的地形相位图和去掉地形相位后的干涉相位图 Figure 4 Topography Phase Diagram Separated from Interferogram and Interferogram Subtracted Topography Phase |
2.5 干涉纹图去噪以及地形改正
采用多视平滑处理的方法去除斑点噪声, 进行多视平滑处理后的干涉强度图和相位图, 如图 5(a)所示。在进行地形改正时, 根据影像元数据中提供的卫星轨道状态向量信息、雷达时间参数、斜距改正参数以及SRTE3提供的外部DEM进行几何校正并投影到WGS84坐标系中, 经过地形改正后的干涉强度图和干涉相位图如图 5(b)所示。
![]() |
| 图 5 多视处理后以及地形改正后的干涉强度图(分贝值)和干涉相位图 Figure 5 Interference Intensity Images (dB) and Interference Phase Images After Multilooking Processing and Terrain Correction |
为了验证经过几何校正以及投影变换的影像与实际地表的符合程度, 可以将干涉相位图及强度图叠加到Google Earth上进行比较分析, 如图 6所示。从图 6可以看出, 经过地形改正后的影像能够很好地与实际地貌吻合, 表明经过S1TBX软件处理后的影像能较好地反映实际地貌, 能够做进一步的数据处理, 如相位解缠、生成DEM以及形变监测等。
![]() |
| 图 6 在Google Earth上叠加干涉相位图和干涉强度图 Figure 6 Overlayed Interference Intensity Image and Interference Phase Image on Google Earth |
3 结束语
本文详细说明了S1TBX软件在InSAR数据处理中的具体流程以及基本原理, 以Etna火山地区的两景ERS-1/2影像为例进行实证研究, 证明了S1TBX软件能够较好地处理1级以及更高的SAR数据。
| [1] | 许才军, 何平, 温扬茂, 等. InSAR技术及应用研究进展[J]. 测绘地理信息, 2015, 40(2): 1–9 |
| [2] | 李少青, 杨武年, 杨彦通, 等. Sentinel-1A_IW雷达影像在地震形变信息提取中的应用[J]. 测绘, 2016, 39(2): 56–59 |
| [3] |
韩桂红. 干旱区盐渍地极化雷达土壤水分反演研究[D]. 乌鲁木齐: 新疆大学, 2013 |
| [4] | 杜顺季, 邢诚. 基于NEST的差分干涉测量及其应用[J]. 地理空间信息, 2013, 11(3): 70–72 |
| [5] | 陈强, 刘国祥, 李永树. 粗/精轨道数据对卫星InSARDEM精度影响的对比分析[J]. 遥感学报, 2006, 10(4): 475–481 DOI: 10.11834/jrs.20060471 |
| [6] | 陈金星, 张波, 王超, 等. 高分辨率SAR辐射定标中两种角反射器响应测量方法的应用对比研究[J]. 遥感技术与应用, 2015, 30(4): 677–683 |
| [7] | 喻小东, 郭际明, 黄长军, 等. 基于SIFT算法的InSAR影像配准方法试验研究[J]. 遥感信息, 2013, 28(2): 66–69 |
| [8] | 朱邦彦, 李建成, 储征伟, 等. 利用时序InSAR反演常州市地表沉降速率[J]. 测绘通报, 2016, (5): 26–31 |
| [9] | 季伟, 戴鑫. InSAR干涉条纹图滤波方法比较与分析—以南屯矿区为例[J]. 测绘地理信息, 2014, 39(5): 14–17 |
| [10] | 于双, 徐佳, 杨英宝, 等. 高分辨率SAR图像杂波模型对比研究[J]. 测绘科学, 2014, 39(1): 133–136 |
| [11] | 朱磊, 韩天琪, 水鹏朗, 等. 一种抑制合成孔径雷达图像相干斑的各向异性扩散滤波方法[J]. 物理学报, 2014, 63(17): 453–463 |
| [12] | 王超, 张红, 刘智. 星载合成孔径雷达干涉测量[M]. 北京: 科学出版社, 2002 |
2018, Vol. 43







