一种边坡雷达影像与地形点云的误匹配校正方法

本发明涉及一种边坡形变监测雷达影像与地形点云基于少量人工控制点的误匹配校正方法，属于露天矿边坡形变监测技术领域。

背景技术：

近年来，全世界地质灾害多发频发，各种地质原因和人类活动等因素均可能会诱发区域形变并衍生灾害。地基边坡形变监测雷达(gb-insar)技术越来越多地应用于露天矿边坡、水利水电工程、自然山体滑坡、崩塌、泥石流灾害监测场景。直线重复轨道雷达更是以其高时空分辨力和高形变测量精度，可连续大面积监测，受雨雪雾恶劣天气影响较小等优势，在各类边坡形变监测中得到日益广泛的应用。雷达在二维极坐标空间构像，对于不熟悉雷达成像原理的监测人员存在不直观的问题，实际使用需变换至熟悉的几何以直观定位形变区。

现有的技术方案用来实现雷达图像和地形点云匹配主要有平面相似变换法和几何映射方法。平面相似变换法为星载雷达方法在地基雷达领域应用的延伸，需要雷达图像精确的正射纠正，但常常难以获得纠正参数。几何映射通过测量轨道信息和成像参数，数据获取简单得到广泛应用。然而由于现场条件的特殊性，单独依靠几何映射方法匹配结果会发生较大偏差。目标边坡多在持续变形极度危险，使得人员难以深入布设大量地面控制点校正偏差，因此亟需依靠少量人工控制点校正几何映射匹配偏差。现已提出的雷达图像和地形点云匹配方法多是测量直线轨道端点然后依据成像参数进行几何映射的粗匹配方法，无法应对边坡监测准确定位异常形变区的目标。

技术实现要素：

为克服现有技术的不足，本发明提供一种基于少量人工目标控制、显著目标识别增选控制点和最小二乘优化的边坡形变监测雷达影像与地形点云误匹配校正方法。

本发明为解决上述技术问题采用以下技术方案：

一种边坡雷达影像与地形点云的误匹配校正方法，包括以下具体步骤：

步骤1，测量雷达轨道两端点坐标，结合雷达的设计参数计算天线行迹矢量

步骤2，使用几何映射方法完成雷达图像和三维地形点云数据粗匹配，得到粗匹配映射表trough；

步骤3，选取雷达图像中若干人工目标，测量得到人工目标的二维坐标集合{d1(r,θ),…,dn(r,θ)}，并定位人工目标在三维地形点云数据中的三维坐标集合{d1(x,y,z)t,…,dn(x,y,z)t}以作为人工目标的真实坐标；

步骤4，查找粗匹配映射表trough获得{d1(r,θ),…,dn(r,θ)}在粗匹配映射表中匹配的三维坐标集合{d1(x，y,z)s,…，dn(x，y，z)s}作为人工目标的偏差坐标；

步骤5，将人工目标的真实坐标和偏差坐标输入空间坐标变换方程，得到空间坐标变换方程的变换参数初值；

步骤6，进行显著地物和周边区域的相关性分析提取控制点对；

步骤7，利用步骤5中的变换参数初值、步骤6中的控制点对及最小二乘迭代估计空间坐标变换方程的变换参数的最优值，对粗匹配映射表trough进行变换完成误匹配校正。

进一步，步骤1中天线行迹矢量为：其中，为雷达轨道中轴向量，δh为俯仰调节旋钮中心距离轨道轴心高度，lantenna为天线中心至俯仰调节主轴距离，θsight为雷达视主轴与水平面夹角，和分别为以雷达合成孔径中心建立的坐标系os-xsyszs中ys和zs轴方向的空间单位向量。

进一步，坐标系os-xsyszs中xsosys为水平面，zs为铅垂线反方向，xs轴方向为合成孔径波束的指向在水平面的投影，ys轴与xs轴垂直并按右手坐标系定出。

进一步，步骤2具体为：

天线行迹矢量作为几何映射算法输入参数，利用距离-多普勒算法计算三维地形点云数据pterrain中各点与雷达的相对斜距和相对方位角；

