一种多部地基SAR联合观测时的图形配准方法

邓云开 袁 泉 胡 程 田卫明 赵 政

(1. 北京理工大学信息与电子学院雷达技术研究所,北京 100081； 2. 湖南省岳阳市气象局,湖南岳阳 414000;3. 卫星导航电子信息技术教育部重点实验室(北京理工大学),北京 100081)

1 引言

滑坡灾害是地质灾害中发生频率最高和危害最大的一种,给人类的生命和财产安全带来巨大的威胁[1]。在边坡宏观失稳之前,其表面都会先发生微小位移。因此,获取监测区域的表面形变信息,对于滑坡灾害的预测预警具有很重要的意义。近二十年来,地基SAR(Synthetic Aperture Radar,合成孔径雷达)在形变监测领域得到了广泛的应用,如对矿区边坡、陡峭山体、冰川、水电站大坝和危险建筑物等的滑坡监测[2]。作为一种高精度的形变测量仪器,地基SAR具有全天时、全天候、大范围覆盖、非接触式测量、高精度和形变测量速度快等优点。

在实际应用中,地基SAR系统的一个典型缺点是,仅能获取雷达视线方向的一维形变信息[3]。如采用地基SAR系统对矿区边坡、陡峭山体等进行形变测量时,由于不同类型的坡体会表现出不同的表面形变模式,且同一坡体的不同区域的表面形变特征也不同。仅采用一部地基SAR获取一维视线方向的形变信息,可能与监测区域的真实形变方向存在较大的偏差[4]。因此,地基SAR形变测量技术还存在较大的改进空间。如果可以获取到可能会发生滑坡区域的三维形变信息,结合相应的力学机理及地质模型,有利于进行预报预警。

在地基SAR领域,为实现三维形变信息的测量,最基本的解决方案是采用至少三部地基SAR从不同位置联合观测,并获取目标区域在不同方向上的一维形变信息,然后实现三维形变信息的解算[5]。但在地基SAR成像时,获取的是观测区域的三维地形在雷达成像平面上的投影,由于几何畸变的影响,不同观测角度下的地基SAR图像差别较大,基于常规的多项式拟合的变换方法无法实现图像配准。因此,需要提出一种新的适合多部地基SAR联合观测下的图像配准方法。

针对上述问题,本文提出了一种几何配准方案。首先获取观测场景的DEM(Digital Elevation Model,数字高程模型)信息以及各部雷达的位置信息,并将场景高程信息及各部雷达统一到同一坐标系,其次基于雷达成像的投影原理,实现各部雷达的图像中同名投影点的关联。

2 雷达图像配准

2.1 常规配准方案

在星载干涉SAR领域,由于获取两幅SAR图像的观测几何存在差异,为保证用于干涉计算的两幅图像像素对应于地面同一散射单元,需要对图像进行配准处理。一般情况下,图像配准的方案是基于某一种配准测度准则,计算两幅SAR图像在一些控制点处的匹配度量值,然后建立多项式模型来求解出其他像素点处的配准偏移量,从而实现像素级或者亚像素级的配准[6]。

(1)

式(1)为图像配准中常用的二阶多项式模型,其中,a0～a5、b0～b5为拟合参数,(x1,y1)、(x2,y2)分别为同名控制点在主辅图像中的坐标。如果可以确定多对同名控制点的坐标,则可以根据式(1)采用最小二乘法求解出拟合系数。

图1 多部雷达联合观测示意图Fig.1 Diagram of joint measurements of multiple GB-SAR

在采用多部地基SAR联合观测进行三维形变测量时,为获取较高的三维形变测量精度,各部雷达系统的观测视角应存在较大的差异,如图1所示。但在地基SAR成像时,获取的是观测区域的三维地形在雷达成像平面上的投影[7],由于几何畸变的影响,不同观测角度下的地基SAR图像差别较大,需要分析基于常规的多项式拟合的配准方法是否依然适用。

