一种基于雷达影像与电子海图信息的岸线配准方法

摘要

本发明公开了一种基于雷达影像与电子海图信息的岸线配准方法，包括以下步骤：对雷达影像进行处理，获取其中的河/海岸线，并转换为全局坐标系下的像素位置坐标；根据S57海图格式标准，从电子海图中提取代表河/海岸线的线条信息，并转换为大地坐标系下的经纬度坐标；通过灰色关联处理手段，针对雷达影像和电子海图所提取出来的岸线坐标进行关联匹配，比较不同位置及尺度下的数据关联度，在这两种不同来源的数据之间实现粗配准；通过ICP配准处理手段，针对雷达影像和电子海图所提取出来的岸线坐标进行精配准。本发明提供的岸线识别方法识别精确度高，为其他目标识别划定范围做铺垫，有效降低雷达目标识别的错误率。

说明书

技术领域

本发明涉及地理信息技术领域，尤其涉及一种基于雷达影像与电子海图信 息的岸线配准方法。

背景技术

在船舶现代化监管中，海事雷达是一种非常重要的监管手段。由于其工作 环境复杂，岸线周围错综复杂的建筑物、行驶车辆、风浪、水面等往往会对海 事雷达产生较大的干扰，出现杂波，有必要针对雷达图像进行处理，生成岸线 的基本轮廓，进而为岸线识别，以及水上目标识别做好铺垫，以便降低雷达目 标识别的错误率。此外，对雷达影像和电子海图的叠加显示，能够将电子海图 的地理信息优势和雷达装备的实时探测能力结合在一起，从而更加准确地获悉 周围水域的障碍物信息，大幅提高船舶的避碰能力。因此，通过图像处理技术， 从雷达影像中提取出岸线数据，并与电子海图进行配准，实现雷达影像和电子 海图的叠加显示，在保障水上交通安全方面具有重要意义。

发明内容

本发明要解决的技术问题在于针对现有技术中的缺陷，提供一种基于雷达 影像与电子海图信息的岸线配准方法。

本发明解决其技术问题所采用的技术方案是：一种基于雷达影像与电子海 图信息的岸线配准方法，包括以下步骤：

1)对雷达影像进行处理，获取其中的河/海岸线，并转换为全局坐标系下的 像素位置坐标；

2)根据S57海图格式标准，从电子海图中提取代表河/海岸线的线条信息， 并转换为大地坐标系下的经纬度坐标；

3)通过灰色关联处理手段，针对雷达影像和电子海图所提取出来的岸线坐 标进行关联匹配，比较不同位置及尺度下的数据关联度，在这两种不同来源的 数据之间实现粗配准；

4)通过ICP配准处理手段，针对雷达影像和电子海图所提取出来的岸线坐 标进行精配准。

按上述方案，所述步骤1)中雷达影像的图像处理使用包括连通域判别、边 缘提取、Gabor小波处理图像处理技术。

按上述方案，所述步骤(1)中基于雷达影像的河/海岸线信息提取方法，包 括以下步骤：

(1.1)对雷达影像进行灰度变换，将彩色图像转化为8位灰度图像；

(1.2)对雷达影像进行连续帧间差分处理，提取雷达检测中的背景影像； 定义第t帧雷达检测影像为：其中i、j代表第i行，第j列的像素点；则连 续三帧雷达影像可以分别表示为：同时定义第t帧的背景 图像为：

(1.3)将这三帧连续视频图像分为两组，第t-1帧和第t帧为一组，第t 帧及第t+1帧为另一组，分别对这两组视频的每一个象素点进行帧间差减，并 将其绝对值保存在中，如公式(1)所示：

$D_{i,j}^{-1}=|I_{i,j}^{t-1}-I_{i,j}^t|,D_{i,j}^{+1}=|I_{i,j}^1-I_{i,j}^{t+1}|---\left(1\right)$

(1.4)将中各像素的值分别与一个预先给定的阈值T₀进行比 较，如果某一点i，j上对应的的值均大于T₀，则可以判断这一点 在当前连续三帧图像中都是运动的；反之，该点即为背景图像中的一部分， 如公式(2)所示，对前30秒的雷达检测影像反复进行背景提取处理，即可构 建出完整的初始背景影像；

