作者中心
专家审稿
责编办公
编辑办公
-
摘要: 小型化、轻量化的无人机(UAV)为合成孔径雷达(SAR)提供了更加灵活、机动的观测平台，无人机载干涉SAR (InSAR)逐步应用于干涉测量领域。无人机小而轻，易受气流扰动的影响，采用多航过模式进行干涉测量时，飞行航迹非线性且不平行。非线性、不平行的飞行轨迹导致两幅图像之间存在几何畸变。在复杂地形条件下，无人机载InSAR的干涉图像对之间的偏移量大且具有明显的空变特性，给图像配准带来了很大的技术挑战，常规的基于多项式拟合的配准方法不再适用。该文提出了一种利用地形辅助分区的图像配准方法。首先基于航迹信息生成高程门限，利用外部地形对测量区域进行图像分区处理，然后对区域内的偏移量构建多项式变换模型，对各区域边界处的偏移量施加约束，并进行联合求解，最后获得连续的全局偏移量拟合面，通过对辅图像进行重采样实现精配准。基于P波段无人机载InSAR获取的实测数据，初步验证了该方法的有效性。
-
关键词:
- 无人机载干涉合成孔径雷达 /
- 图像配准 /
- 复杂地形 /
- 多项式模型 /
- 图像分区
Abstract: Miniaturized and lightweight Unmanned Aerial Vehicles (UAV) provide a flexible platform for Synthetic Aperture Radar (SAR). The application of UAV Interferometric SAR (InSAR) is gradually increasing in interferometric measurement fields. UAVs are small and light, which are easily affected by airflow disturbances. Their trajectories are nonlinear and unparallel when adopting the multipass mode for interferometry. The nonlinear and unparallel trajectories result in geometric distortion between the interferometric image pairs. Under complex topography conditions, the interferometric image pairs of UAVs have large offsets that are obviously space-dependent, thereby resulting in substantial technical challenges during image registration. Conventional image registration methods based on polynomial fitting are no longer applicable. In this study, we proposed an image registration method based on image partition with topography assistance. First, an elevation threshold is generated based on the UAV trajectories, and the measurement area is partitioned using the assisted topography. Then, a polynomial fitting model is constructed for offsets within each partition with constraints applied at the partition boundaries for joint optimization. Finally, continuous global offset fitting surfaces are obtained, and precise image registration is achieved by resampling the slave image. The effectiveness of the proposed method is preliminarily validated using real measurement data obtained from UAV InSAR in the P-band. -
1. 引言
干涉合成孔径雷达(Interferometric Synthetic Aperture Radar, InSAR)广泛应用于滑坡监测和地形测绘等领域[1]。InSAR图像配准是指根据两幅图像同名点之间的坐标映射关系，将辅图像重采样到与主图像相同的像素网格，使配准后两幅图像的同名点对位于同一分辨单元内[2]。InSAR产生的复图像必须经过精确配准才能保证图像对之间具有良好的相干性。因此，配准作为干涉测量中重要的一环，直接影响到干涉图的质量和地形测量的精度[3,4]。
针对InSAR图像的配准，其核心是偏移量的提取与估计。偏移量的提取与估计方法一般分为两类，一类是利用较为精准的外部轨道等信息，结合干涉测量几何直接计算同名点对之间的偏移量[5]。这类方法通常以成像区域的平均高程作为整个区域的统一高程，即需要满足平地假设[6]。相对于飞行高度，当成像场景相对平坦，没有大高程起伏时，配准精度能达到亚像素级。因此，这类方法通常应用于星载InSAR测量[7]。另一类方法则是基于InSAR图像中的信息来估计偏移量。这类方法主要采用最大互相关系数(Maximum Cross Correlation, MCC)[8]、最大谱函数(Maximum Spectrum Function, MSF)[9]或最小平均波动函数(Average Fluctuation Function, AFF)[10]等作为测度函数，基于滑动窗口搜索，实现同名点对之间的偏移量估计[11]。根据提取到的偏移量，基于所选像素点的坐标位置构建多项式拟合模型，通过最小二乘法估计模型参数，建立起从主图像到辅图像之间的映射关系[12]。在实际处理中，基于多项式拟合的方法需根据干涉图像对之间偏移量具体情况选择合适的多项式模型阶数，但对于空变的偏移量易失效。
与传统的InSAR配准算法不同，其他利用SAR图像信息的配准方法不再使用测度函数来匹配同名点，而是基于SAR图像中的点[13]、线[14]和边缘结构[15]等显著几何特征来估计配准偏移量。然而，传统的特征提取方法十分复杂，且难以获得更全面和更有表现力的特征[16]。随着深度学习的发展，卷积神经网络(Convolution Neural Network, CNN)逐步应用于图像配准中的特征提取和匹配[17]。这些方法更适合于具有明显变化的图像处理场景，例如多角度图像配准[18]和异源图像配准[19]等，需要相对固定的成像场景[20]和大量的图像数据集进行训练，不适合InSAR图像配准。
与机载InSAR相比，无人机体积小，飞行易受气流扰动的影响，航迹非线性[21]；载重有限，通常仅负载单天线，需要多次飞行才能实现干涉测量，航迹不平行[22]。非线性、不平行的飞行轨迹导致两幅SAR图像存在几何畸变，偏移量变化复杂。另外，无人机飞行高度有限，测量区域的高程与飞行高度相差不大，星载InSAR或常规机载InSAR测量中的平地近似假设不再成立。在复杂地形条件下，干涉图像对的不同同名像素之间的偏移量具有明显的空变特性，特别是在长基线模式下，平地的区域偏移量较小，而在有高程区域，偏移量可以达到数个像素[23]。当成像区域内高程起伏较大时，偏移量空变明显，传统的全局多项式拟合配准方法不再适用[24,25]。
本文围绕无人机载InSAR图像配准问题开展了研究，分析了无人机载InSAR成像投影几何模型，建立了适用于无人机载InSAR的配准偏移量模型，论证了地形对无人机载InSAR图像配准的影响。为了解决复杂地形条件下的配准难题，本文提出了一种基于图像分区的精配准算法。首先，基于无人机载InSAR偏移量模型，结合航迹信息，通过设定偏移量门限得到高程门限；其次，利用外部地形数据进行图像分区，并进行连通域分析，剔除面积较小的区域；最后，构建带边界约束的子块多项式全局拟合模型，引入Lagrange算子，利用最小二乘法联合求解模型参数，并进行图像重采样获取精配准后的图像。利用P波段无人机载InSAR实测数据验证了该方法的有效性，实验数据结果表明配准后的相干性大大提高，残差点数量减少。在高程较大的区域，干涉相位的质量显著提高。
2. 无人机载SAR成像几何模型
常规线性孔径SAR工作在正侧视条带模式下，天线波束指向始终与平台运动方向垂直[26]。图1所示为常规线性孔径SAR成像几何示意图。常规线性孔径SAR成像时，是将观测区域的三维地形投影在二维成像平面上。以平台运动方向作为x方向，竖直向上方向作为z方向，根据右手准则，构建三维成像直角坐标系。其中，x-y平面为成像平面，雷达运动轨迹的中心位置S为坐标原点，轨迹两端分别为S1和S2[27]。
假设观测场景中一地形点
P(xp,yp,zp) ，其在成像平面上的投影点为P′(x′p,y′p,0) 。雷达沿着一平行于x轴的直线运动，根据后向投影(Back Projection, BP)成像算法，地形点P(xp,yp,zp) 与其对应的投影点P′(x′p,y′p,0) 的距离徙动曲线始终保持完全一致，可以得到投影点P′(x′p,y′p,0) 的坐标为{x′p=xpy′p=√y2p+(H−zp)2−H2=√y2p+z2p−2Hzp (1) 其中，H代表无人机平台相对于成像平面的高度。
在实际的无人机载SAR系统中，平台易受到气流扰动等外界因素的影响，难以实现匀速直线运动，呈现出非线性的航迹特性。图2所示为实际无人机载SAR成像几何。此时，天线相位中心(Antenna Phase Center, APC)与地形点
P(xp,yp,zp) 之间的距离徙动曲线可表示为R(t)=√(vt+Δx(t)−xp)2+(Δy(t)−yp)2+(Δz(t)+H−zp)2 (2) 其中，
(Δx(t),Δy(t),Δz(t)) 代表偏离匀速直线运动轨迹的三维运动误差。雷达APC与投影点
P′(x′p,y′p,0) 之间的距离徙动曲线可表示为R′(t)=√(vt+Δx(t)−x′p)2+(Δy(t)−y′p)2+(Δz(t)+H)2 (3) 在非线性航迹条件下，当地形点
P(xp,yp,zp) 在成像平面上时，即zp=0 时，地形点P(xp,yp,zp) 与投影点P′(x′p,y′p,0) 重合，此时地形点P(xp,yp,zp) 和投影点P′(x′p,y′p,0) 的距离徙动曲线始终保证完全一致，可以精确聚焦成像。当地形点P(xp,yp,zp) 不在成像平面上时，即zp≠0 时，地形点P(xp,yp,zp) 和投影点P′(x′p,y′p,0) 的距离徙动曲线无法保证完全一致，此时，投影点P′(x′p,y′p,0) 的坐标无法利用成像几何，直接通过式(1)进行求解。为确定三维空间中一个地形点
P(xp,yp,zp) 在成像平面上的投影点P′(x′p,y′p,0) ，本文以地形点P(xp,yp,zp) 与投影点P′(x′p,y′p,0) 的距离徙动曲线之差的平方和作为目标函数，(xn,yn,zn) 代表非线性航迹上任意点的坐标，采用迭代数值求解的方式来估计该目标函数的最小值，从而得到该投影点的坐标。构建目标函数如式(4)所示：argminx′p,y′p(∑n((xn−x′p)2+(yn−y′p)2+z2n−(xn−xp)2−(yn−yp)2−(zn−zp)2)2) (4) 实际数值求解时，通过提取无人机非线性航迹的线性分量，基于无人机载SAR成像几何，根据式(1)计算出该地形点
P(xp,yp,zp) 在线性分量下的投影点，作为迭代的初始解。利用仿真地形和实测航迹对投影点位置获取方法进行验证。实测航迹如图3所示，在线性分量附近无规律波动。仿真地形如图4所示，仿真地形中存在地形起伏。
图5所示为迭代解与真实解之差。真实解是指利用BP成像算法，对仿真地形进行成像处理，将三维地形点投影至二维平面上，提取的各地形点的成像幅度峰值点的坐标位置。迭代解是指结合航迹信息，利用式(4)对仿真地形上的各地形点进行迭代处理，得到的投影点位置。可以发现，迭代出的近似解与真实投影点位置的差异在10–13 m量级，完全可以忽略。因此，迭代得到的近似解可以作为地形点在二维平面上的投影点位置。
图 5 迭代解与真实投影点位置误差Figure 5. Error between the iteration solution and the real projection point与此同时，在该条航迹下，对比非线性航迹中的线性分量得到的初始解和迭代后得到的近似解，其二者差异在10–5 m量级，如图6所示。因此，在一定的约束条件下，可以用非线性航迹中线性分量下的投影点位置作为非线性航迹的真实投影点位置。
下面分析以线性分量下投影点位置代替非线性航迹投影点位置的边界约束条件。实际上，在应用式(4)进行投影点位置获取时，假设在成像过程中，地形点
P(xp,yp,zp) 在二维平面上的成像结果不散焦。此时，可利用非线性航迹中线性分量下的投影点位置作为非线性航迹的投影点位置。因此，可将非线性航迹下地形点P(xp,yp,zp) 在二维平面上的成像结果不散焦的条件作为快速获取投影点位置的边界条件。在进行BP成像时，相干叠加时多普勒相位不能完全补偿，存在相位误差，严重时导致成像散焦。斜距误差历程可以表示为
ΔR(t)=R(t)−R′(t) (5) 忽略沿航迹方向的运动误差
Δx(t) ，将式(1)代入式(5)进行近似化简可得ΔR(t)=−Δy(t)Rp(yp−y′p)−Δz(t)Rpzp (6) 其中，
Rp=√y2p+(H−zp)2 ，代表的是天线相位中心到地形点P(xp,yp,zp) 的最短斜距。假设目标处在波束中心位置时，APC与投影点
P′(x′p,y′p,0) 和地形点P(xp,yp,zp) 之间的斜距相等，则有√(Δy(t)−y′p)2+(Δz(t)+H)2=√(Δy(t)−y′p−δ)2+(Δz(t)+H−zp)2 (7) 其中，
δ=yp−y′p ，可得δ=H+Δz(t)Δy(t)−y′pzp (8) 将式(8)代入式(6)，进一步整合可得
ΔR(t)=−zpy′p(Δy(t)cosβ+Δz(t)sinβ)=−zpy′pΔr(t)cos(β−α) (9) 其中，
Δr=√(Δy)2+(Δz)2 ，代表实际雷达APC与参考位置之间的运动误差，α 代表实际雷达APC与参考位置连线与水平轴的夹角，β 代表目标下视角。具体几何关系如图7所示。假设最大斜距误差为
ΔRm ，最大运动误差为Δrm 。通常认为对于点目标成像时，相位误差大于π/4时会影响成像质量。因此，可以假设满足成像质量的条件为zpy′pΔrm≤λ16 (10) 其中，
λ 为发射信号波长。当航迹确定时，对于三维空间中任意目标点，其投影点位置也确定，成像处理点目标不散焦时，能量聚集于投影点位置。因此，当式(10)成立时，目标真实投影点位置等效于航迹线性分量对应的投影点位置。
实际成像处理时，即使不考虑雷达位置的测量误差，对于起伏地形成像时，波长较小的高频段 SAR很难满足式(10)所示的边界条件，而工作在较低频段的SAR则较容易满足。图8所示为地形与最大斜距误差的关系。当最大运动误差和目标地距一定时，最大斜距误差随目标高度的增大而增大；当目标高度一定时，最大斜距误差随着地距的变大而减小，随着最大运动误差的增大而增大。当运动误差越小，目标地距越大时，成像质量受目标高度的影响越小。图8中
λ 为1 m。
图 8 不同最大运动误差和目标地距下最大斜距误差与目标高程之间的关系Figure 8. Relationship between the maximum slant range error and the target height with different the maximum motion error and the target ground distance3. 外部地形辅助分区的无人机载InSAR图像配准方法
3.1 无人机载InSAR偏移量建模
基于前文推导的无人机载SAR成像投影几何模型与边界条件可知，在较低频段下，利用导航系统提供的雷达位置信息，成像投影点位置可以等效于非线性航迹中线性分量对应的投影点位置。
图9所示为无人机载InSAR干涉测量投影几何。对于地形点
P(xp,yp,zp) ，分别用两次航过对其进行成像处理。两次成像的合成孔径中心分别为C1(xc1,yc1,zc1) 和C2(xc2,yc2,zc2) ，两次成像的投影点位置分别为P1(x′p1,y′p1,0) 和P2(x′p2,y′p2,0) 。以其中一次航过为例，提取非线性航迹中的线性分量，假设其单位方向向量为
(m0,n0,p0) ，直线上一点为(x1,y1,z1) ，直线方程可表示为{x=x1+m0ty=y1+n0tz=z1+p0t (11) 其中，t为直线方程的参数。
为了获取
P1(x′p1,y′p1,0) 的位置坐标，将线性分量上任意两点和地形点P(xp,yp,zp) 构成的三角形，绕着该线性分量旋转，与成像平面的交点位置即为投影点的位置。所形成的弧⌢PP1 所在平面的平面方程可写为m0(x−xp)+n0(y−yp)+p0(z−zp)=0 (12) 此时与地形点
P(xp,yp,zp) 对应的合成孔径中心C1(xc1,yc1,zc1) 即为弧⌢PP1 所在平面与直线方程(11)的交点，即C1(xc1,yc1,zc1) 满足m0(xp−x1)+n0(yp−y1)+p0(zp−z1)=t0 (13) 其中，
t0 为C1 在直线方程中对应的参数。因此，
C1(xc1,yc1,zc1) 可表示为{xc1=x1+m0t0yc1=y1+n0t0zc1=z1+p0t0 (14) 弧的半径r可表示为
r=√(xc1−xp)2+(yc1−yp)2+(zc1−zp)2 (15) P1(x′p1,y′p1,0) 位于直线C1P 在x-y平面的投影上，即位于(xp,yp,0) 与(xc1,yc1,0) 的连线上，且与(xc1,yc1,0) 的距离为√r2−z2c1 。因此，P1(x′p1,y′p1,0) 可表示为{x′p1=√r2−z2c1cosθ+xc1y′p1=√r2−z2c1sinθ+yc1 (16) 其中，
cosθ=xp−xc1√(xp−xc1)2+(yp−yc1)2 (17) 将式(15)和式(17)代入式(16)化简可得
{x′p1=(xp−xc1)⋅√1+(zc1−zp)2−z2c1(xp−xc1)2+(yp−yc1)2+xc1y′p1=(yp−yc1)⋅√1+(zc1−zp)2−z2c1(xp−xc1)2+(yp−yc1)2+yc1 (18) 当孔径为线性且平行于飞行方向时，即
xp=xc1 ,yc1=0 ,zc1=H 时，式(18)与式(1)一致。进一步地，对于地形点
P(xp,yp,zp) 两次航过的偏移量模型可以表示为{Δx=xc1−xc2+(xp1−xc1)⋅(1−√1+2zp(zc1−zc2)(xp1−xc1)2+(yp1−yc1)2)Δy=yc1−yc2+(yp1−yc1)⋅(1−√1+2zp(zc1−zc2)(xp1−xc1)2+(yp1−yc1)2) (19) 由式(19)可以看出，无人机飞行高度受限，与实际地形相差不大，星载InSAR或常规机载InSAR测量中的平地近似假设不再成立。干涉图像对的偏移量与地形相关，具有很强的空变性。特别是在长基线模式下，平地的区域偏移量较小，而在有高程区域，偏移量可以达到数个像素。当成像区域内高程起伏较大时，偏移量出现明显的空变现象，传统的配准方法失效。因此，基于多项式的全局偏移量拟合模型无法准确描述真实的偏移量情况，需要对不同地形区域分别处理，才能获取较好的配准效果。
3.2 基于多级高程门限的图像分区
基于BP成像算法，对两次SAR回波数据在同一成像坐标系下进行成像处理[28]。由式(19)可知，当航迹固定时，任意三维空间内的目标，其二维配准偏移量也随之确定。那么已知目标的三维空间位置即可计算出该点的二维配准偏移量。
首先，提取非线性航迹中的线性航迹分量，利用式(14)获取各像素点合成孔径中心位置
(xc,yc,zc) 。然后，根据式(19)，对于任一像素点，其二维配准偏移量与目标高程一一对应，设置偏移量门限，代入成像平面中各像素点的位置，从而获取各像素点的高程门限。通常认为，当偏移量小于1个像素时，使用多项式拟合的方法可以较好地实现全局偏移量的拟合。因此，需设置多级高程门限进行图像分区，而后对子区域内进行多项式拟合，才能更准确地表示区域内的偏移量情况。最后，利用地形数据与获取的多级高程门限相比，将成像平面进行区域划分，划分为低于门限的接近成像平面的高程区域和高于门限的其他高程区域。对分块后的结果进行连通域分析，剔除面积较小的区域，得到最终的图像分区结果。
实际处理过程中，在外部地形的精度不足时，不能保证分区的完全准确。本文对图像分区时采用亚像素级的偏移量作为分区门限，实际各区不可避免地会存在部分像素点的偏移量超过门限，但是本文方法核心思想是将全图的偏移量拟合问题，转化为若干个子区域的偏移量在约束条件下的联合拟合求解问题，分区误差对局部偏移量的多项式拟合影响有限。
3.3 带约束的子块多项式拟合
无人机载InSAR的二维配准偏移量与地形有关。对于整个场景区域，其二维配准偏移量应该是连续的。进行图像分区后，如果直接采用常规的多项式拟合方法来获取各个子区域内的偏移量，由于边界处的拟合偏移量往往过大或过小，易导致各子区域边界连接处偏移量出现较大“跳变”，因此本文提出了带约束的子块多项式拟合方案。
以一阶多项式拟合为例，假设图像分为N个子区域，并且子区域内M个参考点的二维坐标和二维配准偏移量已知，可以得到
{an1xm,n+bn1ym,n+cn1=Δxm,nan2xm,n+bn2ym,n+cn2=Δym,n (20) 其中，
(xm,n,ym,n) 代表第n个子区域第m个参考点在雷达图像中的二维坐标，(Δxm,n,Δym,n) 代表第n个子区域第m个参考点的二维配准偏移量。an1,bn1,cn1,an2,bn2,cn2 分别为待估计参数。对于区域边缘的像素，需要附加约束条件，使二维配准偏移量连续，即需满足：
{an1xm,nq+bn1ym,nq+cn1=aq1xm,nq+bq1ym,nq+cq1an2xm,nq+bn2ym,nq+cn2=aq2xm,nq+bq2ym,nq+cq2 (21) 其中，
(xm,nq,ym,nq) 代表第n个子区域和第q个子区域分界处的第m个参考点的二维坐标，\left( \Delta {x_{m,nq}}, \Delta {y_{m,nq}} \right) 代表第n个子区域和第q个子区域分界处的第m个参考点的二维配准偏移量。{a_{n1}},{b_{n1}},{c_{n1}},{a_{n2}}, {b_{n2}},{c_{n2}}, {a_{q1}},{b_{q1}},{c_{q1}},{a_{q2}},{b_{q2}},{c_{q2}} 分别为待估计参数。参考点选择时，在子区域内选择主图像中信噪比较高的点；同时，选择所有子区域边界处的点作为参考点。
如图10所示，以主图像的参考点
\left( {i,j} \right) 为中心点，确定一定大小的匹配窗口，在辅图像上相同位置为中心选取比匹配窗口大的搜索窗口。沿着行向和列向，在搜索窗中顺序移动匹配窗，同时计算每个位置处的相干系数，选取相干系数最大的位置作为同名点的位置，并得到同名点在主辅图像中的二维偏移量。根据式(20)和式(21)，构建目标函数联合求解，目标函数可表示为
\left\{ \begin{aligned} & \mathop {\min }\limits_{{a_{n1}},{b_{n1}},{c_{n1}}} \sum\limits_{n = 1}^N \sum\limits_{i = 1}^M \bigr( {a_{n1}}{x_{m,n}} + {b_{n1}}{y_{m,n}} \\ & \qquad + {c_{n1}} - \Delta {x_{m,n}} \bigr)^2 \\ & \mathop {\min }\limits_{{a_{n2}},{b_{n2}},{c_{n2}}} \sum\limits_{n = 1}^N \sum\limits_{i = 1}^M \bigr( {a_{n2}}{x_{m,n}} + {b_{n2}}{y_{m,n}} \\ & \qquad + {c_{n2}} - \Delta {y_{m,n}} \bigr)^2 \\ & {\mathrm{s.t}}. \\ & \quad {g_1} = {a_{n1}}{x_{m,nq}} + {b_{n1}}{y_{m,nq}} + {c_{n1}} - {a_{q1}}{x_{m,nq}}\\ & \qquad\quad - {b_{q1}}{y_{m,nq}} - {c_{q1}} = 0 \\ & \quad {g_2} = {a_{n2}}{x_{m,nq}} + {b_{n2}}{y_{m,nq}} + {c_{n2}} - {a_{q2}}{x_{m,nq}} \\ & \qquad\quad - {b_{q2}}{y_{m,nq}} - {c_{q2}} = 0 \end{aligned} \right. (22) 引入Lagrange算子，求解目标函数。引入Lagrange算子后的目标函数可表示为
\left\{ \begin{aligned} & \mathop {\min }\limits_{{a_{n1}},{b_{n1}},{c_{n1}}} \sum\limits_{n = 1}^N \sum\limits_{i = 1}^M \bigr( {a_{n1}}{x_{m,n}} + {b_{n1}}{y_{m,n}} + {c_{n1}} \\ & \qquad - \Delta {x_{m,n}} \bigr)^2 - \sum\limits_{n = 1}^N {{\lambda _n}{g_1}} \\ &\mathop {\min }\limits_{{a_{n2}},{b_{n2}},{c_{n2}}} \sum\limits_{n = 1}^N \sum\limits_{i = 1}^M \bigr( {a_{n2}}{x_{m,n}} + {b_{n2}}{y_{m,n}} + {c_{n2}} \\ & \qquad - \Delta {y_{m,n}} \bigr)^2 - \sum\limits_{n = 1}^N {{\mu _n}{g_2}} \end{aligned} \right. (23) 对于式(23)，可利用最小二乘法联合求解模型参数。根据图像中所有点的二维坐标，即可以实现全局偏移量的估计。最后，对辅图像进行插值，即可获得精配准后的图像。图11所示为基于地形辅助的无人机载InSAR图像分区配准方法流程图。
图 11 基于地形辅助的无人机载InSAR图像分区配准方法流程图Figure 11. Flowchart of topography-assisted UAV InSAR registration method with image partition4. 实测数据验证
4.1 P波段无人机载SAR系统
为验证本文提出的基于图像分区的无人机载InSAR图像配准方法的有效性，采用北京理工大学自研的P波段无人机载SAR原理样机开展验证实验。系统照片如图12所示，系统参数如表1所示。
表 1 系统参数Table 1. System parameters参数 数值 中心频率 400 MHz 信号带宽 60 MHz 脉冲宽度 2 μs 脉冲重复周期 200 μs 4.2 实测数据分析
4.2.1 北京平谷金海湖机场实验
4.2.1.1 实验信息
利用P波段无人机载SAR在北京平谷区金海湖机场(40.19°N, 117.23°E)附近进行了无人机载SAR干涉测量实验。图13所示为实验场景，蓝色实线表示无人机飞行轨迹，红色实线圈出的区域为雷达观测区域。观测区域主要由一些山丘组成。航迹总长480 m，无人机平台飞行速度为8 m/s，单航过飞行时长为60 s。
两次航过的航迹情况如图14所示，可以看到水平向和垂直向偏离直线运动的分量很小，最大运动误差不超过0.2 m。图15为插值前后的外部数字高程模型(Digital Elevation Model, DEM)。原始DEM来自SRTM官方数据，空间分辨率为90 m，需将DEM结果插值到成像像素分辨单元下，这里插值后空间分辨率为1 m。结合图15，其相位误差分量满足式(10)所示条件，成像结果可以聚焦。
图16为金海湖场景中P波段无人机载SAR成像结果，其中以最大幅值为参考对图像进行了归一化处理。利用高精度的位置测量系统，基于BP成像算法能够很好地补偿多普勒相位，实现高质量聚焦成像[29]。
4.2.1.2 配准结果分析
对于垂直基线为12 m的两个航过的雷达图像，获取的相干系数图与滤波后[30]的干涉相位图如图17所示。可以看出，高密度干涉条纹所在区域的相干性较低。这些严重失相干区域的高程较大，即有地形的区域。
图 17 金海湖场景12 m基线配准前相干系数图和干涉相位图Figure 17. Coherence coefficient diagram and interferogram with 12 m baseline before registration of Jinhai Lake结合两个航过的航迹信息，根据式(19)，设置偏移量门限为0.5个像素单元，可以得到如图18(a)所示的门限结果。可以看到，距离向越近的目标，斜距越短，对高程也越敏感，在很低的高程即会出现较大的配准偏移量。图18(b)表示偏移量大于1.5个像素单元时所需的高程门限，可以看到此时的二级高程门限已经高于图15所示的DEM结果，此时无需使用二级高程门限。以高程门限将所示的外部地形结果进行分区，并进行连通域处理，剔除面积较小的部分，得到图19所示的图像分区结果，可以看到，地形区域均有效地被提取出来。
图20所示为图像分区结果与相关系数、干涉相位的对比结果。可以看出，由于地形导致的图像失配，被提取出来的地形区域普遍存在相干性低、干涉相位质量差的问题。
图 20 金海湖场景地形分区结果Figure 20. Comparison between the segmentation results of Jinhai Lake图像分区后，将各区域内高信噪比和区域边界处的像素点作为参考点，利用基于窗口滑动的配准偏移量提取方法对所有参考点的二维配准偏移量进行提取，并进一步筛选出高相干性的参考点。图21为精配准结果，可以看到在地形起伏区域干涉相位质量大幅提升。
图 21 金海湖场景本文所提方法的精配准结果Figure 21. Fine registration results of the proposed method of Jinhai Lake4.2.1.3 比较与评估
为了进一步说明本文所提方法的有效性，对比了几种常见的InSAR配准方法[31]，分别是平均波动函数、最大频谱法和最大互相关法。利用3种方法配准后的相干系数图和干涉相位图如图22所示。
图 22 金海湖场景常见InSAR方法配准后的相干系数图和干涉相位图Figure 22. Coherence coefficient diagrams and interferograms after registration of common InSAR methods of Jinhai Lake从图22可以看出，常见的InSAR配准方法虽然能够提升一定的相干系数和干涉图质量，但对于起伏地形区域提升有限。相比常见的InSAR配准方法，本文所提的配准方法，以地形分区为先导，对地形分区后的不同高程子区域进行分区配准，保证了各子区域的干涉图质量。如图23所示，本文所提方法较常规方法有一定程度的提升，平均相干系数与残差点数量如表2所示，可以看到残差点数量大幅减少，证明了本文所提方法的有效性。图24为所提方法反演出的地形。
图 23 金海湖场景不同方法相干系数分布曲线Figure 23. Coherence coefficient distribution with different methods of Jinhai Lake表 2 金海湖场景干涉图评价Table 2. Evaluation of interferograms of Jinhai Lake方法 平均相干系数 残差点占比(%) 配准前 0.44 23.60 AFF 0.48 19.88 MSF 0.47 20.35 MCC 0.52 19.66 所提方法 0.58 15.33 4.2.2 重庆奉节老林沟实验
4.2.2.1 实验信息
图25所示为重庆奉节县老林沟(30.89°N, 109.48°E)实验场景，蓝色实线表示无人机飞行轨迹，红色实线圈出的区域为雷达观测区域。观测区域主要由一舌头形山体组成，如黄色虚线所示。航迹总长480 m，无人机平台飞行速度为8 m/s，单航过飞行时长为60 s。
两次航过的航迹情况如图26所示，可以看到水平向和垂直向偏离直线运动的分量同样很小，最大运动误差不超过0.2 m。如图27所示为插值前后的 DEM，插值后空间分辨率仍为1 m。结合图27所示的DEM，其相位误差分量满足式(10)所示条件，成像结果可以聚焦。图28为P波段无人机载SAR成像结果，同样以最大幅值为参考对图像进行了归一化处理。从图28可以看到，图像的左边区域有部分阴影区域，信噪比较低。
4.2.2.2 配准结果分析
对于垂直基线为13 m的两个航过的雷达图像，获取的相干系数图与干涉相位图如图29所示。可以看出，由于地形相对较高，除了接近平地的区域有一定的相干性和干涉条纹外，其余区域内几乎完全失相干，看不到干涉条纹。
图 29 老林沟场景13 m基线配准前相干系数图和干涉相位图Figure 29. Coherence coefficient diagram and interferogram with 13 m baseline before registration of Laolin Gou结合两个航过的航迹信息，根据式(19)，设置偏移量门限为0.5个像素单元，可以得到如图30(a)所示的门限结果。图30(b)表示偏移量大于1.5个像素单元时所需的高程门限，可以看到此时的二级高程门限已经包含了绝大部分的DEM结果，无需构建三级高程门限。以高程门限将所示的外部地形结果进行分区，并进行连通域处理剔除面积较小的部分，得到图31所示的图像分区结果。对照图29所示的干涉相位图和相干系数图，可以看出低于一级门限的区域尚能保持一定的相干性并有一定程度的干涉条纹，而高于一级门限区域几乎失相干，说明地形区域均有效地被提取出来。
图像分区后，利用所提方法进行图像配准。图32所示为精配准结果，可以看到在地形起伏区域干涉相位质量大幅提升。
图 32 老林沟场景本文所提方法的精配准结果Figure 32. Fine registration results of the proposed method of Laolin Gou4.2.2.3 比较与评估
同样利用4.2.1.3节所述的3种常见InSAR配准方法进行对比，得到的相干系数图和干涉相位图如图33所示。
图 33 老林沟场景常见InSAR方法配准后的相干系数图和干涉相位图Figure 33. Coherence coefficient diagrams and interferograms after registration of common InSAR methods of Laolin Gou从图33可以看出，在整体信噪比和相干系数较低的情况下，AFF和MSF方法受到噪声的影响严重，灵敏度大幅降低；MCC方法虽然能够提升一定的相干系数和干涉图质量，但对于起伏地形区域提升有限。相比常规方法，本文所提的配准方法更好地提升了干涉图质量。如图34所示，本文所提方法较常规方法有一定程度的提升，平均相干系数与残差点数量如表3所示，可以看到残差点数量大幅减少，证明了本文所提方法的先进性。
图 34 老林沟场景不同方法相干系数分布曲线Figure 34. Coherence coefficient distribution with different methods of Laolin Gou表 3 老林沟场景干涉图评价Table 3. Evaluation of interferograms of Laolin Gou方法 平均相干系数 残差点占比(%) 配准前 0.30 35 AFF 0.31 35 MSF 0.32 34 MCC 0.35 32 所提方法 0.38 28 5. 结语
本文面向无人机平台干涉测量，重点研究了无人机平台本身给图像配准带来的影响和解决方案。对于复杂地形，无人机载InSAR不同同名像素之间的偏移量大且具有明显的空变特性，常规的多项式拟合的变换方法难以取得良好的全局配准效果。针对这一问题，本文提出了一种基于图像分区的无人机载InSAR图像配准方法，利用外部粗地形，并结合航迹信息对图像进行分区处理，将全图的偏移量拟合问题，转化为若干个子区域的偏移量在约束条件下的联合拟合求解问题，实现图像对之间的精配准。
该方法首先基于航迹信息生成高程门限，再利用外部粗地形对测量区域进行图像分区处理。之后对各子区域内的偏移量进行多项式拟合，并对区域之间边界处的偏移量施加约束条件，构建带约束的多项式拟合目标函数，引入Lagrange算子，利用最小二乘法联合求解模型参数，并进行图像重采样获取精配准后的图像。基于P波段无人机载InSAR获取的实测数据，验证了该方法的有效性。
本文还有部分问题尚未解决。首先，本文所用数据集十分有限，无法保证所提方法能够适应不同复杂地形条件。其次，我们暂未获取到对应场景的较为准确的DEM，因此未对生成地形的精度进行评估。最后，对于所提方法自身，由于实验条件所限，也未对配准精度直接进行定量评估。后续，课题组将继续致力于无人机载InSAR测量领域的关键技术研究，获取尽可能多的地形场景验证所提方法的鲁棒性和适应性。
-
图 5 迭代解与真实投影点位置误差
Figure 5. Error between the iteration solution and the real projection point
图 8 不同最大运动误差和目标地距下最大斜距误差与目标高程之间的关系
Figure 8. Relationship between the maximum slant range error and the target height with different the maximum motion error and the target ground distance
图 11 基于地形辅助的无人机载InSAR图像分区配准方法流程图
Figure 11. Flowchart of topography-assisted UAV InSAR registration method with image partition
图 17 金海湖场景12 m基线配准前相干系数图和干涉相位图
Figure 17. Coherence coefficient diagram and interferogram with 12 m baseline before registration of Jinhai Lake
图 20 金海湖场景地形分区结果
Figure 20. Comparison between the segmentation results of Jinhai Lake
图 21 金海湖场景本文所提方法的精配准结果
Figure 21. Fine registration results of the proposed method of Jinhai Lake
图 22 金海湖场景常见InSAR方法配准后的相干系数图和干涉相位图
Figure 22. Coherence coefficient diagrams and interferograms after registration of common InSAR methods of Jinhai Lake
图 23 金海湖场景不同方法相干系数分布曲线
Figure 23. Coherence coefficient distribution with different methods of Jinhai Lake
图 29 老林沟场景13 m基线配准前相干系数图和干涉相位图
Figure 29. Coherence coefficient diagram and interferogram with 13 m baseline before registration of Laolin Gou
图 32 老林沟场景本文所提方法的精配准结果
Figure 32. Fine registration results of the proposed method of Laolin Gou
图 33 老林沟场景常见InSAR方法配准后的相干系数图和干涉相位图
Figure 33. Coherence coefficient diagrams and interferograms after registration of common InSAR methods of Laolin Gou
图 34 老林沟场景不同方法相干系数分布曲线
Figure 34. Coherence coefficient distribution with different methods of Laolin Gou
表 1 系统参数
Table 1. System parameters
参数 数值 中心频率 400 MHz 信号带宽 60 MHz 脉冲宽度 2 μs 脉冲重复周期 200 μs 
下载: 导出CSV
表 2 金海湖场景干涉图评价
Table 2. Evaluation of interferograms of Jinhai Lake
方法 平均相干系数 残差点占比(%) 配准前 0.44 23.60 AFF 0.48 19.88 MSF 0.47 20.35 MCC 0.52 19.66 所提方法 0.58 15.33
下载: 导出CSV
表 3 老林沟场景干涉图评价
Table 3. Evaluation of interferograms of Laolin Gou
方法 平均相干系数 残差点占比(%) 配准前 0.30 35 AFF 0.31 35 MSF 0.32 34 MCC 0.35 32 所提方法 0.38 28
下载: 导出CSV
-
[1] MOREIRA A, PRATS-IRAOLA P, YOUNIS M, et al. A tutorial on synthetic aperture radar[J]. IEEE Geoscience and Remote Sensing Magazine, 2013, 1(1): 6–43. doi: 10.1109/MGRS.2013.2248301 [2] 韦顺军, 唐欣欣, 张晓玲. 基于DFT模型的大场景InSAR图像配准[J]. 遥感学报, 2019, 23(5): 859–870. doi: 10.11834/jrs.20197459WEI Shunjun, TANG Xinxin, and ZHANG Xiaoling. Image registration algorithm for InSAR large scenes via DFT model[J]. Journal of Remote Sensing, 2019, 23(5): 859–870. doi: 10.11834/jrs.20197459 [3] LIN Xue, FANGI D S, and LI Fangfang. Acqusition of high-precision digital terrain model using P-band airborne repeat-pass SAR interferometry[C]. 2016 CIE International Conference on Radar (RADAR), Guangzhou, China, 2016: 1–4. [4] ZENG Tao, LIU Minkun, WANG Yan, et al. Tomographic SAR imaging with large elevation aperture: A P-band small UAV demonstration[J]. Science China Information Sciences, 2022, 65(3): 132303. doi: 10.1007/s11432-021-3391-2 [5] 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 [6] ZOU Weibao and CHEN Libin. Determination of optimum tie point interval for SAR image coregistration by decomposing autocorrelation coefficient[J]. IEEE Transactions on Geoscience and Remote Sensing, 2019, 57(7): 5067–5084. doi: 10.1109/TGRS.2019.2896383 [7] MA Zhangfeng, JIANG Mi, ZHAO Yi, et al. Minimum spanning tree co-registration approach for time-series sentinel-1 TOPS data[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2019, 12(8): 3004–3013. doi: 10.1109/JSTARS.2019.2920717 [8] YAGUE-MARTINEZ N, EINEDER M, CONG Xiaoying, et al. Ground displacement measurement by TerraSAR-X image correlation: The 2011 Tohoku-Oki earthquake[J]. IEEE Geoscience and Remote Sensing Letters, 2012, 9(4): 539–543. doi: 10.1109/LGRS.2012.2196020 [9] SCHEIBER R and 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] ZOU Weibao, LI Yan, LI Zhilin, et al. Improvement of the accuracy of InSAR image co-registration based on tie points— a review[J]. Sensors, 2009, 9(2): 1259–1281. doi: 10.3390/s90201259 [11] ZITOVÁ B and FLUSSER J. Image registration methods: A survey[J]. Image and Vision Computing, 2003, 21(11): 977–1000. doi: 10.1016/S0262-8856(03)00137-9 [12] HE Ke, CHEN Yongguang, JIA Xin, et al. A method based on interpolation mapping function for InSAR image registration[C]. 2016 IEEE International Conference on Electronic Information and Communication Technology (ICEICT), Harbin, China, 2016: 351–355. [13] DELLINGER F, DELON J, GOUSSEAU Y, et al. SAR-SIFT: A SIFT-Like algorithm for SAR images[J]. IEEE Transactions on Geoscience and Remote Sensing, 2015, 53(1): 453–466. doi: 10.1109/TGRS.2014.2323552 [14] PAUL S and PATI U C. SAR image registration using an improved SAR-SIFT algorithm and Delaunay- Triangulation-based local matching[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2019, 12(8): 2958–2966. doi: 10.1109/JSTARS.2019.2918211 [15] WANG Linhui, XIANG Yuming, YOU Hongjia, et al. A robust multiscale edge detection method for accurate SAR image registration[J]. IEEE Geoscience and Remote Sensing Letters, 2023, 20: 4006305. doi: 10.1109/LGRS.2023.3279141 [16] LI Xin, WANG Taoyang, CUI Hao, et al. SARPointNet: An automated feature learning framework for spaceborne SAR image registration[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2022, 15: 6371–6381. doi: 10.1109/JSTARS.2022.3196383 [17] YANG Zhuoqian, DAN Tingting, and YANG Yang. Multi-Temporal remote sensing image registration using deep convolutional features[J]. IEEE Access, 2018, 6: 38544–38555. doi: 10.1109/ACCESS.2018.2853100 [18] XIANG Yuming, PENG Lingxiao, WANG Feng, et al. Fast registration of multiview slant-range SAR images[J]. IEEE Geoscience and Remote Sensing Letters, 2022, 19: 4007505. doi: 10.1109/LGRS.2020.3045099 [19] YE Yuanxin, TANG Tengfeng, ZHU Bai, et al. A multiscale framework with unsupervised learning for remote sensing image registration[J]. IEEE Transactions on Geoscience and Remote Sensing, 2022, 60: 5622215. doi: 10.1109/TGRS.2022.3167644 [20] BOVEIRI H R, KHAYAMI R, JAVIDAN R, et al. Medical image registration using Deep Neural Networks: A comprehensive review[J]. Computers & Electrical Engineering, 2020, 87: 106767. doi: 10.1016/j.compeleceng.2020.106767 [21] FREY O, WERNER C L, and COSCIONE R. Car-borne and UAV-borne mobile mapping of surface displacements with a compact repeat-pass interferometric SAR system at L-band[C]. IGARSS 2019-2019 IEEE International Geoscience and Remote Sensing Symposium, Yokohama, Japan, 2019: 274–277. [22] WANG Yan, DING Zegang, LI Linghao, et al. First demonstration of single-pass distributed SAR tomographic imaging with a P-band UAV SAR prototype[J]. IEEE Transactions on Geoscience and Remote Sensing, 2022, 60: 1–18. doi: 10.1109/TGRS.2022.3221859 [23] WANG Zhen, DING Zegang, SUN Tao, et al. UAV-based P-band SAR tomography with long baseline: A multimaster approach[J]. IEEE Transactions on Geoscience and Remote Sensing, 2023, 61: 1–21. doi: 10.1109/TGRS.2023.3272039 [24] 邓云开, 袁泉, 胡程, 等. 一种多部地基SAR联合观测时的图形配准方法[J]. 信号处理, 2018, 34(11): 1269–1276. doi: 10.16798/j.issn.1003-0530.2018.11.001DENG Yunkai, YUAN Quan, HU CHENG, et al. An image registration method applied for joint measurements of multiple GB-SAR[J]. Journal of Signal Processing, 2018, 34(11): 1269–1276. doi: 10.16798/j.issn.1003-0530.2018.11.001 [25] NITTI D O, HANSSEN R F, REFICE A, et al. Impact of DEM-assisted coregistration on high-resolution SAR interferometry[J]. IEEE Transactions on Geoscience and Remote Sensing, 2011, 49(3): 1127–1143. doi: 10.1109/TGRS.2010.2074204 [26] 黄晓涛, 陈乐平, 范崇祎, 等. 低频叶簇穿透雷达成像技术[J]. 电波科学学报, 2020, 35(4): 469–485. doi: 10.13443/j.cjors.2020041302HUANG Xiaotao, CHEN Leping, FAN Chongyi, et al. Low frequency foliage penetration radar imaging technology[J]. Chinese Journal of Radio Science, 2020, 35(4): 469–485. doi: 10.13443/j.cjors.2020041302 [27] HAN Congrong, TIAN Weiming, and MEI Hongyan. MIMO radar fast imaging algorithm based on sub-image combination[C]. 2019 6th Asia-Pacific Conference on Synthetic Aperture Radar (APSAR), Xiamen, China, 2019: 1–6. [28] LIN Jiahe, LV Xiaolei, and LI Rui. Automatic registered back-projection approach based on object orientation for airborne repeat-track interferometric SAR[J]. IET Radar, Sonar & Navigation, 2018, 12(9): 1066–1076. doi: 10.1049/iet-rsn.2018.5053 [29] 田卫明, 刘富强, 谢鑫, 等. 基于GPU粗细粒度和混合精度的SAR后向投影算法的并行加速研究[J/OL]. 信号处理, 1–12. http://kns.cnki.net/kcms/detail/11.2406.tn.20230509.0955.002.html, 2023.TIAN Weiming, LIU Fuqiang, XIE Xin, et al. Research on parallel acceleration processing technology of SAR back projection algorithm based on two granularities and mixing precision of GPU[J/OL]. Journal of Signal Processing, 1–12. http://kns.cnki.net/kcms/detail/11.2406.tn.20230509.0955.002.html, 2023. [30] 冯丽源, 邓云开, 聂祥飞, 等. 一种基于Non-Local Means的地基差分干涉雷达相位滤波改进方法[J]. 信号处理, 2022, 38(1): 100–108. doi: 10.16798/j.issn.1003-0530.2022.01.012FENG Liyuan, DENG Yunkai, NIE Xiangfei, et al. An improved phase filtering method for Ground-based differential interferometer radar based on Non-Local means[J]. Journal of Signal Processing, 2022, 38(1): 100–108. doi: 10.16798/j.issn.1003-0530.2022.01.012 [31] FANG Dongsheng, LV Xiaolei, YUN Ye, et al. An InSAR fine registration algorithm using uniform tie points based on voronoi diagram[J]. IEEE Geoscience and Remote Sensing Letters, 2017, 14(8): 1403–1407. doi: 10.1109/LGRS.2017.2715189 期刊类型引用(0)
其他类型引用(1)
-
下载:
计量
- 文章访问数: 756
- HTML全文浏览量: 246
- PDF下载量: 196
- 被引次数: 1