如果参考坐标系中存在一点P(xp,yp,zp),首先基于地基SAR成像的投影原理,分析该点P在这两个雷达坐标系下的投影点P1(xp1,yp1)和P2(xp2,yp2)。图2所示为地基SAR成像的几何示意图,雷达1的合成孔径方向沿x轴,孔径中心位于坐标原点,雷达1的成像平面为x-y平面。在雷达1坐标系下,投影点P1的坐标可以表示为

(2)

图2 地基SAR成像几何示意图Fig.2 Diagram of GB-SAR imaging

两组单位矢量之间的关系可以表示为

(3)

其中,ux=cosθcosφ,uy=sinθcosφ,uz=sinφ,νx=-sinθ,νy=cosθ,νz=0,nx=-cosθsinφ,ny=-sinθsinφ,nz=cosφ。

两个雷达坐标系之间的变换关系可以表示为[9]

(4)

其中,(x,y,z)和(x′,y′,z′)分别为一点在两个雷达坐标系下的坐标。则点P在雷达2坐标系下的三维坐标可以表示为

(5)

因此,在雷达2坐标系下,投影点P2的坐标可以表示为

(6)

其中,xp2为投影点P2的方位向坐标,为点P到雷达2的斜距,yp2为距离向坐标。

结合式(2)和式(6),可以看出,点P在两个雷达坐标系下的投影点坐标(xp1,yp1)和(xp2,yp2)之间无法建立起式(1)所示的多项式转换关系。这也就意味着,基于常规的多项式拟合的变换方法无法实现多部地基SAR联合观测下的图像配准。

2.2 几何配准方案

为解决多部(M部)地基SAR联合观测进行三维形变测量时存在的图像配准问题,本文提出了一种几何配准的方案。关键问题包括两个:三维坐标系变换和同名像素点关联。

首先采用激光雷达获取观测场景的DEM信息,然后采用差分GPS获取各部雷达的坐标信息。由于激光雷达和差分GPS分别采用各自的坐标系,因此需要解决三维坐标系变换的问题,将场景DEM和各部雷达统一到同一坐标系下。

将场景DEM及各部雷达统一到同一坐标系后,可以基于雷达成像的投影原理,将表征DEM的每一个地形点在各部雷达的成像视角下进行投影。为实现各部雷达的图像的配准,需要研究同名像素点的关联方法。

2.2.1 三维坐标系变换

为确定两个三维直角坐标系之间的变换关系,需要已知至少三个控制点在这两个坐标系下的三维坐标[10]。在采用激光雷达获取场景DEM时,布设多个棱镜作为控制点,测量这些棱镜的三维坐标,之后采用差分GPS同样获取这些棱镜的三维坐标。

假设一个控制点在两个三维直角坐标系下的坐标分别为(xs,ys,zs)T和(xt,yt,zt)T,可以建立如下变换关系:

(7)

其中,(Δx,Δy,Δz)T为两个坐标系之间的平移参数,k为尺度因子,α、β和γ为三个旋转参数,Rx(α)、Ry(β)和Rz(γ)为三个旋转矩阵,且

建立误差方程V=f(X)-L,其中X为需要估计的参数向量,L=(xt,yt,zt)T为观测向量,f(X)为由未知参数X组成的非线性函数。

(8)

2.2.2 同名投影点关联

由于雷达成像结果中部分像素点对应的是非有效回波区域,可以先对雷达图像中的所有像素点进行筛选,如利用多幅雷达图像,可以基于幅度离差法选择出一些幅度高度稳定的像素点[11],定义为像素点子集Si,i=1,2,...,M。每个子集Si中的像素点表示为Kn,i,n=1,2,...,Ni,Ni为第i个子集中的像素点数量。

(9)

其中,Δx和Δy分别表示各部雷达成像时的x方向和y方向的取样间隔。

3 雷达图像配准

为验证本文提出的方法的有效性,基于北京理工大学自主研发的两部SAR系统,对河北省迁安市马兰庄的一处露天铁矿进行了监测。图3(a)所示为监测场景的照片,红色和绿色五角星分别代表两部雷达的摆放位置。图3(b)所示为采用激光扫描仪获取的雷达监测区域的DEM结果,方位向和距离向上采样间隔均为1 m。该矿坑具有典型的边坡结构,矿坑直径约500 m,矿坑周围存在山体,符合大场景成像的需求。