$B_{i,j}^t=\left(\begin{array}{cc}{B}_{i,j}^{t-1},& \mathrm{if}{D}_{i,j}^{-1}>T_0\mathrm{and}{D}_{i,j}^{+1}>T_0\\ I_{i,j}^t,& \mathrm{otherwise}\end{array}\right)---\left(2\right)$

(1.5)统计背景影像的灰度分布直方图，设定阈值T₁，针对目标信息的灰 度特性进行二值化；如公式(3)所示：

$B_{i,j}^t=\left(\begin{array}{cc}255& \mathrm{if}{B}_{i,j}^t>T_1\\ 0,& \mathrm{otherwise}\end{array}\right)---\left(3\right)$

(1.6)在得到二值化的背景图象后，利用区域生长法将灰度值为255的像 素点按邻域连通性原则合成一个区域，对连通区域像素值大小的统计和排序， 如果区域面积过小或者长宽比不符合岸线的一般特征，则认为是干扰信号而予 以删除，最终选取出包含了岸线有效信息的连通区域；其中，连通区域的选取 采用四领域法，即目标像素p及其周围的(x+1,y)，(x,y-1)，(x-1,y)， (x,y+1)呈现联通状态则认为是同一连通域；

(1.7)采用拉普拉斯算子对包含岸线有效信息的连通域进行边缘提取，获 取可能的岸线边缘；拉普拉斯算子的定义如公式(4)、(5)所示：

${▿}^{2}f=4z_5-(z_2+z_4+z_6+z_8)---\left(4\right)$

${▿}^{2}f=8z_5-(z_1+z_2+z_3+z_4+z_5+z_6+z_7+z_8)---\left(4\right)$

(1.8)对背景边缘图像采用Hough变换获取岸线的方向信息，Hough变换 的极坐标定义为：x cosθ+y sinθ＝ρ，通过对所有边缘像素点的极坐标(ρ,θ)进 行统计，选取最大的θ作为岸线的方向角信息，并以其为依据确定Gabor变换的 角度，代入Gabor小波中实现对岸线点的提取；二维Gabor小波变换以及傅里 叶变换的定义如公式(6)、(7)、(8)所示：

$g(x,y)=\frac{1}{2\pi {\sigma }_x{\sigma }_y}\exp (-\frac{1}{2}(\frac{x^2}{{\sigma }_x^2}+\frac{y^2}{{\sigma }_y^2}))---\left(6\right)$

h(x,y)＝g(x,y)exp(2πjWx) (7)

$H(x,y)=\exp (-\frac{1}{2}(\frac{{(u-W)}^2}{{\sigma }_u^2}+\frac{v^2}{{\sigma }_v^2}))---\left(8\right).$

按上述方案，所述步骤(3)中针对来自雷达影像和电子海图的岸线信息进 行灰色关联分析，确定几何形态最为接近的岸线位置，实现两种岸线数据间的 粗配准，具体包括以下步骤：