通过插值搜索方法获得取雷达图像i(r,θ)中插值点的像素索引，得到粗匹配映射表trough，完成粗匹配。

进一步，pterrain中第i个顶点ai相对于孔径中心os的斜距：

式中表示ai相对于孔径中心os的斜距，||·||2代表向量的2-范数，n为pterrain中顶点总数；

以相对竖直中心平面左侧为方位负角度，pterrain中第i个顶点ai相对于孔径中心os的方位角为：

式中，表示ai相对于孔径中心os的方位角，为ai在天线行迹矢量的垂足，|·|代表向量的模。

进一步，步骤5中空间坐标变换方程为：其中为真实坐标，[txtytz]^t为平移因子，λ为缩放因子，rx(ωx)、ry(ωy)、rz(ωz)为各坐标轴的旋转矩阵，ωx、ωy、ωz为各坐标轴的偏转角，为偏差坐标；空间坐标变换方程的变换参数为平移因子、缩放因子和各坐标轴的偏转角。

进一步，

进一步，步骤6具体为：

提取雷达图像内带有监测场景内的显著地物的子图像i^sub以及三维地形点云数据中包含这些显著地物的子点云p^sub，查找trough内i^sub图像中心点匹配的三维坐标d^sub(x,y,z)t；

计算p^sub各点相对于雷达的斜距和方位角，斜距和方位角二维网格化生成网格，使用雷达入射角的经验散射模型计算网格点反射率，生成p^sub的反射率图i^σ，获取p^sub各点与i^σ中各点的映射表t^pcl；

使用相关性分析方法分析i^sub和i^σ，提取i^sub图像中心坐标与相关性峰值点坐标i^sub(r,θ)；

查找映射表t^pcl内i^sub(r,θ)对应的三维坐标d^sub(x,y,z)s，d^sub(x,y，z)t与d^sub(x，y，z)s形成一组非人工设置的控制点对；

重复上述步骤，获取包括若干控制点对的集合s。

进一步，雷达入射角的经验散射模型为σiδ为二维网格点的值，ηiδ代表雷达波入射角。

进一步，步骤7具体为：

利用步骤5中的变换参数初值、步骤6中的控制点，使用最小二乘迭代估计空间坐标变换方程的变换参数的最优值，将变换参数的最优值代入空间坐标变换方程形成变换方程f2；

将三维地形点云数据pterrain代入变换方程f2得到误匹配校正参考点集ptransform，两者之间的映射关系记为gtransform；

利用gtransform遍历pterrain与ptransform中各点的最近欧式距离，形成最近欧式距离vdist向量及最近欧式距离索引vidx向量；

剔除ptransform内最近欧式距离超过雷达图像像元空间分辨单元大小的点，保证误匹配纠正结果为像元级，将粗匹配结果prough内每点的形变值赋给pterrain内索引为vidx的点，形成pcorrected完成误匹配纠正。

本发明采用以上技术方案与现有技术相比，具有以下技术效果：

(1)利用辅助量测设备结合雷达设计参数计算天线行迹矢量，避免了常规方法仅仅测量轨道两端坐标造成的几何映射较大误差；

(2)通过场景内少量人工目标的匹配偏差确定空间变换方程参数初值，解决了空间变换方程初值选取困难，常规方法参数置零造成迭代难以收敛的问题；

(3)通过少量人工目标控制，利用场景内显著地物与周边区域的相关性分析提取大量控制点，算法简单，易于实现，解决了目标区域持续变形缺少足够的人工目标难题；

(4)通过大量控制点最小二乘迭代优化变换方程参数，并加以匹配距离约束，可达像元级匹配校正精度。

附图说明

图1是以雷达理想孔径中心为原点建立的用于匹配的测量观测坐标系示意图；

图2是边坡形变监测雷达系统与斜坡地形数据俯视视角几何关系模型示意图；

图3是本发明实施例提供的一种边坡形变监测雷达影像与地形点云基于少量人工控制点的误匹配校正方法的流程图。

具体实施方式

下面结合附图对本发明的技术方案做进一步的详细说明：

本发明一种边坡形变监测雷达影像与地形点云基于少量人工控制点的误匹配校正方法，如图3所示，包括以下步骤：

