王洪明 李如仁 陈振炜 张过

(1 沈阳建筑大学 交通工程学院，沈阳 110168)(2 长春建筑学院 交通学院，长春 130604)(3 武汉大学 测绘遥感信息工程国家重点实验室，武汉 430079)

随着国家高分辨率卫星对地观测计划实施完成，高分三号(GF-3)作为“高分”家族中唯一的合成孔径雷达(SAR)卫星，具有高分辨率、大成像幅宽、多成像模式等特点，与欧洲哨兵-1(Sentinel-1)卫星等的免费数据相比具有高分辨率的优势，与陆地合成孔径雷达(TerraSAR)等卫星数据相比具有更高的性价比[1]。因此，分析GF-3卫星在地表形变监测与地质灾害防治领域的应用能力具有深远意义[2]。SAR差分干涉测量(D-InSAR)技术已广泛应用于测绘地表变形，如火山学、城市和基础设施变形分析[3]，以及地质灾害早期识别[4]。其原理是：通过对2幅SAR图像在不同时间的雷达相位进行差分[5]，得到2幅SAR图像在视线方向(LOS)上发生的地表变形。

目前，利用GF-3卫星进行的测量研究如下。文献[1]利用GF-3卫星影像，基于雷达影像强度的偏移量追踪技术，提取了贡嘎冰川表面的运动情况，填补了利用GF-3卫星影像监测冰川运动的空白。文献[6]对GF-3卫星干涉测量进行研究，提出基于航天飞机雷达地形任务(SRTM)数字高程模型(DEM)上控制点信息消除轨道误差的方法，提高了形变提取精度。文献[7]首次利用GF-3卫星数据，基于时序InSAR技术完成了对北京地表的沉降量提取。其中：文献[6]中仅考虑消除轨道误差对测量精度的影响，没有考虑InSAR处理过程中的图像配准和图像滤波过程对干涉测量精度的影响。文献[7]中仅从时序InSAR处理全过程上完成形变提取，没有考虑形变图像配准精度和滤波方法对形变相位保持能力的问题。另外，领域内对GF-3卫星存在干涉测量应用较少和精度较低等问题。为提升GF-3卫星在地表形变提取方面的能力，提高卫星在国家防灾减灾应用的效能，本文在GF-3卫星D-InSAR的处理过程中，提出“基于影像几何与内容特征的配准”和“多尺度相位滤波法”2种方法，利用2景GF-3卫星单视复数影像，基于雷达差分干涉测量技术提取了研究区域地表形变，结合外部数据对比分析，评估地表形变提取的精度，提高了GF-3卫星在地表形变监测应用中的效率与精度，可为今后GF-3卫星用于地表形变监测提供理论依据，促进其在地质灾害早期识别与防治领域中的应用。

1 研究数据与方法

1.1 数据

GF-3卫星的复原轨道精度优于10 m(1σ)，精密轨道精度优于20 cm(1σ)[8]。轨道的不精确性导致较大的轨道误差，进而导致基线参数信息不够精确。因此，轨道误差的处理和消除尤为重要。此外对于平地相位测量，地形相位的误差和地形相位的不连续性，也限制着合成孔径雷达干涉测量InSAR形变提取。经过影像筛选，本文研究数据为2景GF-3卫星单视复数(SLC)影像，其具体参数如表1所示，时间基线为29天，空间垂直基线为1 106.913 430 3 m。辅助InSAR处理的数据为30 m空间分辨率的Versions 4 SRTM DEM数据，用其进行地形相位去除。

表1 GF-3卫星影像参数

1.2 研究方法

本文采用D-InSAR方法对2景雷达单视复数影像进行处理。差分干涉测量方法要求2景SAR影像在相同的地理坐标系下具备相同的入射角，并且成像区高度重叠。在2次不同的回波信号中记录着条纹在距离向SAR天线和地面点的变化轨迹，除去地形信息，可得到差分干涉相位值，这个值记录着形变的信息。差分干涉测量技术以雷达干涉测量为理论基础，相位可用式(1)表示。

φint=φfla+φter+φdef

(1)

式中：φint为干涉相位；φfla为平地相位；φter为地形相位；φdef为形变相位。

将φint中φfla平地相位与φter地形相位减去之后即可得到形变相位φdef，再由相位转形变公式即可得到形变量Δd，如式(2)所示。

(2)

式中：λ为雷达波长。