(3.1)针对雷达影像和电子海图的特征，确定灰色关联分析中的参考数据 列和比较数据列：由于雷达影像为位图形式，而电子海图为矢量图形式，因此 将雷达影像中提取的岸线坐标设定为参考数据列，定义为：R＝{r(k)|k＝ 1,2,…,n，将电子海图中的岸线信息设定为比较数据列，分别定义为： E_i＝{e_i(k)|k＝1,2,…,n}。之后比较R和若干个E_i之间的几何形态相似程度，判 断其关联程度；

(3.2)雷达影像中的岸线信息为位图坐标，电子海图中的岸线信息为经纬 度信息，数据量纲不同，需根据公式(9)进行数据的无量纲化处理如下：

$e_i\left(k\right)=\frac{e_i\left(k\right)}{e_i\left(1\right)},r\left(k\right)=\frac{r\left(k\right)}{r\left(1\right)},i=\mathrm{1,2},...,m;k\mathrm{1,2},...,n---\left(9\right)$

(3.3)根据公式(10)和(11)，计算比较数据列R与参考数据列E_i中对应 数据的关联系数如下：

${\xi }_i\left(k\right)=\frac{{\min }_i\underset{k}{\min }{\Delta }_i\left(k\right)\rho {\max }_i{\max }_k{\Delta }_i\left(k\right)}{{\Delta }_i\left(k\right)+\rho {\max }_i{\max }_k{\Delta }_i\left(k\right)}---\left(10\right)$

Δ_i(k)＝|r(k)-e_i(k)|  (11)

(3.4)取所有关联系数的平均值为数据列之间的关联度，并以此表示比较 数据列R与参考数据列E_i之间的关联程度，如公式(12)如示：

$r_i=\frac{1}{n}{\Sigma }_{k=1}^n{\xi }_i\left(k\right),k=\mathrm{1,2},...,n---\left(12\right)$

(3.5)将关联度r_i按大小排序，取关联度最大的比较数据列R与参考数据列 E_i为同一条岸线的粗配准结果。

按上述方案，所述步骤(4)中针对来自雷达影像和电子海图的岸线信息采 用了改进的ICP点云配准方法，实现了两种岸线数据间的精配准，具体包括以 下步骤：

(4.1)将粗配准后的雷达影像岸线数据存储在点集P中，定义为： {P|P_i∈R³,i＝1,2,…,N}，p_i为岸线上某点的坐标，定义为：(x_pi，y_pi)；同时， 将对应的电子海图岸线数据存储在点集Q中，定义为：{Q|q_i∈R³,i＝1,2,…,N}， q_i为岸线上某点的坐标，定义为：(x_qi,y_qi)；作为同一段岸线的非同源数据，P、 Q两个点集之间应符合刚体变换，假定该刚体变换是由一个2*2的旋转矩阵R 和一个2*2的平移矩阵T实现的，如公式(13)所示：

$\left(\begin{array}{c}{x}_{\mathrm{pi}}\\ y_{\mathrm{pi}}\end{array}\right)=R\left(\begin{array}{c}{\lambda }_{\mathrm{qi}}\\ y_{\mathrm{pi}}\end{array}\right)+T---\left(13\right)$

(4.2)根据最优匹配原则设计数据配准的目标函数，如公式(14)所示， 点集R、T应使该目标函数的值最小：

$f(R,T)={\Sigma }_{i=1}^n{\left|\right|P_i^k-(R{Q}_i+T)\left|\right|}^2=\min ---\left(14\right)$

(4.3)计算点集P中的每一个点在Q点集中的对应最近点；

(4.4)计算上述对应点对平均距离最小的刚体变换，求得旋转矩阵R和平 移矩阵T；

(4.5)对Q使用上一步的平移矩阵T和旋转矩阵R，得到新的变换点集， 并替换获得新的点集Q；

(4.6)如果新的点集Q与参考点集P之间满足公式(14)中目标函数的要 求，即平均距离小于给定的阀值，则停止迭代，完成精配准；否则，以新的点 集Q与参考点集P继续迭代下去，直到达到目标函数的要求。

本发明产生的有益效果是：本发明提供的岸线识别方法识别精确度高，为 其他目标识别划定范围做铺垫，有效降低雷达目标识别的错误率。

附图说明

下面将结合附图及实施例对本发明作进一步说明，附图中：

图1是本发明实施例的方法流程图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例， 对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解 释本发明，并不用于限定本发明。

如图1所示，一种基于雷达影像与电子海图信息的岸线配准方法，包括以 下步骤：

1)使用包括连通域判别、边缘提取、Gabor小波处理图像处理技术对雷达 影像进行处理，获取其中的河/海岸线，并转换为全局坐标系下的像素位置坐标；

(1.1)对雷达影像进行灰度变换，将彩色图像转化为8位灰度图像；

(1.2)对雷达影像进行连续帧间差分处理，提取雷达检测中的背景影像； 定义第t帧雷达检测影像为：其中i、j代表第i行，第j列的像素点；则连 续三帧雷达影像可以分别表示为：同时定义第t帧的背景 图像为：