步骤1，利用辅助量测设备测量雷达轨道两端点坐标p1、p2(即雷达行进起始点和终止点限位坐标)，结合雷达的设计参数计算天线行迹矢量，使用几何映射方法完成雷达图像和地形数据粗匹配。所述雷达的设计参数包括合成孔径长度ls、步进间隔δl、方位向采样点数m、俯仰调节旋钮中心距离轨道轴心高度δh、天线中心至俯仰调节主轴距离lantenna。

天线行迹矢量：其中，为雷达轨道中轴向量，θsight为雷达视主轴与水平面夹角，分别为以雷达合成孔径中心建立的坐标系各轴空间单位向量。

步骤2，使用几何映射方法完成雷达图像和地形数据粗匹配，得到粗匹配映射表trough。

天线行迹矢量作为几何映射算法输入参数，利用常规距离-多普勒算法计算地形点云数据pterrain中各点与雷达的相对斜距r3d和相对方位角θ3d，通过常规最邻近搜索获得取雷达形变图像i(r,θ)中最近邻近像素索引(r2d,θ2d)，得到粗匹配映射表trough，完成粗匹配。得到粗匹配结果prough，prough包括坐标点和每点对应形变值。

步骤3，依据场景内的雷达和地形点云均可精确量测的少量人工目标的匹配偏差确定空间变换方程参数迭代初值，并进行显著地物和周边区域的相关性分析提取更多控制点。

通过识别雷达图像中为高亮度像元定位场景内具有少量人工的雷达波的强反射目标，获得二维坐标d1(r,θ),…,dn(r,θ)。

通过辅助量测设备定位人工目标在三维地形点云数据中的三维坐标d1(x，y，z)t，…，dn(x，y,z)t作为真实坐标。

查找粗匹配映射表trough，获得二维坐标d1(r,θ),…,dn(r,θ)在粗匹配映射表中匹配的d1(x,y,z)s,…,dn(x,y,z)s作为人工目标的偏差坐标。

将人工目标真实坐标和偏差坐标输入空间坐标变换方程f：获得七个空间坐标变换方程参数的初值λ⁰,方程中为真实坐标，λ为缩放因子，[txtytz]^t为平移因子，rx(ωx)、ry(ωy)、rz(ωz)为坐标系间各坐标轴的旋转矩阵，ωx、ωy、ωz为偏转角，为偏差坐标。

提取雷达图像内带有监测场景内的显著地物的子图像i^sub以及地形点云内包含这些显著地物的子点云p^sub，查找trough内i^sub图像中心点匹配的三维坐标d^sub(x，y,z)t。

计算p^sub各点相对于雷达的斜距和方位角，网格化生成二维网格，使用雷达入射角的经验散射模型计算网格点反射率生成各子点云区的反射率图i^σ，获取p^sub各点与i^σ图像网格点的映射表t^pcl。

使用常规相关性分析方法分析i^sub和i^σ，提取i^sub图像中心坐标与相关性峰值点坐标i^sub(r，θ)。

查找映射表t^pcl内i^sub(r，θ)对应的三维坐标d^sub(x，y,z)s，集合d^sub(x，y，z)t与d^sub(x，y，z)s形成一组识别的非人工设置的控制点对。

重复上述步骤，获取识别大量控制点对集合s。

步骤4，利用空间坐标变换方程参数初值、控制点对及最小二乘迭代估计空间变换方程最优参数，对几何映射方法的粗匹配结果进行变换完成误匹配校正。

利用大量控制点对s，变换参数初值使用最小二乘迭代寻优变换方程参数，获得优化的最终变换参数带回f形成变换方程f2。

将pterrain带入变换方程f2得到误匹配校正参考点集ptransform，映射关系记为gtransform。

利用映射gtransform遍历pterrain与ptransform各点的最近欧式距离，形成最近欧式距离vdist向量及最近欧式距离索引vidx向量。

查找vidx，剔除vdist超过雷达图像像元空间分辨单元大小的点，保证误匹配纠正结果为像元级，将prough每点形变量赋给pterrain内索引为vidx点，形成pcorrected完成误匹配纠正。

