基于栅格投影的地基合成孔径雷达三维地形匹配俯仰角模糊处理方法

王彦平 刘航 李洋 林赟

（北方工业大学，北京 100144）

引 言

地基合成孔径雷达(ground-based synthetic aperture radar, GB-SAR)系统能够进行全天时、全天候的高分辨率成像[1]. 通过测量图像的相位，可以实现亚毫米级精度的形变反演，具有非接触、高精度、大区域连续监测的技术优势[2]. 作为一种主动式微波成像传感器，GB-SAR能在几分钟内获取数平方公里形变信息[3]. 通过在不同时间点对目标区域的重复观测，获取时间序列干涉数据用于形变监测，可获得精度为毫米级的形变监测结果，是进行区域性监测、地表形变监测以及定点连续测量的重要手段[4]，在露天矿监测[5-6]、山体滑坡监测[7-8]、大坝安全监测[9-10]等领域具有广阔的应用前景.

GB-SAR成像投影方式不同于传统地形图的正射投影以及摄影测量的中心投影，图像中的监测目标形态与实际有较大的差异，使用人员往往要花费较大精力进行图像的解读，造成使用上的不便. 三维激光扫描仪能够高速度、高精度和高密度地获取物体表面的三维空间坐标，且无需接触被测物体[11]，可以真实描述目标的整体结构及形态特性，并通过扫描测量点云来重现目标的完整原形及矢量化数据结构[12]. 将GB-SAR图像与三维激光扫描仪获取的三维地形进行匹配，可以将不直观的二维GB-SAR形变监测数据转变为更加贴近人眼视觉习惯的空间三维表达形式，方便GB-SAR图像解译人员迅速定位出形变区域.

近年来有不少研究者致力于GB-SAR三维地形匹配方法的研究. 文献[13]提出了一种数字高程模型(digital elevation model, DEM)数据与SAR图像匹配算法. 文献[14]提出一种利用三维激光扫描技术提供外部高程信息，依据GB-SAR成像几何投影原理，在投影面上完成点云和像元匹配的三维匹配方法. 文献[15]基于距离向和方位向条件的几何映射，提出了一种GB-SAR图像与地形数据三维匹配算法.

近距离对桥梁等高大建筑物的形变监测是GBSAR一个新的应用领域，对比常规的矿山或滑坡监测，其监测区域内不同目标之间常常存在较为特殊的几何关系. 在基于栅格投影的GB-SAR三维地形匹配的过程中，这种特殊的几何关系会造成俯仰角模糊现象，影响三维匹配结果的准确性，导致匹配结果的解译困难，因此需要研究一种处理俯仰角模糊现象的方法. 然而上述研究都没有针对近距离高大建筑物的监测，以及对基于栅格投影的GB-SAR三维匹配过程中俯仰角模糊现象的产生原因及处理方法进行讨论.

本文分析了基于栅格投影的GB-SAR三维地形匹配方法应用过程中俯仰角模糊现象的产生原因，并提出一种俯仰角模糊现象的处理方法. 实测结果表明，投影结果中虚假影像去除效果明显，说明本方法对俯仰角模糊现象具有良好的处理效果，有利于对GB-SAR三维地形匹配结果进行解译及分析.

1 GB-SAR距离-方位坐标系与东北天坐标系的转换关系

三维激光扫描仪能够快速获取高精度地面和地物三维信息，生成以东北天坐标系为基准的三维地形数据，而GB-SAR图像则是以其系统自身位置为原点建立距离-方位坐标系，要实现GB-SAR图像到三维地形的投影需保证两种数据处于同一坐标系.因此，本节分别建立东北天坐标系下的GB-SAR系统模型，以及距离-方位坐标系下GB-SAR平面几何模型，并对两种坐标的转换关系进行推导，实现坐标系的统一.

1.1 东北天坐标系下GB-SAR系统模型

相位中心点S的斜距为RSP，与雷达运动方向的夹角为θAP.x轴指向正北方，y轴指向正东方，z轴为高度向.

首先，使用三维激光扫描仪扫描监测场景，得到东北天坐标系下监测场景的三维地形. 东北天坐标系下GB-SAR系统三维模型如图1所示，使用实时动态(real-time kinematic，RTK)载波相位差分技术对GB-SAR系统导轨的两个端点进行坐标定位，取两端点坐标平均值作为雷达天线相位中心坐标，进一步可计算得到三维地形内各个目标点与雷达天线相位中心基于东北天坐标系的斜距以及与雷达运动方向的夹角. 如图1所示，三维地形中任意一点P与天线

图1 东北天坐标系下GB-SAR系统三维模型Fig. 1 3D model of GB-SAR system in east-north-up coordinate system

