CN113516606B - 倾斜摄影数据辅助边坡雷达影像与光学照片匹配融合方法 - Google Patents

Abstract

本发明公开了一种倾斜摄影数据辅助边坡雷达影像与光学照片匹配融合方法，首先，利用地基雷达设站信息及倾斜影像的空间关系完成几何映射；然后，在倾斜摄影数据中依据照片拍摄点坐标及朝向透视投影变换得到瞬时视角图像，提取瞬时视角图像和拍摄照片公共点后使用最小二乘优化的仿射变换匹配，获得雷达测量图像和拍摄照片像元间映射关系；最后，引入图像融合算法，将实时形变与拍摄照片颜色空间融合得到拍摄照片‑形变融合图。本发明利用倾斜影像的三维信息辅助完成地基雷达图像和光学拍照匹配，避免了常规的方法直接利用二维图像变换方法忽略成像原理差异造成的巨大匹配误差问题。

Description

倾斜摄影数据辅助边坡雷达影像与光学照片匹配融合方法

技术领域

本发明涉及一种倾斜摄影数据辅助边坡雷达影像与光学照片匹配融合方法，属于露天矿边坡形变监测技术领域。

背景技术

滑坡是仅次于地震、发生最频繁、造成损失最严重的一种地质灾害。滑坡灾后现场应急抢险搜救人员及工程车辆众多，此时若发生次生灾害，损失将无法估量。对滑坡的残余岩体持续监测，分析变形特征，专家研判定位危险变形位置并及时预警是目前滑坡灾后现场应急监测广泛认可的技术路线。地基合成孔径干涉雷达(GB-InSAR)是近些年的热点技术，已被证明是区域性形变监测的有力工具，具有较短的时间间隔(重访周期可达分钟量级)。已有的滑坡监测案例表明实测精度可达到亚毫米级，GB-InSAR在滑坡应急监测预警领域有巨大前景。雷达在极坐标系成像，需映射至三维空间才便于人眼感知定位危险区域，不适合安全监测人员直接使用，因此，测量结果常被要求与监测人员拍摄照片相匹配。

现有匹配方式匹配精度低，将地基雷达图像与拍摄的影像通过二维图像变换进行匹配，二者成像原理差距大无法精确对应，难以作为判读参考。边坡现场常常具有三维激光扫描(LIDAR)和无人机倾斜摄影三维数据。三维数据一方面可以通过距离-多普勒分析通过几何映射与雷达图像匹配，另一方面倾斜摄影的瞬时视角图像又可以和拍摄照片进行匹配。因此可以通过倾斜摄影辅助完成地基边坡形变监测雷达图像和拍摄照片匹配融合。

发明内容

为克服现有技术的不足，本发明提供一种倾斜摄影作为中间数据，分别与边坡形变监测雷达影像及拍摄照片匹配，实现边坡雷达形变影像和光学照片匹配融合的方法。

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

倾斜摄影数据辅助边坡雷达影像与光学照片匹配融合方法，包括以下具体步骤：

步骤1，连续测量地基雷达一步一停移动过程中的每个停止位的天线中心坐标，得到点集P_rail，对点集P_rail进行直线拟合得到天线行迹矢量

步骤2，计算三维倾斜摄影数据点集P_map中各点与雷达的相对斜距和相对方位角，形成二维点集P_map3D；

步骤3，计算P_map3D中各点与雷达图像中各像素点的欧式距离，最小欧式距离对应的像素点构成获得最邻近像素点集P_I，得到P_map3D与P_I间一一映射表T_I；

步骤4，从光学照片I_shot内读取光学照片拍摄点的地理坐标P_shot，将P_map的观察点设置在P_shot处，使用投影变换法调整三维倾斜摄影数据的观察视角得到的瞬时图像I_temp，保存P_map与I_temp各点间的映射表T_temp；

步骤5，目视解译I_temp与I_shot中显著地物目标的二维坐标，得到粗匹配公共点集P_co，其中，P_co由两个子点集构成，一个为I_temp内显著地物目标的坐标子点集，另一个为I_shot内与I_temp相同的显著地物目标的坐标子点集；

步骤6，将P_co输入图像仿射变换方程f，得到图像仿射变换方程的变换参数初值，其中变换参数包括旋转因子、平移因子以及为缩放因子；

步骤7，将步骤6中的变换参数初值代回f形成变换方程f1，将I_shot中所有像素点代入f1，得到粗匹配图像I_rough，形成I_shot与I_rough各像素间的粗匹配映射表T_rough，完成粗匹配；