实施例

图1为本发明实施例提供的以雷达理想孔径中心为原点建立的用于匹配的测量观测坐标系示意图

利用全站仪、三维激光扫描、静态全球卫星导航定位系统(gps/gnss)等常用辅助量测设备测量雷达轨道两端坐标p1、p2，p1、p2分别为雷达行进起始点和终止点限位坐标，共观测np次，np≥3，第一次测量结果为第np次测量结果为平均测量结果

所述雷达的设计参数，包括合成孔径长度ls、步进间隔δl、方位向采样点数m、俯仰调节旋钮中心距离轨道轴心高度δh、天线中心至俯仰调节主轴距离lantenna。

天线行迹矢量计算：其中，为雷达轨道中轴向量。具体实施时通过对雷达方位向采样点坐标xm线性拟合，然后在向量投影计算，xm＝mδl-ls/2,m＝0,1,2,...,m-1，m为雷达方位向采样总点数。θsight为雷达视主轴与水平面夹角，分别为以雷达合成孔径中心建立的坐标系各轴空间单位向量，所述坐标系os-xsyszs为地基雷达架设施工使用的测量坐标系，如图1所示，孔径中心os，xsosys为水平面，zs为铅垂线反方向，xs方向为合成孔径波束的指向在水平面的投影，ys轴与xs垂直并按右手坐标系定出。

天线行迹矢量作为几何映射算法输入参数，利用常规距离-多普勒算法计算地形点云数据pterrain中各点与雷达的相对斜距r3d和相对方位角θ3d。

如图2，为本发明实施例提供的边坡形变监测雷达系统与斜坡地形数据俯视视角几何关系模型示意图。借助斜坡地形点云pterrain的各点量化几何映射方法。遍历pterrain中每个顶点相对于孔径中心os的斜距：

式中表示第i个点相对孔径中心的斜距，||·||2为取空间向量的2-范数，n为总点数。

遍历地形点云pterrain中每个顶点相对于孔径中心os的方位角，此处以相对竖直中心平面左侧为方位负角度，则方位角为：

式中，为第i个点相对孔径中心的方位角，为第i个点在天线行迹矢量的垂足，|·|为向量的模。

然后通过常规双线性插值方法即可获得三维空间点在二维雷达图像中的像素位置，将对应像素的形变测量结果赋值给三维空间点即实现了本实施例的几何映射粗匹配方法，得到粗匹配映射表trough，完成粗匹配。

通过识别雷达图像中为高亮度像元定位场景内具有少量人工目标如三面角反射器、金属靶球、金属箱、混凝土体等雷达波的强反射目标，获得二维坐标集合{d1(r，θ)，…，dn(r，θ)}，人工目标数量n≥3。

通过全站仪、三维激光扫描、静态gps/gnss等辅助量测设备定位人工目标在三维地形点云数据中的三维坐标集合{d1(x，y，z)t,…,dn(x,y,z)t},作为真实坐标。

进一步，查找粗匹配映射表trough，获得二维坐标{d1(r,θ),…,dn(r，θ)}在粗匹配映射表中匹配的{d1(x,y,z)s,…,dn(x,y,z)s}作为人工目标的偏差坐标。

将人工目标真实坐标和偏差坐标输入空间坐标变换方程f，本发明实施例以空间右手直角坐标系为例，首先将坐标绕xs轴旋转ωx，得到旋转矩阵rx(ωx)；再将坐标轴绕ys轴旋转ωy，得旋转矩阵ry(ωy)；最后将坐标轴绕zs轴旋转ωz，得到旋转矩阵rz(ωz)。三维坐标转换模型方程为：

三个人工目标的真实坐标集合和偏差坐标集合带入方程f，可求得七个变换参数的初值λ⁰，方程f中为真实坐标，λ为缩放因子，[txtytz]^t为平移因子，rx(ωx)、ry(ωy)、rz(ωz)为坐标系间各坐标轴的旋转矩阵，ωx、ωy、ωz为偏转角，为偏差坐标。