为方便表达，将图1投影到二维平面，建立东-北坐标系下GB-SAR系统二维模型，如图2所示. 原点O为测站，点S(xS,yS,zS)为雷达天线相位中心，直线AB表示GB-SAR系统导轨，点A(xA,yA,zA)、B(xB,yB,zB)为 导轨两个端点，点P(xP,yP,zP)为三维地形中任意一点，RSP为天线相位中心S到点P的距离，θSP为直线SP与y轴正方向的夹角，θSA为直线SA与y轴正方向夹角.

图2 东-北坐标系下GB-SAR系统二维模型Fig. 2 2D model of GB-SAR system in east-north coordinate system

在GB-SAR三维模型中，雷达天线相位中心S到点目标P的距离为

式中：xP、yP、zP和xS、yS、zS分别表示点P和点S的北坐标、东坐标及高度.

直线SP与y轴正方向的夹角θSP为

同理，直线SA与y轴正方向的夹角θSA为

则直线SP与SA夹角θAP为

1.2 距离-方位坐标系下GB-SAR平面几何模型

GB-SAR的图像是监测场景中雷达波束所覆盖的目标在雷达视线向的投影，各像元以雷达天线相位中心为原点，以到雷达天线相位中心的斜距为半径沿圆弧线投影到成像投影面上[16].

如图3所示，用距离-方位坐标系下GB-SAR平面几何模型(矩形abcd区域)表示GB-SAR图像，它是图1中的监测场景投影到斜距平面的斜距图像.点S为天线相位中心投影到斜距平面所对应的点.坐标轴x、y分别为雷达方位向和距离向，点P′为图1中点P投影到斜距平面所对应的点，点P′在GB-SAR图像中坐标为(RP,AP)，RP、AP可以表示为：

将式(4)代入式(5)和(6)中计算，可以将点P的东北天坐标转换到距离-方位坐标，这一过程称为GB-SAR距离-方位坐标系与激光三维地形东北天坐标系的转换.

图3 距离-方位坐标系下GB-SAR平面几何模型Fig. 3 GB-SAR plane geometry model in range-azimuth coordinate system

2 基于栅格投影的GB-SAR三维地形匹配方法

2.1 GB-SAR图像的栅格化处理

如图4所示为距离-方位坐标系下栅格化处理后的GB-SAR平面几何模型，Rmin和Rmax、Amin和Amax分别为雷达距离向和方位向上图像坐标的最小值与最大值. 图3中矩形abcd内的GB-SAR图像分别以分辨率dr、da为距离向和方位向采样间隔，且0

图4 距离-方位坐标系下栅格化处理后的GB-SAR平面几何模型Fig. 4 Gridded GB-SAR plane geometry model in rangeazimuth coordinate system

2.2 三维地形的栅格化处理

如图5所示为东-北坐标系下栅格化处理后的GB-SAR系统二维模型. 基于东北天坐标系，GBSAR系统以及所扫描三维地形处于二维平面xoy内，xoy为地球椭球面在原点处的切平面，x方向指向正北，y方向指向正东，xmin和xmax、ymin和ymax分别为三维地形北坐标和东坐标的最小值与最大值. 分别在x、y方向上以dx、dy为采样间隔建立三维地形栅格矩阵，且0

2.3 GB-SAR图像与三维地形的匹配

三维地形栅格矩阵中目标点P位于1行、列，且点P所在栅格带有点P高度信息zP. 故可根据前文内容，根据xP、yP、zP计算求得GB-SAR图像中与点P相匹配的点P′的距离-方位坐(RP,AP)，进一步可求得点P′在GB-SAR图像栅格矩阵中位于行、列，由此建立两个矩阵的联系. 将GB-SAR图像栅格矩阵中P点包含的GB-SAR图像信息，提取出并赋值在三维地形栅格矩阵中相应栅格. 将三维地形栅格矩阵中所有目标点按以上过程逐个进行投影，生成GB-SAR图像信息-高度信息对应关系矩阵，对应关系矩阵与三维地形栅格矩阵大小一致，且栅格排列顺序基于东-北坐标系. 提取对应关系矩阵中栅格的行列值、高度信息建立三维地形模型，然后提取栅格的GBSAR图像信息匹配到三维地形模型上，便实现了GB-SAR图像与三维地形的匹配.

图5 东-北坐标系下栅格化处理后的GB-SAR系统二维模型Fig. 5 2D model of gridded GB-SAR system in East-North coordinate system

3 基于栅格投影的三维地形匹配俯仰角模糊处理方法

3.1 俯仰角模糊现象