(1.3)将这三帧连续视频图像分为两组，第t-1帧和第t帧为一组，第t 帧及第t+1帧为另一组，分别对这两组视频的每一个象素点进行帧间差减，并 将其绝对值保存在中，如公式(1)所示：

$D_{i,j}^{-1}=|I_{i,j}^{t-1}-I_{i,j}^t|,D_{i,j}^{+1}=|I_{i,j}^1-I_{i,j}^{t+1}|---\left(1\right)$

(1.4)将中各像素的值分别与一个预先给定的阈值T₀进行比 较，如果某一点i，j上对应的的值均大于T₀，则可以判断这一点 在当前连续三帧图像中都是运动的；反之，该点即为背景图像中的一部分， 如公式(2)所示，对前30秒的雷达检测影像反复进行背景提取处理，即可构 建出完整的初始背景影像；

$B_{i,j}^t=\left(\begin{array}{cc}{B}_{i,j}^{t-1},& \mathrm{if}{D}_{i,j}^{-1}>T_0\mathrm{and}{D}_{i,j}^{+1}>T_0\\ I_{i,j}^t,& \mathrm{otherwise}\end{array}\right)---\left(2\right)$

(1.5)统计背景影像的灰度分布直方图，设定阈值T₁，针对目标信息的灰 度特性进行二值化；如公式(3)所示：

$B_{i,j}^t=\left(\begin{array}{cc}255& \mathrm{if}{B}_{i,j}^t>T_1\\ 0,& \mathrm{otherwise}\end{array}\right)---\left(3\right)$

(1.6)在得到二值化的背景图象后，利用区域生长法将灰度值为255的像 素点按邻域连通性原则合成一个区域，对连通区域像素值大小的统计和排序， 如果区域面积过小或者长宽比不符合岸线的一般特征，则认为是干扰信号而予 以删除，最终选取出包含了岸线有效信息的连通区域；其中，连通区域的选取 采用四领域法，即目标像素p及其周围的(x+1,y)，(x,y-1)，(x-1,y)， (x,y+1)呈现联通状态则认为是同一连通域；

(1.7)采用拉普拉斯算子对包含岸线有效信息的连通域进行边缘提取，获 取可能的岸线边缘；拉普拉斯算子的定义如公式(4)、(5)所示：

${▿}^{2}f=4z_5-(z_2+z_4+z_6+z_8)---\left(4\right)$

${▿}^{2}f=8z_5-(z_1+z_2+z_3+z_4+z_5+z_6+z_7+z_8)---\left(4\right)$

(1.8)对背景边缘图像采用Hough变换获取岸线的方向信息，Hough变换 的极坐标定义为：x cosθ+y sinθ＝ρ，通过对所有边缘像素点的极坐标(ρ,θ)进 行统计，选取最大的θ作为岸线的方向角信息，并以其为依据确定Gabor变换的 角度，代入Gabor小波中实现对岸线点的提取；二维Gabor小波变换以及傅里 叶变换的定义如公式(6)、(7)、(8)所示：

$g(x,y)=\frac{1}{2\pi {\sigma }_x{\sigma }_y}\exp (-\frac{1}{2}(\frac{x^2}{{\sigma }_x^2}+\frac{y^2}{{\sigma }_y^2}))---\left(6\right)$

h(x,y)＝g(x,y)exp(2πjWx)  (7)

$H(x,y)=\exp (-\frac{1}{2}(\frac{{(u-W)}^2}{{\sigma }_u^2}+\frac{v^2}{{\sigma }_v^2}))---\left(8\right).$

2)根据S57海图格式标准，从电子海图中提取代表河/海岸线的线条信息， 并转换为大地坐标系下的经纬度坐标；

3)通过灰色关联处理手段，针对雷达影像和电子海图所提取出来的岸线坐 标进行关联匹配，比较不同位置及尺度下的数据关联度，在这两种不同来源的 数据之间实现粗配准；

具体包括以下步骤：