方程f参数优化需要更多控制点，本发明实施例提取雷达图像内带有监测场景内的显著地物如建筑、金属房、巨型岩石、挡墙等的子图像i^sub以及地形点云内包含这些显著地物的子点云p^sub，查找trough内i^sub图像中心点匹配的三维坐标d^sub(x,y,z)t。

计算p^sub各点相对于雷达的斜距和方位角，网格化生成二维网格，使用雷达入射角的经验散射模型，常用的经验散射模型有余弦模型、余弦平方模型，本实施例选取经验模型：σiδ为二维网格点的值，ηiδ代表雷达波入射角。本实施例可通过对p^sub各点与邻近的k个点拟合小面元平面d，计算os与p^sub各点所成向量与小面元d的法向量夹角求取雷达波入射角，k≥6。

计算网格点反射率生成各子点云区的反射率图i^σ，获取p^sub各点与i^σ图像网格点的映射表t^pcl。

使用归一化互相关性分析方法分析i^sub和i^σ，：

式中，c为互相关系数，在0.1～1之间；i^σ(x,y)为图像i^σ网格点(x,y)的灰度值；i^sub(x,y)为i^sub图像像点(x，y)的灰度值；μt为i^σ灰度均值；μi为图像i^sub的灰度均值。提取i^sub图像中心坐标与相关性峰值点坐标i^sub(r,θ)。

查找映射表t^pcl内i^sub(r,θ)对应的三维坐标d^sub(x,y,z)s，集合d^sub(x,y，z)t与d^sub(x,y,z)s形成一组识别的非人工设置的控制点对。重复上述步骤，获取识别大量控制点对集合s。

利用大量控制点对s及变换参数初值λ⁰,使用最小二乘迭代寻优变换方程参数，本实施例具体如下：

对公式(3)在λ⁰，参数初值处进行泰勒展开，保留一阶项得到：

a1＝-sinωzcosωydωz-cosωzsinωydωy，b1＝-cosωzcosωydωz+sinωzsinωydωy，c1＝cosωydωy。

线性化变换公式(5)得到误差方程：

式中，为变换参数的改正数，系数矩阵常数项矩阵

最小二乘迭代令公式(6)误差最小。通过迭代计算控制舍入误差，舍入误差选取单位权中误差自由度为误差方程数量与实际求解方程所需变量数的差，各控制点视为独立观测量权阵p取单位矩阵，求出七个变换参数的最佳估值。实施例迭代计算过程：

将参数初值带入xt,r′，组成式6误差方程。n个控制点将组成3n个误差方程。

利用最小二乘方法求取变换参数改正数x^(∈+1)(∈代表迭代次数)④检核参数的改正数是否小于给定的限差，本实施例取平移量mm,偏转角10^-3″和缩放系数10^-7为限差，若未达到限差则将x^(∈)＝x^(∈-1)+x^(∈)作为新的初值，重复迭代；符合限差，则计算结束，将x^∈作为参数最佳估值，获得优化的最终变换参数带回f形成变换方程f2。

将pterrain带入变换方程f2得到误匹配校正参考点集ptransform，映射关系记为gtransform。

利用映射gtransform遍历pterrain与ptransform各点最近欧式距离。形成距离vdist向量及最近邻索引vidx向量。本实施例引入常用的点云结果化处理方法k-d树(kdimensionaltree)和k近邻(k-nearestneighbor,knn)聚类方法。在处理过程中，利用knn查找k个最接近的点，令k＝1，即求得ptransform与pterrain中最邻近点。

查找vidx，剔除ptransform内最近欧式距离超过雷达图像像元空间分辨单元大小的点，本发明实施例空间分辨单元大小指雷达二维图像像元匹配至三维空间距离向和方位向覆盖范围的平均长度，小于该长度的匹配点将保留，保证误匹配纠正结果为像元级，将prough每点形变量赋给pterrain内索引为vidx点，形成pcorrected完成误匹配纠正。

以上所述，仅为本发明中的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉该技术的人在本发明所揭露的技术范围内，可理解想到的变换或替换，都应涵盖在本发明的包含范围之内，因此，本发明的保护范围应该以权利要求书的保护范围为准。