Co-registration of Image Stacks for Sentinel-1A TOPS Mode Based on Dijkstra's Algorithm
-
摘要: 目前,哨兵卫星被广泛应用于监测地表形变,然而其默认成像模式TOPS(terrain observation with progressive scanning)受粗配准的精度限制,方位向的频谱混叠会导致重叠观测区域出现相位跳变,需要通过精配准进行纠正。多数开源软件都采用几何配准和增强谱分集结合的方式对哨兵影像进行精配准。增强谱分集的精度通常受到时间去相干、几何去相干和方位向形变等诸多因素的影响,但其直接影响因素是相干性,研究了增强谱分集的精度与相干性的关系。同时,为了提高时序影像的配准精度,提出在精确估计相干性的基础上提取高相干点用于增强谱分集,并应用迪杰斯特拉最短路径算法生成最优配准像对,完成时序像对的精配准。利用传统的单主影像增强谱分集和目前精度最高的序贯网络增强谱分集(network-based enhanced spectral diversity,NESD)对所提出的方法进行了精度验证。实验证明,所提方法能够达到对哨兵影像进行有效精配准的目的,并且弥补了NESD方法的不足。Abstract: Today, Sentinel-1A data with terrain observation with progressive scanning (TOPS) imaging mode is increasingly used in earth observation aiming at consistent monitoring of surface change and its deformation. However, due to the limited accuracy of coarse co-registration, spectral aliasing along the azimuth direction enables the presence of phase jumping in overlapping area of neighboring bursts. Although geometrical co-registration in conjunction with enhanced spectral diversity (ESD) is proven to be a feasible strategy to correct such error and has been widely used in some open source softwares (e.g. DORIS (Delft object-oriented radar interferometric software), SNAP (sentinel applications platform), ISCE (InSAR scientific computing environment)), the performance of ESD is still not satisfactory and relies strongly on spatiotemporal decorrelation. Given that decorrelation is generally quantized by interferometric coherence, this paper presents and assesses a new methodology to improve the accuracy of ESD by fully exploring the high coherent targets. Specifically, this method focuses on mitigating the spatiotemporal decorrelation in fine co-registration procedures by:(1) Selecting stable targets with moderate and high coherence using accurate coherence estimation. (2) Maximizing coherence magnitude by optimal interferometric subset chosen from Dijkstra algorithm. We compare and test this method against current favorites based on single master image and network-based enhanced spectral diversity (NESD), and the experimental results demonstrate the value of our method. It can make up for the shortage of NESD method.
-
Keywords:
- Sentinel-1A /
- image co-registration /
- enhanced spectral diversity /
- Dijkstra algorithm /
- TOPS mode
-
2014年4月,欧洲空间局成功发射了哨兵卫星Sentinel-1A[1]。Sentinel-1A凭借其宽幅成像的特点,受到越来越广泛的关注[2]。Sentinel-1A的宽幅成像模式由3个子带组成,每个子带又分为多个雷达波束影像。3个子带在沿袭TerraSAR-X卫星TOPS(terrain observation with progressive scanning)[3]模式的同时,进一步提高了卫星的姿态与轨道的控制精度。Sentinel-1A不仅兼顾了干涉性能,而且能够获得约250 km幅宽的影像,可以达到方位向约20 m、距离向约5 m的影像分辨率。
TOPS成像模式下,哨兵卫星的雷达天线周期性地逆时针旋转,使得重叠区域观测目标的方位向频谱产生混叠。由于影像配准方法的精度限制,辅影像残余的配准偏移作用于观测目标的频谱,会使干涉图出现相位跳变。为了使重叠区域的干涉相位差小于3°(以避免干涉图出现相位跳变),方位向的精配准精度要求达到约0.001像元。
几何配准方法[4]与实数互相关方法[5]被广泛应用于TOPS模式的粗配准,其能够有效地将方位向配准偏移导致的相位跳变控制在[-π, π]之间,从而避免缠绕效应对残余方位向偏移量估计造成影响[6]。但受限于多项式拟合的精度,实数互相关方法通常难以描述地形对辅影像配准偏移产生的影响,对于地形复杂的区域,难以达到期望的精度。几何配准方法的偏移值获取过程与实数互相关方法所采用的多项式拟合不同,该方法利用零多普勒效应对影像所有点进行偏移量估计,因此消除了地形所产生的影响[7];且该方法主要依赖于卫星轨道精度[8],方位向配准偏移误差通常表现为系统误差,因此认为在几何配准后,每组雷达波束影像中的残余配准误差为固定值,易被后续精配准方法所校正。本文采用几何配准方法作为TOPS影像的粗配准方法。
为了精化辅影像所存在的残余配准误差,目前常采用谱分集(spectral diversity,SD)[9]与增强谱分集(enhanced spectral diversity,ESD)进行TOPS模式下方位向残余配准的估计。Prats-Iraola等[10]针对SD与ESD在TerraSAR-X TOPS模式下的表现进行了对比,相对于SD采用的约4 kHz的脉冲采样频率[11],ESD采用约5 kHz的多普勒中心频率差进行方位向配准偏移的估计,精度提高了约1.25倍。针对Sentinel-1A,Yague-Martinez等[12]也对ESD方法在TOPS模式下的方位向配准表现进行了评估,实验证明ESD方法能够有效地提取Sentinel-1A的方位向配准偏移。ESD方法的核心思想为:TOPS模式下,相邻雷达波束影像之间的多普勒质心频率差远远大于方位向带宽,对方位向残余配准误差十分敏感,采用主辅影像的相邻雷达波束影像视为上下频段,将上频段与下频段各自获取的干涉图继续进行干涉来获得因方位向配准偏移的存在而导致的干涉相位差,通过估计出的多普勒质心频率差即可获得固定的方位向残余配准值。
受去相干效应的影响,为了准确获取重叠区域的干涉相位差峰值,多数开源软件(如SNAP(sentinel applications platform)[13]、DORIS(Delft object-oriented radar interferometric software)[14]及ISCE(InSAR scientific computing environment)[15])只采用高于某相干性阈值(如0.15)的像元进行计算。用此方法进行相干性估计时,忽略了窗口中残余条纹对相干性造成的影响,且在存在方位向配准偏移的情况下,直接利用低估的相干性进行像元的选取,容易纳入粗差点,使得偏移计算过程出现离群值。为了避免离群值影响精配准结果,吴文豪等[6]采用粗差探测方法去除ESD离群值,提高了单像对配准精度,但此方法仍未考虑相干性估计精度对ESD造成的影响。时序影像配准过程中,受时间去相干与噪声的影响,低相干性像对的选点精度降低,配准偏移值的估计存在偏差。为了减小相干性对ESD的影响,应避免利用单主影像方式进行影像精配准。Fattahi等[16]、Yague-Martinez等[17]及Xu等[18]采用小基线集式[19]的时序影像配准方式,减小了时间去相干对ESD的影响,也增强了ESD计算偏移值的鲁棒性,但该方法仍存在两个问题:(1)多像对造成了配准效率的下降;(2)干涉对噪声特征各异,导致方位向配准偏移量的估计精度受限于低相干性区域。
为解决上述问题,进一步增强TOPS模式时序合成孔径雷达(synthetic aperture radar,SAR)影像的配准精度和计算效率,本文利用精确估计的相干性提取高相干点并进行雷达波束影像方位向配准偏移的计算,同时提出利用迪杰斯特拉算法[20]寻找最优配准集合,在提高方位向配准精度的基础上进一步提高运算效率。
1 基于ESD的Sentinel⁃1A精配准方法
1.1 ESD原理
目前,Sentinel-1A的主流配准方法分为两个步骤:先将辅影像经过几何配准重采样至主影像(粗配准),再利用ESD方法估计辅影像残余的方位向配准偏移值并进行纠正(精配准)。在进行ESD时,需要获得雷达波束影像重叠区域的干涉相位之差${\rm{\Delta }}{\phi _{{\rm{error}}}}$(如图 1中的虚线标注区域):
$${\rm{\Delta }}{\phi _{{\rm{error}}}} = \angle (\left( {{m_i} \cdot s_i^{\rm{*}}} \right) \cdot \left( {{m_{i + 1}} \cdot s_{i + 1}^{\rm{*}}{)^{\rm{*}}}} \right)$$ (1) 式中,*为共轭算子;$\angle $()为复数取相角操作;${m_i}$与${m_{i + 1}}$分别为主影像相邻的雷达波束影像;${s_i}$与${s_{i + 1}}$分别为辅影像相邻的雷达波束影像。在没有配准偏差和噪声影响的理想状态下,${\rm{\Delta }}{\phi _{{\rm{error}}}}$为0。然而,辅影像经过几何配准后通常包含残余配准偏移,需要结合目标频谱信息来估计残余的偏移像元并进行纠正。
简化TOPS成像模型,假设主辅影像聚焦至相同多普勒中心频率${f_c}$和方位向带宽$B$, 且辅影像已粗配准至主影像,目标的脉冲回波可概括为自变量为方位向时间$t$的函数[4, 21],表达式如下:
$$M\left( t \right) = {{\rm{e}}^{{\rm{j}}{\phi _m}}}{{\rm{e}}^{{\rm{j}}2{\rm{\pi }}{f_c}\left( {t - {t_0}} \right)}}{\rm{sinc}}\left( {{\rm{\pi }}B\left( {t - {t_0}} \right)} \right)$$ (2) $$S\left( t \right) = {{\rm{e}}^{{\rm{j}}{\phi _s}}}{{\rm{e}}^{{\rm{j}}2{\rm{\pi }}{f_c}\left( {t - {t_0}} \right)}}{\rm{sinc}}\left( {{\rm{\pi }}B\left( {t - {t_0}} \right)} \right)$$ (3) 式中,M(t)、S(t)分别为主、辅影像成像目标的雷达脉冲回波函数;sinc( )为辛格函数; ${\phi _m}$、${\phi _s}$分别为主、辅影像成像目标的相位; ${t_0}$为成像目标的零多普勒时间; ${\rm{\Delta }}t = t - {t_0}$为辅影像误配准带来的方位向时间偏移。目标相位差表示为 ${\phi _t} = {\phi _m} - {\phi _s} + {\phi _{{\rm{error}}}}$, ${\phi _{{\rm{error}}}}$即为方位向配准误差带来的相位偏差。由于多普勒中心频率已知,相位偏差可以表示为 ${\phi _{{\rm{error}}}} = 2{\rm{\pi }}{f_c}{\rm{\Delta }}t$,重叠区域的干涉相位差即为 ${\phi _{{\rm{error}}}} = 2{\rm{\pi }}{f_c}{\rm{\Delta }}t$,其中 ${f_{{\rm{ovl}}}}$为多普勒中心频率差。实际上,在TOPS成像过程中,多普勒中心频率并不是固定的,而是以一定的线性速率随方位向时间变化[3, 10, 16],系统性的方位向配准误差 ${\rm{\Delta }}t$带来的相位偏差 ${\rm{\Delta }}{\phi _{{\rm{error}}}}$保持着相同的线性速率变化,表现为方位向多普勒历程斜坡的形式。TOPS模式下,观测目标的多普勒历程如图 2所示。
图 2中,${k_a}$为目标的多普勒调频率,${k_{{\rm{rot}}}}$为多普勒中心调频率,${t_c}$为雷达波束影像1中该目标的方位向时间,${t_0}$为目标的零多普勒时间,目标在相邻雷达波束影像中的方位向时间为${t_c}^{\rm{'}}$,tcycle为多普勒时间历程周期时间。TOPS模式与其他成像模式不同,观测过程中,天线始终以固定速度$\omega $周期性地逆时针扫描。由于天线的旋转与卫星的移动赋予了观测目标复杂的速度分量,导致${k_{{\rm{rot}}}}$以一定的规律不断变化[10]。计算方位向配准误差对应的方位向时间${\rm{\Delta }}t$需计算精确的${f_{{\rm{ovl}}}}$,而从图 2中的几何关系可获得${f_{{\rm{ovl}}}}$。方位向时间${t_c}$可以通过如下关系得到,由
$${k_{{\rm{rot}}}} \cdot {t_c} = {k_a} \cdot \left( {{t_c} - {t_0}} \right)$$ (4) 可得到:
$${t_c} = \frac{{{k_a}{t_0}}}{{{k_a} - {k_{{\rm{rot}}}}}}$$ (5) 其中,
$${k_{{\rm{rot}}}} = \frac{{2V}}{\lambda }\omega $$ (6) 式中,V为卫星飞行切向速度;λ为信号波长。
同样,方位向时间${t_c}^{\rm{'}}$也可以通过相似的几何关系计算得到[10, 16],由
$$\left( {{t_c}^{\rm{'}} - {t_{{\rm{cycle}}}}} \right) \cdot {k_{{\rm{rot}}}} = {k_a} \cdot \left( {{t_c}^{\rm{'}} - {t_0}} \right)$$ (7) 可得到:
$${t_c}^{\rm{'}} = \frac{{{k_a} \cdot {t_0} - {k_{{\rm{rot}}}} \cdot {t_{{\rm{cycle}}}}}}{{{k_a} - {k_{{\rm{rot}}}}}}$$ (8) 由式(5)、式(8)可得到${f_{{\rm{ovl}}}}$:
$${f_{{\rm{ovl}}}} = \left| {{k_a} \cdot \left( {{t_c}^{\rm{'}} - {t_c}} \right)} \right| = \left| {\frac{{{k_a} \cdot {k_{{\rm{rot}}}}}}{{{k_a} - {k_{{\rm{rot}}}}}}} \right| \cdot {t_{{\rm{cycle}}}}$$ (9) 传统的ESD方法是先将雷达波束影像重叠区域等分为多个子块,提取各子块中相干性估值高于经验阈值(如0.15)的像元,然后结合多普勒中心频率差${f_{{\rm{ovl}}}}$与方位向采样间隔时间$\tau $来获得方位向残余配准像元${\rm{\Delta }}x$,最后利用粗差探测方法纠正离群值,从而计算出较为可靠的方位向偏移。${\rm{\Delta }}x$计算如下:
$${\rm{\Delta }}x = \frac{{{\rm{\Delta }}{{\bar \phi }_{{\rm{error}}}}}}{{2{\rm{\pi }}{f_{{\rm{ovl}}}}\tau }}$$ (10) 式中,${\rm{\Delta }}{\bar \phi _{{\rm{error}}}}$为各子块干涉相位之差的平均值。
在校正辅影像方位向偏移的过程中,需要对TOPS影像进行去斜操作,去斜后根据${\rm{\Delta }}x$及头文件中的多普勒调频参数,通过快速傅里叶变换补偿相位部分,从而达到精配准的目的。单像对ESD影像配准偏移量中误差可以表示为[11]:
$${\sigma _{{\rm{\Delta }}x}} = \frac{1}{{2{\rm{\pi }}{f_{{\rm{ovl}}}}\tau }}\sqrt {\sigma _{{\phi _i}}^2 + \sigma _{{\phi _{i + 1}}}^2 - 2{\sigma _{{\phi _{i, i + 1}}}}} $$ (11) 式中,${\sigma _\phi }$通常采用干涉相位差的克拉美罗下界${\sigma _\varphi } = \frac{1}{{\sqrt {2N} }} \cdot \frac{{\sqrt {1 - {\gamma ^2}} }}{\gamma }$来表示[5, 10, 16],其中$\gamma $为计算出的理论相干性;N为计算相干性时窗口中的独立像元个数。由于雷达波束影像重叠区域干涉相位的噪声部分互不相关,忽略相邻雷达波束影像之间的干涉相位协方差${\sigma _{{\phi _{i, i + 1}}}}$,则ESD方位向配准偏移量中误差${\sigma _{{\rm{\Delta }}x}}$可简化为$\frac{{\sqrt 2 {\sigma _\varphi }}}{{2{\rm{\pi }}{f_{{\rm{ovl}}}}\tau }}$。ESD影像配准偏移量中误差与相干性的关系如图 3所示。
1.2 时序影像精配准方法
粗配准后,为了纠正所有辅影像残余的方位向配准偏移,常采用的方法有单主影像增强谱分集(single-master enhanced spectral diversity,SM-ESD)、序贯网络增强谱分集(network-based enhanced spectral diversity,NESD)。SM-ESD即采用单幅影像作为主影像,所有辅影像经方位向配准、偏移量估计并纠正偏移至主影像, 如图 4(a)所示。SM-ESD完成N幅时序影像的精配准只需经过N-1次计算,但随着时间基线的增加,主辅影像间的相干性降低,该方法的精度也在降低。
为了解决时间去相干对时序像对精配准的影响,Fattahi等[16]提出了NESD方法,它是利用小基线集的方式组成配准对并计算出所有组合之间的方位向配准偏移${\rm{\Delta }}x$,通过加权最小二乘的方式获得辅影像与几何配准主影像的方位向偏移$x$, 如图 4(b)所示。为了削弱像对间的时间去相干,NESD方法通常采用时间节点较近的N幅影像作为一个小基线集,其原理描述为[16]:先利用线性方程式$Ax = {\rm{\Delta }}x$描述图论中各影像连接线代表的方位向偏移矩阵${\rm{\Delta }}x$,A矩阵通常称为图论中的入射角设计矩阵,再通过$x = {({A^{\rm{T}}}PA)^{ - 1}}{A^{\rm{T}}}P{\rm{\Delta }}x$获得各辅影像与主影像间的方位向偏移矩阵,P由图 4中各影像连接线之间的权$1/\sigma _{{\rm{\Delta }}{x_{i, j}}}^2$组成,$\sigma _{{\rm{\Delta }}{x_{i, j}}}^2$为式(11)中依据克拉美罗下界[22]计算出的方位向配准偏移方差。
2 利用迪杰斯特拉算法进行ESD精配准
2.1 高相干点ESD
在传统的ESD精配准过程中,相干性未进行相位补偿的情况下,直接选取重叠区域中相干性高于经验阈值(threshold)的点,叠加上述高于相干性阈值点的干涉相位差获得${\rm{\Delta }}{\phi _{{\rm{error}}}}$,即:
$${\rm{\Delta }}{\phi _{{\rm{error}}}} = \angle {(\left( {{m_i} \cdot s_i^{\rm{*}}} \right) \cdot {({m_{i + 1}} \cdot s_{i + 1}^{\rm{*}})^{\rm{*}}})_{\gamma > {\rm{threshold}}}}$$ (12) 获得干涉相位差${\rm{\Delta }}{\phi _{{\rm{error}}}}$后,利用雷达波束影像重叠区域的目标多普勒中心频率差能够计算出方位向配准偏移。这类方法在计算过程中既未考虑相位梯度引起的相干性低估,也未考虑配准误差带来的相干性损失,导致像元的选取精度受到影响。为了提高像元的选取精度,在相干性估计过程中,需要保证滑动窗口中的相位条纹已去除[23-25],因此在相干性估计时需补偿地形[26],即:
$$\gamma = \frac{{\sum\limits_{i = 1}^N {{s_{1i}}s_{2i}^{\rm{*}}{{\rm{e}}^{{\rm{j}}{\phi _{{\rm{topo}}}}}}} }}{{\sqrt {\sum\limits_{i = 1}^N {{{\left| {{s_{1i}}} \right|}^2}} } \sqrt {\sum\limits_{i = 1}^N {{{\left| {{s_{1i}}} \right|}^2}} } }}$$ (13) 式中,N为滑动窗口的大小;${s_{1i}}$和${s_{2i}}$分别代表主、辅影像滑动窗口中第i个像元的复数信号;$i$代表地形相位。
利用上述较为精确的相干性来提取整个重叠区域的高相干点作为相位贡献点,以进行干涉相位差的计算。计算过程未采用分块形式,是因为重叠区域的子块中如果含有水体等严重失相干区域,这时由于缺少分布均匀的相干点,估计出的方位向偏移值与实际偏移值之间存在较大偏差。如果整个重叠区域的多个子块都出现严重失相干情况,则粗差探测方法将无法适用于整体区域的偏移值估计。
利用高相干点ESD和传统ESD方法计算干涉相位差的过程是一致的,但不同的是:(1)在精确估计相干性的基础下,高相干点ESD方法在相位贡献点的选择过程中,采用整个雷达波束影像重叠区域相干性更高的点(本文采用$\gamma > 0.6$),理论上减小了方位向配准偏移中误差;(2)高相干点ESD不采用分块形式,避免了在研究区域富含水体的情况下误估了配准偏移。
2.2 迪杰斯特拉算法
精配准过程中,为选取分布均匀的相位贡献点,传统ESD方法和高相干点ESD方法都依赖于相干性较高的像对,但时序像对中通常会存在相干性较低的主辅影像,如采用SM-ESD方法进行精配准时,时间基线较长的像对很难获得准确的配准偏移值。为了获得分布均匀的高相干点,从而得到较精确的方位向配准偏移值,可采用规划配准像对的方式,利用中间影像进行配准偏移量传递。为了尽可能保证配准偏移量的准确,必须确保传递路径中节点影像间具有较高的相干性,而传递过程中需要解决的是如何量化像对间优劣性的“长度”指标。Refice等[27]利用$1 - \gamma $作为路径的距离,将最小生成树算法在SAR图像配准过程中进行了理论延伸,取得了较好的实验效果。$1 - \gamma $作为路径距离,实际上是为了最大化路径的相干性,但并不能直接说明根据此路径计算出的配准偏移值精度最高。而将克拉美罗下界推导出的方位向偏移值方差$\sigma _{{\rm{\Delta }}x}^2$作为路径长度进行叠加,获得的最短路径距离总和却能够代表方位向偏移值的最高精度。由于各辅影像的头文件和轨道姿态的差异及影像间噪声特征的各异性,计算出的方位向偏移值是相互独立的。根据误差传播理论能够计算出主辅影像方位向偏移值的最小方差,即:
$$\sigma _{{\rm{\Delta }}x}^2 = \mathop \sum \limits_{i = 1}^n \sigma _{{\rm{\Delta }}{x_i}}^2$$ (14) 式中,n为路径中所有节点的数量; $\sigma _{{\rm{\Delta }}{x_i}}^2$为节点i的配准像元方差。
因此,本文提出利用克拉美罗下界推导出的方位向配准偏移方差作为像对间的“长度”,以规划辅影像与主影像间的最小配准方差路径,同时达到缩小配准方差与提高配准效率的目的。针对最短路径的搜索,提出利用迪杰斯特拉算法进行配准路径的规划。
迪杰斯特拉算法采用贪心策略,利用连接点之间的距离松弛整体路径长度,获得主辅影像间的最小配准方差。本文对18幅影像节点进行了最短路径搜索,生成了如图 4(c)所示的迪杰斯特拉时空基线图,各辅影像节点与主影像节点10之间的连线代表配准偏移传递的最短路径。
搜索辅影像节点与主影像节点之间最短路径的核心思想如下:网络路径距离为像对间理论相干性计算出的ESD偏移量方差,影像起点和终点间的最短路径${D_{i, j}}$定义为:
$${D_{i, j}} = \left\{ {\begin{array}{*{20}{l}} {0, i = j}\\ {{\rm{arg \;min}}\left( {{D_{i, k}} + {D_{k, j}}} \right), i \ne j} \end{array}} \right.$$ (15) 式中,k节点为与i节点相连的节点,根据贪心策略,当前最短路径的搜索过程不需要考虑整体的距离总和,只需在当前节点搜索起始点的所有连接点。获得${D_{i.j}}$的同时将i节点置入标记集合L,继续从k节点开始搜索所有的连接节点,如节点出现在集合L中,则跳过搜索,此时搜索原则也应满足式(15),搜索完毕后,同样将k节点置入集合L中。在加入贪心策略获得节点间距离后,总是保持${D_{i.j}}$最短的过程实际为利用节点间距离对整体距离进行松弛的过程。重复上述步骤,直至集合L包含所有的节点,则最短路径搜索过程完毕。
3 实验与分析
3.1 目视检译
为了对比时序像对配准过程中高相干点ESD与传统ESD方法的精度,选择四川盆地作为实验区域。该地区地势起伏较大,且因气候原因,植被覆盖率较高,可以作为比较有代表性的配准实验区。影像拍摄时间为2016-01-04与2016-01-16。首先基于精轨数据与数字高程模型对影像进行几何配准,再利用传统ESD与高相干点方法进行精配准。如果直接对粗配准后的影像进行拼接操作,由于几何配准的精度限制,出现了相位跳变的情况,如图 5(a)所示,其中A-A'表示干涉图相位截面。图 5(b)为高相干点ESD校正后的干涉图,图 5(c)为A-A'的相位横截面(虚线方框代表干涉图相位截面中的相位跳变所在区域)。
可以看到,相比传统方法,高相干点方法弥补了相位跳变,说明利用精确估计的相干性提取的高相干点进行ESD时,其效果比未经相位补偿的相干性效果要好。传统ESD中,受水体影响,大面积包含水体的子块计算出的方位向偏移绝对值通常都会大于0.01像元,结果表现为粗差,从而对整体偏移的计算产生影响;而高相干点通常不受水体影响,且时间相干性较高,比较稳定。
3.2 时序精配准实验对比
结合效果更好的高相干点方法进行时序配准实验,分别利用苏州地区18幅升轨影像和太行山脉东部地区(简称太行地区)18幅升轨影像进行精配准。选择这些研究区域有以下原因:
1)苏州地区作为平原研究区域,富含水体,对传统ESD单像对配准造成的影响很难用传统粗差探测方式解决,可以作为传统ESD配准方法无法达到理论精度的典型案例。
2)虽然影像获取时间包含树木繁盛的季节,但城市区域包含较多的高相干点,植被对ESD的影响较小,苏州地区作为高相干点丰富的理想实验区,用以验证迪杰斯特拉算法在时序影像ESD精配准(本文简称为“迪杰斯特拉ESD”)中的可行性。
3)太行地区平均海拔高于1 000 m,地势起伏较大,受亚热带季风气候影响,山脉植被多为乔木林和针叶林,植被覆盖率较大,对相干性的影响比较严重。此地区可以作为较恶劣条件下的配准实验区,用以验证迪杰斯特拉ESD在山区的可行性。
利用SM-ESD、NESD、迪杰斯特拉ESD 3种方法分别对苏州地区和太行地区的时序影像进行精配准,各辅影像传递后的配准偏移量对比如图 6所示。图 6(a)中,迪杰斯特拉ESD生成的配准像对传递后的配准偏移量与NESD及SM-ESD方法对比,没有超过0.001像元的偏移误差,这说明迪杰斯特拉ESD在平原地区是有效的。针对苏州地区的辅影像16,迪杰斯特拉ESD方法估计出的偏移量与NESD的结果偏差较大,利用迪杰斯特拉ESD传递的配准偏移量纠正前后的干涉图如图 7所示。图 7(a)为苏州地区辅影像16与主影像生成的干涉图,图 7(b)为应用迪杰斯特拉ESD纠正后的干涉图,图 7(c)为相位横截面(虚线方框代表干涉图相位截面中的相位跳变所在区域)。可见相位跳变已消除,证明了迪杰斯特拉ESD在此区域的有效性。图 6(b)中,太行地区辅影像10和辅影像13处迪杰斯特拉ESD估计出的影像配准偏移量与NESD方法估计出的偏移量存在较大偏差,对比应用两种方法纠正后重叠区域的残余干涉相位差均值,结果如图 8所示。
图 8展示了重叠区域偏移纠正后各子块残余的干涉相位差均值,可以看出,应用NESD偏移纠正后的相位残差均值约为迪杰斯特拉ESD的5倍,证明迪杰斯特拉ESD在太行地区辅影像10、13处的偏移估计精度高于NESD。实际上,相干性较高的主辅影像经单次计算即可获得较为精确的偏移量,而NESD对小基线集中所有偏移量进行加权最小二乘,实际上降低了较高相干性像对间的偏移计算精度。迪杰斯特拉ESD根据配准偏移方差进行配准像对的规划,在主辅影像相干性较好的情况下,直接根据较可靠的单次偏移估计值进行精配准,弥补了NESD的这个缺点。
应用图 6中迪杰斯特拉ESD的方位向配准偏移估计值纠正后,实验区域的去平(de-flatten,去除平地相位)干涉图雷达波束影像拼接结果如图 9所示。
从图 9可以看出,所有去平干涉图均未出现相位跳变现象,为进一步证明迪杰斯特拉ESD的有效性,定量说明迪杰斯特拉ESD的配准精度,本文设计了与NESD方法的对比实验。针对低相干性的太行地区,在前一次计算出的配准偏移量$\Delta x$基础上,利用$\Delta x$进行校正,校正后序贯估计出配准偏移量$\Delta x$,重复进行此迭代实验100次。为方便与NESD的对比,将偏移量传递至时间节点的第1幅影像,估计出的误配准偏移量箱式图如图 10所示。
箱式图能够以统计特征的形式表示两种不同配准方法估计误配准像元的迭代收敛能力。从图 10可以看出,NESD估计配准像元偏移的能力随着相干性的降低超出至0.001像元,而迪杰斯特拉ESD则始终能够保持在0.001像元以内的精度,此实验证明了迪杰斯特拉ESD的精度高于NESD。
除上述精度优势外,计算效率也是Sentinel-1A TOPS模式影像配准过程中需要关注的热点问题。本文实验配置为Windows 10、i7 3.2 GHz处理器,软件环境为MATLAB 2016a。表 1给出了苏州地区和太行地区的影像配准计算效率, 其中,苏州地区和太行地区的影像大小(方位向×距离向)分别为4353×23362像素、4596×22130像素。
表 1 影像精配准计算效率Table 1. Computational Efficiency of Image Co⁃registration实验区域 NESD 迪杰斯特拉ESD 偏移计算次数 耗时/s 偏移计算次数 耗时/s 苏州地区 153 692.4 34 144.2 太行地区 153 715.3 22 102.5 由表 1可知,17幅辅影像精配准至主影像的过程中,相比NESD,迪杰斯特拉ESD在苏州地区(平原地区)的偏移计算次数和耗时减少了约5倍,在太行地区(山地地区)减少了约7倍,能够满足快速监测背景下的配准效率要求。
4 结 语
ESD是目前哨兵卫星影像使用最广泛的精配准方式,相干性作为ESD的重要参数,其估计的优劣程度直接决定了哨兵影像的配准精度。除干涉基线外,地形、形变等相位贡献导致的相位梯度也是ESD精度的重要影响因素。本文提出在进行相位补偿之后,利用精确估计的相干性提取高相干点,提高了ESD精配准过程中的选点精度,这有益于高相干区域的影像配准。时序影像精配准中采用高相干点ESD,相比传统ESD方法,精度有所提升。针对时序影像精配准过程,提出应用迪杰斯特拉最短路径算法,将克拉美罗下界推导出的方位向配准偏移方差作为路径长度,以此规划配准像对。高相干点ESD结合迪杰斯特拉算法的优势在于不仅避免了繁杂的配准数,同时也提高了配准精度与配准效率。就迪杰斯特拉ESD配准精度而言,对比小基线集NESD配准方式,高相干性像对间只需采用单次配准偏移计算就能获得精确的结果,避免了加权最小二乘对精确结果的扰动。对比应用两种方法后雷达波束影像重叠区域的残余干涉相位差,迪杰斯特拉ESD方法的相位残差均值比NESD方法减少了约400%。就配准效率而言,迪杰斯特拉ESD的配准次数约为NESD的1/5(仅在平原地区),计算时间能够降低400%以上,为应急背景下的快速监测提供了良好的基础条件。
致谢: 感谢欧洲空间局通过哥白尼计划提供的Sentinel⁃1系列卫星数据。 -
表 1 影像精配准计算效率
Table 1 Computational Efficiency of Image Co⁃registration
实验区域 NESD 迪杰斯特拉ESD 偏移计算次数 耗时/s 偏移计算次数 耗时/s 苏州地区 153 692.4 34 144.2 太行地区 153 715.3 22 102.5 -
[1] Torres R, Snoeij P, Geudtner D, et al. GMES Sentinel-1 Mission[J]. Remote Sensing of Environment, 2012, 120:9-24 doi: 10.1016/j.rse.2011.05.028
[2] Geudtner D, Torres R. Sentinel-1 System Overview and Performance[C]. 2012 IEEE International Geoscience and Remote Sensing Symposium, Munich, Germany, 2012
[3] De Zan F, Guarnieri A M. TOPSAR:Terrain Observation by Progressive Scans[J]. IEEE Transactions on Geoscience and Remote Sensing, 2006, 44(9):2352-2360 doi: 10.1109/TGRS.2006.873853
[4] Sansosti E, Berardino P, Manunta M, et al. Geometrical SAR Image Registration[J]. IEEE Transactions on Geoscience and Remote Sensing, 2006, 44(10):2861-2870 doi: 10.1109/TGRS.2006.875787
[5] De Zan F. Accuracy of Incoherent Speckle Tracking for Circular Gaussian Signals[J]. IEEE Geoscience and Remote Sensing Letters, 2013, 11(1):264-267 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=820f951d82b3f68db785abb6bc0be039
[6] 吴文豪, 周志伟, 李陶, 等. 精密轨道支持下的哨兵卫星TOPS模式干涉处理[J]. 测绘学报, 2017, 46(9):1156-1164 http://d.old.wanfangdata.com.cn/Periodical/chxb201709011 Wu Wenhao, Zhou Zhiwei, Li Tao, et al. A Study of Sentinel-1 TOPS Mode Co-registration[J]. Acta Geodaetica et Cartographica Sinica, 2017, 46(9):1156-1164 http://d.old.wanfangdata.com.cn/Periodical/chxb201709011
[7] 李振洪, 宋闯, 余琛, 等. 卫星雷达遥感在滑坡灾害探测和监测中的应用:挑战与对策[J]. 武汉大学学报·信息科学版, 2019, 44(7):967-979 doi: 10.13203/j.whugis20190098 Li Zhenhong, Song Chuang, Yu Chen, et al. Application of Satellite Radar Remote Sensing to Landslide Detection and Monitoring:Challenges and Solutions[J]. Geomatics and Information Science of Wuhan University, 2019, 44(7):967-979 doi: 10.13203/j.whugis20190098
[8] Tian X, Malhotra R, Xu B, et al. Modeling Orbital Error in InSAR Interferogram Using Frequency and Spatial Domain Based Methods[J]. Remote Sensing, 2018, 10(4):508 doi: 10.3390/rs10040508
[9] Scheiber R, Moreira A. Coregistration of Interferometric SAR Images Using Spectral Diversity[J]. IEEE Transactions on Geoscience and Remote Sensing, 2000, 38(5):2179-2191 doi: 10.1109/36.868876
[10] Prats-Iraola P, Scheiber R, Marotti L, et al. TOPS Interferometry with TerraSAR-X[J]. IEEE Transactions on Geoscience and Remote Sensing, 2012, 50(8):3179-3188 doi: 10.1109/TGRS.2011.2178247
[11] Bamler R, Eineder M. Accuracy of Differential Shift Estimation by Correlation and Split-Bandwidth Interferometry for Wideband and Delta-k SAR Systems[J]. IEEE Geoscience and Remote Sensing Letters, 2005, 2(2):151-155 doi: 10.1109/LGRS.2004.843203
[12] Yague-Martinez N, Prats-Iraola P, Gonzalez F R, et al. Interferometric Processing of Sentinel-1 TOPS Data[J]. IEEE Transactions on Geoscience and Remote Sensing, 2016, 54(4):2220-2234 doi: 10.1109/TGRS.2015.2497902
[13] Veci L, Lu J, Foumelis M, et al. ESA's Multi-mission Sentinel-1 Toolbox[C]. The EGU General Assembly Conference, Vienna, Austria, 2017
[14] Kampes B, Usai S. Doris:The Delft Object-Oriented Radar Interferometric Software[C]. The 2nd International Symposium on Operationalization of Remote Sensing, Enschede, Netherlands, 1999
[15] Rosen P A, Gurrola E, Sacco G F, et al. The InSAR Scientific Computing Environment[C]. The 9th European Conference on Synthetic Aperture Radar, Nuremberg, Germany, 2012
[16] Fattahi H, Agram P, Simons M. A Network-Based Enhanced Spectral Diversity Approach for TOPS Time-Series Analysis[J]. IEEE Transactions on Geoscience and Remote Sensing, 2016, 55(2):777-786 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=f2005014546cd683023c17e807ec04c0
[17] Yague-Martinez N, de Zan F, Prats-Iraola P. Coregistration of Interferometric Stacks of Sentinel-1 Tops Data[J]. IEEE Geoscience and Remote Sensing Letters, 2017, 14(7):1002-1006 doi: 10.1109/LGRS.2017.2691398
[18] Xu B, Li Z, Feng G, et al. Continent-Wide 2-d Co-seismic Deformation of the 2015 Mw 8.3 Illapel, Chile Earthquake Derived from Sentinel-1A Data:Correction of Azimuth Co-registration Error[J]. Remote Sensing, 2016, 8(5):376 doi: 10.3390/rs8050376
[19] Berardino P, Fornaro G, Lanari R, et al. A New Algorithm for Surface Deformation Monitoring Based on Small Baseline Differential SAR Interferograms[J]. IEEE Transactions on Geoscience and Remote Sensing, 2002, 40(11):2375-2383 doi: 10.1109/TGRS.2002.803792
[20] Dijkstra E W. A Note on Two Problems in Connexion with Graphs[J]. Numerische Mathematik, 1959, 1(1):269-271 doi: 10.1007/BF01386390
[21] Rocca F. Modeling Interferogram Stacks[J]. IEEE Transactions on Geoscience and Remote Sensing, 2007, 45(10):3289-3299 doi: 10.1109/TGRS.2007.902286
[22] Guarnieri A M, Tebaldini S. Hybrid Cramér-Rao Bounds for Crustal Displacement Field Estimators in SAR Interferometry[J]. IEEE Signal Processing Letters, 2007, 14(12):1012-1015 doi: 10.1109/LSP.2007.904705
[23] Jiang M, Ding X, Li Z, et al. InSAR Coherence Estimation for Small Data Sets and Its Impact on Temporal Decorrelation Extraction[J]. IEEE Transactions on Geoscience and Remote Sensing, 2014, 52(10):6584-6596 doi: 10.1109/TGRS.2014.2298408
[24] Jiang M, Ding X, Li Z. Hybrid Approach for Unbiased Coherence Estimation for Multitemporal InSAR[J]. IEEE Transactions on Geoscience and Remote Sensing, 2013, 52(5):2459-2473 https://www.researchgate.net/publication/261296082_Hybrid_Approach_for_Unbiased_Coherence_Estimation_for_Multitemporal_InSAR
[25] 蒋弥, 丁晓利, 李志伟, 等. 基于时间序列的InSAR相干性量级估计[J]. 地球物理学报, 2013, 56(3):799-811 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=dqwlxb201303009 Jiang Mi, Ding Xiaoli, Li Zhiwei, et al. InSAR Coherence Magnitude Estimation Based on Data Stack[J]. Chinese Journal of Geophysics, 2013, 56(3):799-811 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=dqwlxb201303009
[26] Wang T, Liao M, Perissin D. InSAR Coherence-Decomposition Analysis[J]. IEEE Geoscience and Remote Sensing Letters, 2009, 7(1):156-160 http://d.old.wanfangdata.com.cn/Periodical/hwyhmb201604013
[27] Refice A, Bovenga F, Nutricato R. MST-Based Stepwise Connection Strategies for Multipass Radar Data, with Application to Coregistration and Equalization[J]. IEEE Transactions on Geoscience and Remote Sensing, 2006, 44(8):2029-2040 doi: 10.1109/TGRS.2006.872907
-
期刊类型引用(3)
1. 冯莉,吕修凯,崔生乐,杨春梅,徐晓燕. 基于双向迪杰斯特拉算法移栽机补栽路径规划及仿真. 中国农机化学报. 2023(03): 177-182 . 百度学术
2. 艾强,曾春琴. 时序InSAR干涉组合网络优化方法及适应性分析. 工程勘察. 2023(08): 67-71+78 . 百度学术
3. 马张烽,蒋弥,丁琪瑄. 时序Sentinel-1TOPS模式SAR数据精配准. 测绘学报. 2021(05): 634-640 . 百度学术
其他类型引用(6)