步骤8，使用特征提取方法提取I_temp和I_rough内的图像特征，得到若干I_temp和I_rough中的相同特征点，形成精匹配公共点对集P_fine，P_fine中元素为I_temp和I_rough中相同特征点的像元坐标；

步骤9，利用步骤8中的P_fine、步骤6中的变换参数初值，使用最小二乘迭代估计图像仿射变换方程f的最优变换参数，将最优变换参数代回f形成变换方程f2；

步骤10，将I_rough代入变换方程f2得到I_fine，形成I_rough与I_fine各像素间的精匹配映射表T_fine，通过重采样与插值方法即可得到I_shot与雷达图像P_I间的映射关系T_final；

步骤11，查找T_final映射表，通过图像融合方法将雷达图像中各像素点的形变值与I_shot中RGB颜色通道进行融合，获得融合图像I_fusion，完成匹配融合过程。

进一步，步骤1中的方向由雷达运行起始点指向终止点。

进一步，步骤2中利用距离-多普勒算法计算P_map中各点与雷达的相对斜距和相对方位角。

进一步，P_map中第i个顶点A_i与雷达的相对斜距：

式中，N表示P_map中的顶点数，(x_i,y_i,z_i)为A_i的三维坐标，(x_s,y_s,z_s)为雷达合成孔径中心O_s的三维坐标；

P_map中第i个顶点A_i与雷达的相对方位角为：

式中，为A_i在天线行迹矢量/>上的垂足，|·|代表向量的模，P₂为雷达行进终止点。

进一步，步骤6中图像仿射变换方程f为：

其中为目标坐标，/>为待匹配坐标，R为旋转因子，T_x、T_y为平移因子，ρ为缩放因子。

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

(1)利用倾斜影像的三维信息辅助完成地基雷达图像和光学拍照匹配，避免了常规的方法直接利用二维图像变换方法忽略成像原理差异造成的巨大匹配误差问题；

(2)通过识别场景内显著地物先进行倾斜影像瞬时图像与光学拍照图粗匹配后利用特征识别大量公共点优化仿射变换方程参数，解决了常规变换方法初值选择困难问题；

(3)通过特征识别大量公共点最小二乘迭代优化变换方程参数，方法易于实现，运算速度较快，可实现像元级匹配融合精度；

(4)边坡形变监测雷达的形变测量图像实时地与地质灾害勘察人员拍摄照片匹配，易用性强，可为滑坡风险评估，灾害预警提供数据支撑。

附图说明

图1是倾斜摄影数据辅助的边坡形变监测雷达影像与光学照片匹配融合示意图；

图2是本发明实施例提供的一种倾斜摄影数据辅助的边坡形变监测雷达影像与光学照片匹配融合方法流程图。

具体实施方式

下面根据说明书附图对本发明的发明内容做进一步详细表述。显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部实施例。基于本发明的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明的保护范围。

本发明倾斜摄影数据辅助边坡雷达影像与光学照片匹配融合方法，首先，利用地基雷达设站坐标、俯仰角信息及倾斜影像的空间关系完成几何映射；然后，在倾斜摄影数据中依据照片拍摄点坐标及朝向透视投影变换得到瞬时视角图像，提取瞬时视角图像和拍摄光学照片公共点后使用最小二乘优化的仿射变换匹配，获得雷达图像和拍摄光学照片像元间映射关系；最后，引入图像融合算法，将边坡形变监测雷达图像与拍摄光学照片颜色空间融合得到拍摄照片-形变融合图。

实施例

实施例提供的一种倾斜摄影数据辅助的边坡形变监测雷达影像与光学照片匹配融合方法流程图，如图2所示。

利用全站仪、三维激光扫描、静态全球卫星导航定位***(GPS/GNSS)等常用的辅助量测设备连续测量地基雷达一步一停移动过程中的每个停止位的天线中心坐标。雷达初次运行测量点集为第n次运行测量点集为/>集合/>得到点集P_rail。

利用常规的直线拟合方法对点集P_rail进行拟合得到雷达天线行迹矢量矢量/>方向由雷达运行起始点指向终止点。

天线行迹矢量三维倾斜摄影数据点集P_map作为几何映射算法输入参数，利用常规距离-多普勒算法计算点集P_map中各点与雷达的相对斜距R_3D和相对方位角θ_3D。

令雷达合成孔径中心为O_s，遍历P_map中每个顶点(x_i,y_i,z_i)相对于孔径中心O_s(x_s,y_s,z_s)的欧式距离：