需要注意的是，在二轨法处理形变时应尽量挑选空间基线较小的影像对。而且，由式(1)和式(2)可知：长波段的SAR数据(例如日本先进陆地观测卫星(ALOS)数据)，其干涉条纹稀疏。

在差分干涉测量方法中，二轨法和三轨法最为常用，二轨法基于外部DEM去除地形相位，三轨法基于2次干涉，2张干涉图相位相减去除地形相位。本文采用二轨法，基于SRTM数据，在GAMMA软件支持下完成地表形变提取。主要实现过程如图1所示。单视复数(SLC)图像配准、生成干涉图、干涉图滤波、地形相位模拟、生成差分干涉图、相位解缠、形变提取与地理编码。单视复数SAR图像配准是基于影像几何与内容特征的配准，过程为基于轨道参数的影像粗配准，采用5×5多视平均SAR影像的距离向分辨率和方位向分辨率；影像精配准从主影像和辅影像的子窗体(256像素×256像素)中，自动计算同名地物像素之间的交叉相关函数，进而保证方位向和距离向上配准精度达到1/10像素；基于拟合函数将辅助影像重采样到主影像下。干涉图生成与相干性计算，即重采样后的辅影像与主影像共轭相乘，得到干涉图与相干图。干涉图滤波，是基于多尺度滤波法对干涉图进行滤波，减少干涉条纹中由时空基线引起的噪声。基于DEM的地形相位模拟与地形相位去除，就是将DEM文件由小端模式改为大端模式，利用查找表精确地模拟地形相位，进而将其去除。相位解缠，即采用最小费用流解缠方法，为了尽量保证相位解缠的连续性，不对差分干涉图分块处理，利用狄洛尼(Delaunay)三角剖面，进行相位解缠，提高研究区干涉条纹不连续区域的相位解缠精度，设定的相位解缠阈值范围是0.10～0.15。

图1 本文D-InSAR实现过程

2 试验与分析

2.1 形变提取与分析

波长、入射角和像元分辨率都影响着形变的提取。C波段的GF-3卫星对地表微小形变的灵敏性要优于其他长波段雷达卫星；其相对较小影像的入射角间接增加了差分干涉图中的条纹数；较高的像元分辨率提高了其最大形变梯度值。

与传统方法不同的是，本文采用大小为256×256的配准窗口，为保证配准窗口不重叠地覆盖整景影像，距离向和方位向上的窗口数目设置为64×48，一共需要计算3068个偏移量。本文提出的雷达影像配准方法为：首先，利用轨道数据对影像进行配准；然后，通过初步配准后的主辅影像的信杂比(SCR)阈值，结合阈值再选取一系列配准点，得到3852个配准点，即偏移量的计算量比传统方法高约25%。传统方法的二项式拟合均方根误差在距离向上为0.064，方位向上高达0.286，而本文配准方法在2个方向上则分别是0.027和0.082，明显优于传统方法。图2表示2种方法的多项式差异，方位向上的最大差异超过了0.35像素。

图2 本文方法与传统方法影像方位向偏移量差异

对于滤波，本文提出的多尺度相位滤波方法基于2个尺度的，对干涉图进行22和24这2个尺度进行分割。首先，将干涉图分解成实部信号与虚部信号，两者在自适应阈值情况下进行信号掩膜生长，其中，自适应阈值S2是指动态的目标分块所有均值，具体筛选像素信号由干涉像素的质量S1和阈值S2进行对比，当S1大于S2时为信号，否则为噪声。其次，将掩膜之后的干涉图信号部分的系数实部、虚部都乘以2，噪声系数不变。然后，将分离出来的信号通过变换处理计算干涉图的干涉相位值。最后，按照上述3步进行2次即可。与传统滤波方法的仿真对比，如图3所示。评价指标为目视效果判读与滤波前后相位信噪比，其信噪比公式为

图3 本文滤波方法与传统滤波方法对比

(3)

式中：σM为滤波前后均方差。

本文方法的信噪比为42.11 dB，明显高于传统滤波方法的39.63 dB，也代表着多尺度滤波有更好的去噪声和保留形变相位信息的能力。

基于SRTM进行反向地理编码模拟的SAR影像，如图4所示。模拟SAR影像与主影像在方位向和距离向均达到0.01像素，进而确定精确查找表，得到地形相位图如图5所示，保证地形相位的精确去除。