近距离对桥梁等高大建筑物进行形变监测时，其监测区域内不同目标之间常常存在较为特殊的几何关系. 在基于栅格投影的GB-SAR三维匹配的过程中，这种特殊的几何关系会造成俯仰角模糊现象.

图6所示为GB-SAR天线波束的方位向切面，M点和N点、O点和L点分别处于以雷达天线相位中心为中心点，以到雷达天线相位中心的斜距为半径的等距离弧线上，L点和M点所在区域存在高位散射体. 将四个点分别投影到地距上，理论上相对雷达的地距出现虚假投影值的栅格可能比含有真实投影值的栅格更远，如O点和L点. 但实际观测中往往不能穿透目标，因此出现虚假投影值的栅格一般比含有真实投影值的栅格在地距上更近，如M点与N点. 投影到地距后，近地距的N点、M点到雷达的斜距及与雷达形成的方位角相同，所投影的GBSAR图像信息来自GB-SAR图像中同一像元，近地距的N点出现与M点相似的模糊影像，此现象即俯仰角模糊现象. 很明显，N点所包含信息是虚假的，需要去除. 下面将对投影过程中产生的俯仰角模糊现象产生原因进行分析并提出相应的处理方法.

图6 GB-SAR天线波束的方位向切面Fig. 6 The azimuthal section of GB-SAR antenna beam

3.2 产生原因

在常规的GB-SAR矿山或滑坡监测中，GB-SAR面对边坡地形进行扫描，场景内监测目标的高度由近到远缓慢递增，很少存在斜距及方位角相同的监测目标. 将获取的GB-SAR图像与三维地形匹配，对应关系矩阵中不同栅格所带有的GB-SAR图像信息来自图像中的不同像元，不会造成投影值的误匹配，也不会产生俯仰角模糊现象. 而在近距离扫描如桥梁之类的高大建筑物的时候，容易出现特殊的几何关系，场景近地距地面投影栅格与远地距建筑物处的投影栅格有相同的斜距及方位角，进行三维匹配时，恰巧近距地面的栅格不具有有效激光雷达数据，导致GB-SAR图像中同一像元投影到不同栅格上，造成近地距地面投影栅格也虚假地出现了投影值，进而产生俯仰角模糊现象.

俯仰角模糊现象如图6所示. 雷达处于天线相位中心S点，M、N两点到雷达天线相位中心点S的斜距相同，且在同一方位角度上，但处于不同俯仰角上.

M、N分别为GB-SAR图像中(RM,AM) 、(RN,AN)两个像元的投影，由于M、N两点到雷达天线相位中心点S的距离相同，且在同一方位角度上，故θSM=θSN、RSM=RSN. 由式(5)和(6)可计算出，RM=RN、AM=AN.

经投影后，M、N两点由GB-SAR图像中同一像元投影得到，故所包含的图像信息(如反射强度、形变信息)也相同，造成N点俯仰角模糊现象.

3.3 修正方法

从以上分析可以看出，将如图6中点N这一类存在虚假投影值点进行处理，即可消除俯仰角模糊现象.

根据本文第2节内容，在东-北坐标系下，将三维地形进行栅格化处理，建立三维地形栅格矩阵，并进一步生成GB-SAR图像信息-高度信息对应关系矩阵.M点在对应关系矩阵中位于,，N点在对应关系矩阵中位于,，S点在对应关系矩阵中位于,，根据M、N、S三点在对应关系矩阵中的位置，提取三维地形栅格矩阵中相同位置栅格所带有的高度信息，M、S分别获取高度值zM、zS，三维坐标为 (xM,yM,zM)、(xS,yS,zS)；而N点无法索引有效激光点云数据，对应栅格高度值为λ，因此三维坐标为(xN,yN,λ). 由于M、N两点到天线相位中心点S的距离相同，即RSM=RSN，有

由式(7)可知，由于高度值λ的存在，导致N点到雷达天线相位中心S点的距离计算错误，进而产生错误的斜距、方位角，导致N点所包含的图像信息与M点相同，产生俯仰角模糊现象.

因此在对三维地形栅格矩阵索引的过程中，筛选所有高度值为λ的栅格，将不干扰匹配结果的GBSAR图像信息(如雷达反射强度信息可设置为0 dB，形变信息设置为0 mm)填入对应关系矩阵的对应栅格，即可完成修正. 根据前文内容可建立俯仰角模糊处理流程图，如图7所示.