式中表示第i个点相对孔径中心O_s的斜距，N表示P_map模型总点数。

遍历P_map中每个模型顶点相对于孔径中心O_s的方位角，以相对竖直中心平面左侧为雷达方位负角度，则方位角为：

式中，为第i个点相对孔径中心的方位角，/>为第i个点在行迹轴线/>的垂足。|·|代表向量的模。形成二维点集P_map3D:{(R_3D,θ_3D)₁,(R_3D,θ_3D)₂,...,(R_3D,θ_3D)_N}，计算P_map3D各点与雷达图像I(r,θ)各网格点的欧式距离取最小欧式距离点获得最邻近像素点集P_I:{(R_2D,θ_2D)₁,(R_2D,θ_2D)₂,...,(R_2D,θ_2D)_N}，得到P_map3D与P_I间一一映射表T_I。

本实施例光学照片I_shot拍摄地点的地理坐标通过GPS/GNSS获取，读取拍摄照片点的地理坐标P_shot，倾斜摄影数据将观察点设置在P_shot处，再确定视线方向即可确定倾斜摄影数据的瞬时二维图像I_temp。实施例/>为观察点P_shot与场景中心P_R0所成向量。投影变换调整倾斜摄影数据观察视角至瞬时图像I_temp与拍摄照片I_shot内容初步一致。本实施例中投影变换，倾斜摄影初始为世界坐标系O-XYZ，先变换至以观察点为中心的目坐标系P_shot-X_eY_eZ_e，负Z_e轴指向观察方向/>世界坐标系和目坐标系均取右手三维空间直角坐标系。然后，将P_map中各点投影至屏幕空间，所述投影指目标物描述从世界坐标系到目坐标系的变换等价与运用平移、旋转运算将目坐标系叠加到世界坐标系上的变换，再使用透视投影。将平行与X_eO_eY_e平面且离视点的距离等于f的平面作为瞬时图像I_temp投影面，那么，目坐标系中的一点在瞬时图像I_temp投影面上的网格点坐标(X_m,Y_m)可由下式进行计算：

式中，(X_s,Y_s,Z_s)为P_shot在世界坐标系O-XYZ中的坐标，瞬时图像I_temp投影面相对于平面XOY方位角θ和俯仰角α，观察点P_shot与瞬时图像I_temp投影面距离f，世界坐标系中任意目标点(X_M,Y_M,Z_M)，保存P_map与I_temp各网格点间的映射表T_temp。

目视解译瞬时图像I_temp与拍摄照片I_shot图像中显著地物目标的二维坐标得到粗匹配公共点集P_co，点集P_co包括为I_temp内显著地物目标的坐标子点集，还包括/>为I_shot内与I_temp相同的地物目标的坐标子点集，共np对公共点。

将粗匹配公共点集P_co输入图像仿射变换方程得到参数初值ρ⁰、R⁰、/>方程f中/>为目标坐标，/>为待匹配坐标，R为旋转因子，T_x、T_y为平移因子，ρ为缩放因子。

将参数初值带回方程f形成变换方程f1，将I_shot的所有图像网格点带入f2，得到粗匹配图像网格I_rough，得到粗匹配映射表T_rough，完成粗匹配。

在粗匹配结果控制下，使用特征提取方法提取I_temp和I_rough内的点、线、面特征，得到大量相同特征点作为精匹配公共点集P_fine。

利用精匹配公共点集P_fine，变换参数初值ρ⁰、R⁰、使用最小二乘迭代寻优仿射变换方程f1参数，本实施例中最小二乘仿射变换使用迭代法优化，以平面右手直角坐标系为例，本实施例仿射变换后再进行双线性插值变换，变换模型为：

建立误差方程：

其中，n为公共点P_fine个数，则模型参数的一种最小二乘解为/>P为等权重矩阵。最小二乘迭代,令公式(4)误差最小。通过迭代计算控制舍入误差，舍入误差选取单位权中误差/>自由度/>为误差方程数量与实际求解方程所需变量数的差，各控制点视为独立观测量权阵P取单位矩阵，迭代求出变换参数的最佳估值/>带回f形成变换方程f2。

将I_rough带入变换方程f2得到I_fine，映射关系记为T_fine，此时已获得边坡形变监测雷达图像P_I与P_map3D间一一映射表T_I，P_map3D和P_map为一一对应关系，P_map与I_temp各点间的映射表T_temp，I_shot与I_rough粗匹配映射表T_rough，I_rough与I_fine精匹配映射表T_fine，通过常规重采样与插值方法即可得到I_shot光学照片与雷达图像P_I间的映射关系T_final。