(3.1)针对雷达影像和电子海图的特征，确定灰色关联分析中的参考数据 列和比较数据列：由于雷达影像为位图形式，而电子海图为矢量图形式，因此 将雷达影像中提取的岸线坐标设定为参考数据列，定义为：R＝{r(k)|k＝ 1,2,…,n，将电子海图中的岸线信息设定为比较数据列，分别定义为： E_i＝{e_i(k)|k＝1,2,…,n}。之后比较R和若干个E_i之间的几何形态相似程度，判 断其关联程度；

(3.2)雷达影像中的岸线信息为位图坐标，电子海图中的岸线信息为经纬 度信息，数据量纲不同，需根据公式(9)进行数据的无量纲化处理如下：

$e_i\left(k\right)=\frac{e_i\left(k\right)}{e_i\left(1\right)},r\left(k\right)=\frac{r\left(k\right)}{r\left(1\right)},i=\mathrm{1,2},...,m;k\mathrm{1,2},...,n---\left(9\right)$

(3.3)根据公式(10)和(11)，计算比较数据列R与参考数据列E_i中对应 数据的关联系数如下：

${\xi }_i\left(k\right)=\frac{{\min }_i\underset{k}{\min }{\Delta }_i\left(k\right)\rho {\max }_i{\max }_k{\Delta }_i\left(k\right)}{{\Delta }_i\left(k\right)+\rho {\max }_i{\max }_k{\Delta }_i\left(k\right)}---\left(10\right)$

Δ_i(k)＝|r(k)-e_i(k)| (11)

(3.4)取所有关联系数的平均值为数据列之间的关联度，并以此表示比较 数据列R与参考数据列E_i之间的关联程度，如公式(12)如示：

$r_i=\frac{1}{n}{\Sigma }_{k=1}^n{\xi }_i\left(k\right),k=\mathrm{1,2},...,n---\left(12\right)$

(3.5)将关联度r_i按大小排序，取关联度最大的比较数据列R与参考数据列 E_i为同一条岸线的粗配准结果。

4)通过ICP配准处理手段，针对雷达影像和电子海图所提取出来的岸线坐 标进行精配准；

具体包括以下步骤：

(4.1)将粗配准后的雷达影像岸线数据存储在点集P中，定义为： {P|p_i∈R³,i＝1,2,…,N}，p_i为岸线上某点的坐标，定义为：(x_pi,y_pi)；同时， 将对应的电子海图岸线数据存储在点集Q中，定义为：{Q|q_i∈R³,i＝1,2,…,N}， q_i为岸线上某点的坐标，定义为：(x_qi,y_qi)；作为同一段岸线的非同源数据，两 个点集之间应符合刚体变换，假定该刚体变换是由一个2*2的旋转矩阵R和一 个2*2的平移矩阵T实现的，如公式(13)所示：

$\left(\begin{array}{c}{x}_{\mathrm{pi}}\\ y_{\mathrm{pi}}\end{array}\right)=R\left(\begin{array}{c}{\lambda }_{\mathrm{qi}}\\ y_{\mathrm{pi}}\end{array}\right)+T---\left(13\right)$

(4.2)根据最优匹配原则设计数据配准的目标函数，如公式(14)所示， 点集R、T应使该目标函数的值最小：

$f(R,T)={\Sigma }_{i=1}^n{\left|\right|P_i^k-(R{Q}_i+T)\left|\right|}^2=\min ---\left(14\right)$

(4.3)计算点集P中的每一个点在Q点集中的对应最近点；

(4.4)计算上述对应点对平均距离最小的刚体变换，求得旋转矩阵R和平 移矩阵T；

(4.5)对Q使用上一步的平移矩阵T和旋转矩阵R，得到新的变换点集， 并替换获得新的点集Q；

(4.6)如果新的点集Q与参考点集P之间满足公式(14)中目标函数的要 求，即平均距离小于给定的阀值，则停止迭代，完成精配准；否则，以新的点 集Q与参考点集P继续迭代下去，直到达到目标函数的要求。

应当理解的是，对本领域普通技术人员来说，可以根据上述说明加以改进 或变换，而所有这些改进和变换都应属于本发明所附权利要求的保护范围。