图3 露天矿坑Fig.3 Open-pit mine

图4 雷达照片Fig.4 Photos of the radar systems

雷达1为一部轨道SAR,基于天线沿着高精密滑轨的机械移动来形成合成孔径,如图4(a)所示[12]。雷达2为一部MIMO-SAR(Multiple-Input Multiple-Output SAR),采用多输入多输出技术,通过多个发射天线和接收天线的特殊排列来等效成一个大的合成孔径,可以视为一种特殊工作体制的地基SAR,如图4(b)所示[13]。

表1所示为两部雷达系统的参数,除了雷达1的合成孔径长度为1.2 m,雷达2的合成孔径长度为1.138 m,以及相应的方位角分辨率有所不同外,其他参数均相同。相比于轨道SAR,MIMO-SAR的优势在于不需要采用高精密滑轨来实现合成孔径,可以避免天线在滑轨上运动时,由于天线的振动所带来的轨道相位误差。

表1 雷达系统参数

基于上述两部雷达系统,在2017年12月8日,对马兰庄矿坑进行了持续一小时的同步监测,分别获取了30幅雷达图像。在图像获取前,采用激光扫描仪获取观测场景的DEM信息,然后在其扫描范围内布设4个棱镜,并获取这4个棱镜在激光扫描仪坐标系下的三维坐标,之后再采用差分GPS对这4个棱镜及两部雷达进行三维坐标的测量。以这4个棱镜为控制点, 基于2.2.1小节的三维坐标变换方法,获取将差分GPS坐标系转换到激光扫描仪坐标系的转换矩阵,则可以计算出两部雷达在激光扫描仪坐标系下的位置信息和角度信息。计算可得,雷达1孔径中心三维坐标为(0.18 m,2.21 m,-0.27 m),方位角和俯仰角分别为θ1=0°和φ1=1.621°,雷达2孔径中心三维坐标为(-346.73 m,88.45 m,31.2 m),方位角和俯仰角分别为θ2=23.10°和φ2=-3.71°。

图5(a)和5(b)所示分别为这两部雷达的成像结果(dB图),由于两部雷达的观测视角差别较大,成像结果的几何特征很明显存在差别。在成像时,方位向及距离向的采样间隔均为1 m,即图像中像素点的大小为1 m×1 m。

图5 成像结果Fig.5 SAR imaging results

由于雷达成像结果中部分像素点对应的是非有效回波区域,尤其是对幅值很弱的像素点。基于幅度离差法,对这两部雷达的30幅图像分别进行了永久散射体的选择,选择门限设置为0.25。图6(a)和图6(b)所示分别为基于雷达1和雷达2的30幅图像选择出的永久散射体,对应的即为小节2.2.2中的像素点子集S1和S2。

基于小节2.2.2中的同名投影点关联方法即可以实现对这两部雷达的图像配准,图7所示为雷达图像配准结果。图中,红色的地形点表示其在两部雷达的观测视角下均为永久散射体,蓝色的地形点仅在雷达1的观测视角下为永久散射体,黄色的地形点则仅在雷达2的观测视角下为永久散射体。

图6 像素点子集Fig.6 Pixel subsets

图7 雷达图像配准结果Fig.7 Image registration result

4 结论

本文提出了一种新的适合多部地基SAR联合观测下的几何配准方案。通过一部轨道地基SAR和一部MIMO-SAR对一个露天矿坑从不同的位置进行了观测,首先基于多个棱镜实现了三维坐标系的变换,将两部雷达与场景地形信息统一到同一坐标系下,其次基于永久散射体方法分别选择出了两部雷达的图像中的像素点子集,并基于本文提出的同名投影点关联方法,实现了两部雷达的图像配准。文章基于实测数据初步验证了该方法的可行性,但还存在一些不足,如需要在观测场景中布设多个参考点验证配准精度、需要验证基于该配准方法实现多部雷达联合观测时三维形变测量的准确性等。