以N点的处理为例，首先，在激光三维点云地形栅格矩阵中，N点位于矩阵处，xmin≤xN≤xmax、ymin≤yN≤ymax. 以m=1、n=1为起始逐点处理，分别以m、n为行列值到三维地形栅格矩阵中进行索引，提取对应栅格带有的高度值并反馈. 如果返回值不为λ，则索引GB-SAR图像对应像元的图像信息并填入对应关系矩阵；如果返回值为λ，则将不干扰匹配结果的GB-SAR图像信息填入对应栅格. 当处理N点时，，n=，在三维点云地形栅格矩阵索引高度值为λ，进行判断，将不干扰匹配结果的GB-SAR图像信息填入对应栅格，完成虚假投影值的去除. 逐点处理其余点，当且时处 理流程结束.

图7 俯仰角模糊处理流程图Fig. 7 Flowchart of the correction of pitch angle ambiguity

4 实测数据分析

新首钢大桥位于北京市境内，本试验使用GBSAR设备对在建中的大桥主体进行监测，雷达设备安装于距大桥中心约350 m处，雷达导轨安装稳定且基本处于水平状态，雷达监测现场如图8所示.

利用GB-SAR设备获取场景基于距离-方位坐标系的监测图像，结果如图9所示.

使用三维激光扫描仪获取GB-SAR监测区域的激光点云数据，数据基于东北天坐标系，结果如图10所示.

分别以dr=0.5 m、da=0.5 m为距离向和方位向采样间隔建立GB-SAR图像栅格矩阵，分别以dx=0.5 m 、dy=0.5 m为采样间隔建立三维地形栅格矩阵，然后进行GB-SAR图像与三维地形的匹配，建立GB-SAR图像信息-高度信息对应关系矩阵. 为方便展示俯仰角模糊现象，将对应关系矩阵的GBSAR图像信息显示，得到了东-北坐标系下的GBSAR图像，如图11(a)所示. 左上角红点为雷达位置，延伸出的红色箭头指向雷达视线方向. 图像中诸多区域存在虚假投影值，如图像中所标记红框内区域.将图9距离-方位坐标系下的原始GB-SAR图像进行局部放大和旋转，得到如图11(b)所示的距离-方位坐标系下的GB-SAR图像，与图11(a)相比，红框内目标差异较大.

图8 近距离桥梁监测现场Fig. 8 Close range of bridge monitoring

图9 距离-方位坐标系下的GB-SAR图像Fig. 9 GB-SAR images in range-azimuth coordinate system

图10 GB-SAR监测区域的激光点云数据Fig. 10 Data of laser point cloud at GB-SAR monitoring area

图11 两坐标系下GB-SAR图像对比(红框中的强散射点)Fig. 11 Comparison of GB-SAR images in 2 coordinate systems (the strong scatter points in red frame)

在图11(a)红框区域内，提取两个存在虚假投影值的点A、C，并分别与同一雷达视线方向上的B、D两点进行对比分析.A、C点处于近地距地面的投影栅格，B、D点处于远地距高位散射体所在的投影栅格. 经计算得知，A、C分别与B、D有相同的斜距及方位角，A、B两点同为距离-方位坐标系下GBSAR图像(315.5 m，− 67 m)处像元的投影，C、D两点同为(329 m，−74 m)处像元的投影. 另外，图10黄框内区域对应图11中红框区域，黄框区域激光点云数据缺失，因此与之对应的红框区域栅格不具有有效激光雷达数据. 根据前文内容可以判断A、C两点具有虚假投影值，存在俯仰角模糊现象，需要进行修正.

根据前文方法，对俯仰角模糊现象进行修正，处理结果如表1所示.B、D两点投影值为真，因此处理前后灰度值不变.A、C两点由于投影值是虚假的，将灰度值替换为0，使模糊影像变为黑色，从而消除对投影结果的影响. 修正后的结果如图12所示，可以看出A、C两点的虚假投影值被去除，且图像中诸多区域如图像中所标记红色区域内的模糊影像也得到了去除，与图11(b)所示的距离-方位坐标系下的GBS AR图像结果接近，效果理想.

表1 俯仰角模糊现象处理结果Tab. 1 Result of the correction of pitch angle ambiguity

图12 对俯仰角模糊现象进行修正后的图像Fig. 12 The image of corrected pitch angle ambiguity

5 结 论

本文针对近距离高大建筑物的形变监测中基于栅格投影的GB-SAR三维地形匹配出现的俯仰角模糊现象进行了分析，对俯仰角模糊现象产生原因进行了分析，通过判断三维地形栅格矩阵上每个栅格是否存在有效激光点云数据的方法，消除了投影结果中的虚假影像，对俯仰角模糊现象完成了修正. 选用北京新首钢大桥的实测数据对所提出的处理方法进行验证，投影结果中虚假影像去除效果明显，说明本方法对俯仰角模糊现象具有良好的处理效果，有利于对GB-SAR三维地形匹配结果进行解译及分析.