图4 模拟SAR影像

图5 模拟地形相位图

基于本文滤波方法处理后的差分干涉图，如图6所示。接下来，在相干图中设置0.3的相位解缠掩膜参数，将相干性低于0.3视为影像失相干[9]，最大限度保证相位解缠的精度。由于最小费用流(MCF)解缠方法既可减少解缠信息孤岛的出现，又可满足计算效率与精度要求，因此，相位解缠方法选择最小费用流方法[10]。

图6 差分干涉图

相位解缠后的图中的相位信息为形变信息，形变信息代表雷达视线向形变，指向雷达方向为正(如图7所示)，可得到城市周边山区地表存在明显的沉降漏斗，城市地表形变较为均匀，无异常沉降，形变量相对于周围地形是一个相对值，形变量±0.01 m表示1组条纹周期。雷达对地表形变探测的敏感度与波长成反比，GF-3卫星采用C波段SAR，因此相比与其他卫星SAR数据在同样的地区可监测出更多地方明显的形变。在研究区，除了城区外，共发现7处形变区域，1处为明显的形变漏斗。城市周边的形变漏斗区域与国土资源部门相关信息相符，此区域为一个大型矿场，随着不断开采，造成不同地区地表的沉降，这也符合研究结果。

图7 雷达视线向形变图

2.2 误差分析与验证

差分干涉测量主要误差源为SAR影像的配准、地形相位模拟、相位解缠等。其中：基于外部轨道信息配准消除由于轨道不准确引起的配准误差；基于外部DEM数据配准可以保证模拟的地形相位精度较高，也可提升影像的配准精度。本文使用30 m分辨率的SRTM数据辅助SAR影像配准与地形相位模拟，GF-3卫星影像对配准精度最终优于1/10个像元。在像对配准时考虑了像对的几何与内容特征，以增加一定计算量为代价获取较高的配准精度，偏移量的计算量比传统配准方法高约25%，二项式拟合均方根误差明显低于传统配准方法。利用多尺度相位滤波将干涉图进行2个尺度分割，与经典滤波方法相比提高了形变相位保持能力，且信噪比较高。

由于本文采用的GF-3卫星数据时间基线为29天，为了验证监测结果的精度，用研究区域中的7处形变区域与同期的GPS点监测结果进行比较，如表2所示，形变速率最小互差为0.76 mm/a，最大互差为5.64 mm/a，形变提取精度达到毫米，形变提取精度有所提高。

表2 GPS监测结果与同期InSAR监测结果比较

3 结论

本文为提高GF-3卫星干涉测量精度与效率，以2景GF-3卫星单视复数影像作为干涉处理研究对象，得出结论如下。

(1)提出基于影像几何与内容特征的配准，配准精度达到0.1像素，偏移量的计算量比传统配准方法高约25%，利用配准的拟合函数将辅助影像重采样到主影像下，保证图像的干涉质量。为更好去除相位噪声，本文针对影像采用多级滤波方法，与传统滤波方法相比，其信噪比明显提高，可更好地保留出形变相位。

(2)基于MCF相位解缠方法，设置0.3的相干性阈值与0.10～0.15解缠阈值，对差分干涉图进行了相位解缠，最后将相位解缠图转为形变图。最终结果显示：城市地表形变较为均匀，无异常沉降，研究区域除城区外，共发现7处形变区域，1处为明显的形变漏斗。形变量相对于周围地形是一个相对值，形变量为±0.01 m表示一组条纹周期。明显的形变漏斗与国土资源部门的相关数据得到印证，其位置为一处大型矿场，其大规模开采导致地表出现大量不均匀沉降，符合研究结果。此外，通过外部数据对比分析可知：二者形变区域总体变化趋势相同，形变提取精度够高。

上述研究结论表明：本文方法适用于GF-3卫星地表形变提取且精度较高，为GF-3卫星在地表形变监测应用中提供了一种技术手段。值得注意的是，GF-3卫星影像对需要精心挑选，除了在较短的重访周期中其空间基线长度尽量小，还要满足像对入射角大致相等和影像覆盖面积大致相同的要求。不过，随着GF-3卫星组网的完成，其数据积累量日益增加，卫星的重访周期变短，可用的影像对越来越多，其空间基线长度也随之达到干涉测量要求，GF-3卫星将在地表形变监测领域、地质灾害早期识别等领域应用更加广泛。