查找T_final映射表，将边坡形变雷达的形变图像网格点形变值与I_shot光学照片红绿蓝(RGB)颜色通道通过图像融合方法进行融合。本实施例采用常用融合方法使雷达测量形变量替换光学照片红色(R)通道

式中，为经过变换后的图像，有红、绿、蓝三通道，红色通道使用了雷达图像P_I形变值IR，绿色通道和蓝色通道使用I_shot光学照片的绿色通道I_VG和蓝色通道I_VB，获得融合图像I_fusion，完成匹配融合过程。倾斜摄影数据辅助的边坡形变监测雷达影像与光学照片匹配融合示意图，如图1所示。

Claims (5)

1.倾斜摄影数据辅助边坡雷达影像与光学照片匹配融合方法，其特征在于，包括以下具体步骤：

步骤1，连续测量地基雷达一步一停移动过程中的每个停止位的天线中心坐标，得到点集P_rail，对点集P_rail进行直线拟合得到天线行迹矢量

步骤2，计算三维倾斜摄影数据点集P_map中各点与雷达的相对斜距和相对方位角，形成二维点集P_map3D；

步骤3，计算P_map3D中各点与雷达图像中各像素点的欧式距离，最小欧式距离对应的像素点构成获得最邻近像素点集P_I，得到P_map3D与P_I间一一映射表T_I；

步骤4，从光学照片I_shot内读取光学照片拍摄点的地理坐标P_shot，将P_map的观察点设置在P_shot处，使用投影变换法调整三维倾斜摄影数据的观察视角得到的瞬时图像I_temp，保存P_map与I_temp各点间的映射表T_temp；

步骤5，目视解译I_temp与I_shot中显著地物目标的二维坐标，得到粗匹配公共点集P_co，其中，P_co由两个子点集构成，一个为I_temp内显著地物目标的坐标子点集，另一个为I_shot内与I_temp相同的显著地物目标的坐标子点集；

步骤6，将P_co输入图像仿射变换方程f，得到图像仿射变换方程的变换参数初值，其中变换参数包括旋转因子、平移因子以及缩放因子；

步骤7，将步骤6中的变换参数初值代回f形成变换方程f1，将I_shot中所有像素点代入f1，得到粗匹配图像I_rough，形成I_shot与I_rough各像素间的粗匹配映射表T_rough，完成粗匹配；

步骤8，使用特征提取方法提取I_temp和I_rough内的图像特征，得到若干I_temp和I_rough中的相同特征点，形成精匹配公共点对集P_fine，P_fine中元素为I_temp和I_rough中相同特征点的像元坐标；

步骤9，利用步骤8中的P_fine、步骤6中的变换参数初值，使用最小二乘迭代估计图像仿射变换方程f的最优变换参数，将最优变换参数代回f形成变换方程f2；

步骤10，将I_rough代入变换方程f2得到I_fine，形成I_rough与I_fine各像素间的精匹配映射表T_fine，通过重采样与插值方法即可得到I_shot与雷达图像P_I间的映射关系T_final；

步骤11，查找T_final映射表，通过图像融合方法将雷达图像中各像素点的形变值与I_shot中RGB颜色通道进行融合，获得融合图像I_fusion，完成匹配融合过程。

2.如权利要求1所述的倾斜摄影数据辅助边坡雷达影像与光学照片匹配融合方法，其特征在于，步骤1中的方向由雷达运行起始点指向终止点。

3.如权利要求1所述的倾斜摄影数据辅助边坡雷达影像与光学照片匹配融合方法，其特征在于，步骤2中利用距离-多普勒算法计算P_map中各点与雷达的相对斜距和相对方位角。

4.如权利要求3所述的倾斜摄影数据辅助边坡雷达影像与光学照片匹配融合方法，其特征在于，P_map中第i个顶点A_i与雷达的相对斜距：

式中，N表示P_map中的顶点数，(x_i,y_i,z_i)为A_i的三维坐标，(x_s,y_s,z_s)为雷达合成孔径中心O_s的三维坐标；

P_map中第i个顶点A_i与雷达的相对方位角为：

式中，为A_i在天线行迹矢量/>上的垂足，|·|代表向量的模，P₂为雷达行进终止点。

5.如权利要求1所述的倾斜摄影数据辅助边坡雷达影像与光学照片匹配融合方法，其特征在于，步骤6中图像仿射变换方程f为：

其中为目标坐标，/>为待匹配坐标，R为旋转因子，T_x、T_y为平移因子，ρ为缩放因子。