CN116797619A - 一种基于偏振图像的河道边界线检测方法 - Google Patents

Abstract

本发明公开了一种基于偏振图像的河道边界线检测方法，包括以下步骤：步骤1、通过偏振图像采集设备获取无人船前方河道场景图像数据；步骤2、对河道场景图像提取偏振特征图像，对水面偏振特征图像进行处理并且进行镜面反射、高光区域消除，得到降噪后的图像；步骤3、对步骤2中降噪后的图像进行处理获得河道边界线的拟合曲线方程。在阳光比较强的情况下，水面上有很多光线和反射；当太阳光从水面入射时在镜面方向形成强烈的反射辐射，形成水体的镜面反射。

Description

一种基于偏振图像的河道边界线检测方法

技术领域

本发明涉及无人船的智能感知领域，具体是一种基于偏振图像的河道边界线检测方法。

背景技术

近年来，无人船艇(USV)技术在执行巡逻、监测、安防等任务中发挥着越来越重要的作用。无人船在内陆河道需要检测河道边界线以确定水面可行驶区域。狭窄河道边界检测是无人船智能感知的重要组成部分。由于水面边界线未知，很难检测河道水域的可行驶区域。对于计算机来说河道没有明确的界限来区分河道边界线内的区域和河道边界线外的区域。

结构化河道指比较标准的河道，城市化结构特征较好的河道，这类河道具有清晰的水岸边界，河道的背景环境比较单一，河道边缘的几何特征也比较明显。因此，针对它的检测问题可以使用廉价的相机作为传感器，通过图像处理算法进行有效边界线的检测，可以为无人船自主巡航提供一种低成本的感知方案。

现有河道边界线检测技术中，普通的可见光相机工作在可见光波段，在雾天或者夜间与低照度环境下图像不清晰，此外，水面在晴天的阳光照射下反射光会非常强烈，相机大面积曝光导致获取的图像失效，这些条件都会导致图像处理算法难以实现河道的准确检测；采用激光雷达的方法中激光雷达价格较贵，且在雾霾条件受到强烈干扰会导致检测数据误差较大。因此常规的传感器不利于恶劣天气条件下的使用。

在实际应用中，USV获得的图像中存在大量噪声干扰，尤其是存在以下情况下：

(1)在阳光比较强的情况下，水面上有很多光线和反射；当太阳光从水面入射时在镜面方向形成强烈的反射辐射，形成水体的镜面反射；

(2)水面存在波浪，波浪顶端形成闪光点，形成大面积高光区域；

(3)河岸和水面之间的边界不清楚；

(4)天空的颜色与水面的颜色相似，仅根据位置判断该区域不可靠。

虽然目前存在针对上述极端行驶条件的图像处理方法，但是在现有的这些方法中，基于普通图像的处理方法计算过程复杂，难以满足实时性要求；基于偏振图像的处理方法多为特定场景，难以满足水面环境的要求；常规的河道边界检测方法准确度、鲁棒性低，难以满足无人感知系统的高精度要求。

发明内容

本发明的目的在于提供一种基于偏振图像的河道边界线检测方法，以解决上述背景技术中提出的问题。

为实现上述目的，本发明提供如下技术方案：

一种基于偏振图像的河道边界线检测方法，其特征在于，包括以下步骤：

步骤1、通过偏振图像采集设备获取无人船前方河道场景图像数据；

步骤2、提取不同角度的偏振特征图像，对水面偏振特征图像进行处理并且进行镜面反射、高光区域消除，得到降噪后的图像；

步骤3、对步骤2中降噪后的图像进行处理获得河道边界线的拟合曲线方程。

作为本发明进一步的方案：所述步骤2中镜面反射、高光区域消除包括以下步骤：

步骤2.1、获取偏振光角度在θ＝0,π/4,π/2,3π/4的4个偏振方向的偏振图像，

其中S₀、S₁、S₂是斯托克斯向量参数，θ是偏振角；

步骤2.2、通过时间域镜面反射抑制，生成时间域偏振辐射图，得到时间域高光消除结果I_time；选择偏振方向与偏振角相垂直的方向生成偏振辐射图I_⊥，取垂直方向对镜面反射光抑制作用最强的偏振辐射图

在成像环境和探测设备不变的情况下，反射光数量随时间变化波动不大，利用这个特性，将取得的I_⊥组成图像序列，进行时域融合，即将时间序列图像中对应像素的最小灰度值赋给最终的图像I_time(x,y)

I_time(x,y)＝min{I_m(x,y,t₁),I_m(x,y,t₂),…,I_m(x,y,t_nt)}

其中t₁,t₂,…,t_nt是连续的时间序列；

步骤2.3、通过空间域镜面反射抑制，生成空间域偏振辐射图，得到空间域高光消除结果I_space；

步骤2.4、将时空域镜面反射抑制图像融合，得到最终镜面反射消除结果的融合图像。

作为本发明进一步的方案：所述步骤2.3中的空间域高光消除处理包括以下步骤：

步骤2.3.1、在这4个方向偏振辐射图中选择对镜面反射光抑制作用最强的偏振辐射图，通过聚类算法在其中选择具有强光的区域R_k，

其中Ig_θ(x,y)是图像像素灰度值，δ_g是高光强度阈值；

步骤2.3.2、计算偏振光的线偏振度DOP和偏振角AOP，

其中S₀、S₁、S₂、S₃是斯托克斯向量参数S＝[S₀,S₁,S₂,S₃]，选用0°、45°、90°、135°四个偏振角度，用其所获得的光强信息来计算相关的偏振信息，其中，其中S₀表示光波总强度图像即I；S₁是水平方向和垂直方向线偏振光强度差；S₂是45°方向和135°两方向线偏振光强度差；S₃表示右旋与左旋圆偏振光强度差；

步骤2.3.3、分别在高光区域R_k(k＝1,2…m_r)对偏振度和偏振角图像作梯度特征融合：

步骤2.3.4、对镜面反射光求和得到空间域偏振图像处理结果

其中，m_r是高光区域R_k的个数。

作为本发明进一步的方案：所述步骤2.3.3中的梯度融合处理包括以下步骤：

步骤2.3.3.1、对偏振度图像和偏振角图像进行归一化操作；

步骤2.3.3.2、对偏振度图像和偏振角图像计算梯度特征；

两幅图像分别在R_k区域内的梯度特征计算如下所示

其中M，N分别代表图像的宽度与高度，ΔI_x(x,y)和ΔI_y(x,y)分别代表图像f(x,y)在x，y方向上的一阶差分，计算方法如下

ΔI_x(x,y)＝f(x,y)-f(x-1,y)

ΔI_y(x,y)＝f(x,y)-f(x,y-1)

步骤2.3.3.3、对偏振度图像以及偏振角图像进行加权融合，偏振度图像以及偏振角图像的融合公式如下

其中F_p,k表示高光区域R_k的偏振信息梯度特征融合后的图像，G_AOP和G_DOP分别为偏振角和偏振度图像对应的梯度特征，μ为权值系数，偏振角图像与偏振度图像相比含有的信息量更大，因此设置0.5<μ<1。

作为本发明进一步的方案：所述步骤2.4中的融合包括以下步骤：

步骤2.4.1选取I_time和I_space的R_k区域对应的子图像和/>进行融合，计算协方差

和协方差矩阵

其中n为子图像和/>的像素总数，/>为对应子图像的第i个像素，/>为对应子图像的像素均值，＊表示复共轭；

步骤2.4.2计算协方差矩阵C的特征值λ₁和λ₂和特征向量v₁和v₂，得到对角矩阵D和矩阵V，V是由特征向量组成的可逆矩阵；

步骤2.4.3取最大特征值对应的特征向量v＝[v₁,v₂]^T，计算子图像的加权系数ξ：

步骤2.4.4计算融合后的子图像

步骤2.4.5重复以上步骤，直到所有的高光区域都融合完成，再将这些处理后的子图像添加到原图像中得到完整的融合图像F(x,y)

式中表示仅仅子图像替换叠加，K_f是子图像的个数，f(x,y)是原偏振图像。

作为本发明进一步的方案：所述步骤3中的河道边界的拟合曲线的获取包括以下步骤：

步骤3.1、对图像进行预处理，降噪处理；

步骤3.2、通过双通道阈值异步并行提取河道边缘特征，并对双通道提取结果进行融合得到河道边界特征像素候选点C_cand；

步骤3.3、通过两阶段曲线段提取方式提取河道边界线特征，并拟合曲线方程。

作为本发明进一步的方案：所述步骤3.2双通道阈值异步并行是指边缘特征提取、图像分割特征两个通道异步并行计算河道边缘，所述河道边界特征像素候选点C_cand是将边缘特征提取、图像分割特征得到的二值图像边缘像素值合并，即：

C_cand＝C_edge∪C_seg；

其中：C_edge为边缘特征提取得到的候选点集合；C_seg为图像分割特征得到的候选点集合。

作为本发明进一步的方案：所述候选点集合C_edge是通过利用高斯可变向滤波器提取任意方向的河道边缘特征，得到候选点集合C_edge；

为了提取边缘信息，使用二阶高斯可变向滤波器检测任意可能方向的河道边缘，所述的方法基于高斯函数G(x,y)的方向导数算子是可变向控制的；

其中x，y是笛卡尔坐标，σ是高斯核的方差；

首先，通过计算高斯方向导数，得到基本滤波器；

河道边缘提取的二阶高斯可变向滤波器基于基本滤波器计算得到，基本滤波器通过x和y方向的方向导数计算：

其次，通过方向角将基本滤波器组合起来，得到能够计算任意方向的滤波响应，组合公式为：

其中θ代表方向角，在算法中应用两个可变向滤波器分别负责左右河道边缘特征的提取；

最后，为了突出河道边缘线特征，利用二值化方法设定阈值得到边缘提取结果：

得到候选点集合C_edge＝{(x,y)|I_edge(x,y)＝255}，其中G_thresh是阈值；

作为本发明进一步的方案：所述候选点集合C_seg首先通过高斯混合模型GMM进行图像分割，获得水面所在的区域图像I_water(u,v)；

GMM输入数据的分布是一组多元高斯分布的混合，概率密度函数为：

其中π_k是第k个高斯分量的权重，0≤π_k≤1且总和等于1，x代表数据点，K_G是分量总数，μ_k和∑_k是第k个分量多元高斯函数的均值和协方差参数；

每个分量都遵循多元正态分布，其概率密度函数形式为：

其中d是维度数，μ是d维均值向量，Σ是协方差矩阵，|Σ|是Σ的行列式；

GMM的特征向量是一个四维向量，其中包含红绿蓝(RGB)值和使非边缘噪声最小化的Sobel算子

x(k)＝[R(u,v),G(u,v),B(u,v),S(u,v)]

其中x是裁剪后的ROI图像的特征向量，k是x的索引，并且(u,v)是图像坐标；

其次，对图像分割所得河道和非河道区域确定提取水域边缘作为河道边界候选点集合C_seg；

选择I_water(u,v)最大的连接区域作为可能的河道区域，计算可以覆盖河道区域的凸包，然后找到并分离需要拟合的河道边界点区域C_seg；

Δi(u,v)＝I(u,v)-I(u-1,v)

得到候选点集合C_seg＝{(u,v)|I_bd(u,v)＝255}，其中W_thresh是阈值。

作为本发明进一步的方案：所述拟合曲线方程的获取包括以下步骤：

步骤3.3.1、第一阶段对C_cand边界点提取小曲线段集合；

通过聚类算法得到小线段集合Q，小线段集合由n_s个连续的边缘像素p(x,y)点组成，当两个相邻点p_i(x_i,y_i)和p_j(x_j,y_j)满足以下条件时认为是连续的，将其聚类在一起；

dist(p_i,p_j)＝|x_i-x_j|+y_i-y_j|<d₀

p_i(x_i,y_i)∈Q，p_j(x_j,y_j)∈Q，i≠j，

其中d₀是连续点距离阈值；

步骤3.3.2、对小曲线段内的边缘点进行二次曲线模型拟合；

将河道形状由二次曲线进行体现，河道边界线的二次方程定义如下

x＝f(y)＝a₀+a₁y+a₂y²

当a₂＝0，代表边界线模型为直线，对应河道模型为直线的场景；

将二次多项式拟合到m个点(x_i,y_i)，i＝1,2,…m，通过求解m个方程组和3个未知数a₀、a₁、a₂得到：

该模型求解过程可以写成一个矩阵方程：

AX＝B

定义模型的残差e_k用于求解最佳模型，定义如下：

求解模型的残差e_k最小化的过程，可以通过求解方程组计算：

(A^TA)X＝A^TB

考虑到矩阵A的特定形式，将矩阵代入方程通过显式执行矩阵乘法，求解方程简化为：

由上式可以看出求解的矩阵具有对称性，矩阵中元素是重复的，求解过程计算结果可以重复使用，从而可以降低计算成本，为了避免多项式的数值稳定性受到点集发散的影响，我们通过对点集归一化来提高稳定性；

步骤3.3.3、第二阶段再对小曲线段进行聚类，对小线段进行分组，将相同方向和共线的小线段分为同一组；

通过计算曲线相似性置信度，将相近且共线的小线段聚类在一起，具体通过系数差和线段之间的最短距离进行判断，计算公式如下：

Δd表示曲线l_i和l_j之间的最短距离，可通过离散采样的方法进行简化求解，Δc是系数之差，α、β表示Δd、Δc在判断相似性时的权重，λ是共线程度的阈值，代表曲线方程的系数；

步骤3.3.4、在上一步聚类的小线段结果中选择最优边界线，并确定河道左右边界；

左右边界必须满足如下关系

dist(L_left,L_right)≥T_d

len(L_left))≥T_len,len(L_right))≥T_len

nums(L_left,L_right)≥T_n

nums(L_left,L_right)＝argmax{NUM(cond{f′_left(y_k)·f′_right(y_k)<0})}

其中L_left和L_right是聚类为河道左右边界线区域小线段聚类结果，dist(L_left,L_right)是两条线段的距离，len(·)是线段的长度，T_d和T_len是距离和长度阈值，nums(L_left,L_right)代表同时满足左右边界水平方向对应点的切线方向相反(正负相反)的点的最大数量，cond{·}代表满足的条件关系，NUM(·)代表点的数量，argmax{·}代表满足条件的最大值；

步骤3.3.5、将上一步获得的左右边界小线段集合所在的点集合并，通过改进的RANSAC算法分别精确拟合左右边界线曲线方程，具体步骤如下：

步骤3.3.5.1、初始化河道边界线模型参数向量Param＝0，迭代次数k＝0，获取上一步河道边界特征点集P，内点集H＝NULL，最优参数向量best_Param＝NULL，最优分数best_score＝0，最大迭代次数iterations根据迭代次数上限值K设定，iterations＝INT(ξ*K)(ξ>1)，INT代表向上取整：

其中p_exp是期望迭代能够得到正确结果的概率，w是内点的占比，n_inliers和n_outliers分别是内点和外点的个数，m是河道边界线模型点的个数；

步骤3.3.5.2、从P中随机选择n个点，n＝3，并且第一个点是随机的，随后选择的点与已选择的点y值不重复，通过这些点计算Param＝[a₀,a₁,a₂]^T，并将这个n个点加入内点集H；

步骤3.3.5.3、对于任意点p_i(x_i,y_i)，p_i∈P，计算err(p_i)，如果err(p_i)<E_th，则将p_i加到内点集H上，重复执行此步骤直到所有点都判断完成，其中err(p_i)是误差：

err(p_i)＝|x_i-a₀-a₁y-a₂y_i ²|

步骤3.3.5.4、如果内点数量count(H)<T_th或Param不符合要求，H＝NULL，则转到步骤3.3.5.2，否则执行下一步；

步骤3.3.5.5、取内点集T的所有点，使用最小二乘法拟合重新计算参数Param，计算当前获得内点集得分score_k，如果score_k>best_score则best_score＝score_k,best_Param＝Param，H＝NULL，迭代次数k＝k+1，其中拟合得分计算如下：

步骤3.3.5.6、如果k<iterations重复步骤3.3.5.2-3.3.5.5，直到迭代次数满足要求。

与现有技术相比，本发明的有益效果是：

1、本发明通过对获取的河道场景图进行偏振处理，通过偏振处理完成对图片的镜面反射消除、高光区域消除，降低水面光反射对河道边界检测算法的影响；

2、本发明提出了基于时间域、空间域和梯度特征偏振图像融合方法，使得偏振特征图像提取对高光抑制作用最强，可以增强水面边界的特征信息，降低图像中无效信息干扰；

3、本发明提供了一种异步并行的双通道阈值方法提取河道边缘像素特征，使得边缘特征提取更加准确，对各种复杂河道环境场景的河道边缘像素提取具有鲁棒性，并且通过异步并行操作加快了处理速度；

4、本发明提供了一种基于两阶段曲线段提取方式提取河道边界线特征方程的方法，两个阶段的设计过滤掉了冗余和有噪声的信息，同时保留了关键特征信息，增强了鲁棒性，提高了边界线拟合精度。

附图说明

图1为本发明河道边界检测及评价方法的主要流程图；

图2为本发明通过偏振图像处理去除镜面反射引起的水面高光区域的流程图；

图3为本发明基于异步并行的双通道阈值法提取河道边缘点的流程图；

图4为本发明基于两阶段曲线段提取方式提取河道边界线曲线方程的流程图；

图5为本发明河道左右边界线上的点的切线方向示意图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

请参阅图1-5，本发明实施例中，一种基于偏振图像的河道边界线检测方法，包括以下步骤：

步骤1、获取无人船前方河道场景图像，本实施例通过安装于无人船正前方的偏振相机采集河道场景图像；

步骤2、对河道场景图像提取不同角度偏振特征图像，对水面偏振特征图像进行处理并且进行镜面反射、高光区域消除，得到降噪后的图像；

所述步骤2中镜面反射、高光区域消除包括以下步骤：

步骤2.1、获取偏振光角度在θ＝0,π/4,π/2,3π/4的4个偏振方向的偏振图像，

其中S₀、S₁、S₂是斯托克斯向量参数，θ是偏振角；

步骤2.2、通过时间域镜面反射抑制，生成时间域偏振辐射图，得到时间域高光消除结果I_time；选择偏振方向与偏振角相垂直的方向生成偏振辐射图I_⊥，取垂直方向对镜面反射光抑制作用最强的偏振辐射图

在成像环境和探测设备不变的情况下，反射光数量随时间变化波动不大，利用这个特性，将取得的I_⊥组成图像序列，进行时域融合，即将时间序列图像中对应像素的最小灰度值赋给最终的图像I_time(x,y)

I_time(x,y)＝min{I_m(x,y,t₁),I_m(x,y,t₂),…,I_m(x,y,t_nt)}

其中t₁,t₂,…,t_nt是连续的时间序列；

步骤2.3、通过空间域镜面反射抑制，生成空间域偏振辐射图，得到空间域高光消除结果I_space；

空间域高光消除处理包括以下步骤：

步骤2.3.1、在这4个方向偏振辐射图中选择对镜面反射光抑制作用最强的偏振辐射图，通过聚类算法在其中选择具有强光的区域R_k，

其中Ig_θ(x,y)是图像像素灰度值，δ_g是高光强度阈值；

步骤2.3.2、计算偏振光的线偏振度DOP和偏振角AOP，

其中S₀、S₁、S₂、S₃是斯托克斯向量参数S＝[S₀,S₁,S₂,S₃]，选用0°、45°、90°、135°四个偏振角度，用其所获得的光强信息来计算相关的偏振信息，其中，其中S₀表示光波总强度图像即I；S₁是水平方向和垂直方向线偏振光强度差；S₂是45°方向和135°两方向线偏振光强度差；S₃表示右旋与左旋圆偏振光强度差；

偏振相机传感器中的偏振镜是在传感器上每个透镜下通过线栅偏振镜微观实现的，偏振镜在四个像素组中具有0°、45°、90°和135°偏振角，分别对应相机输出图像的I₀、I₄₅、I₉₀和I₁₃₅，分别表示能透过的不同极化角度的偏振光；

偏振度(I_DOP)图像主要用来表示反射光中线偏振成分的多少，偏振角(I_AOP)图像主要用来表示辐射的两分量之间的相位差。I_DOP和I_AOP图像在将两者进行有效的融合，可以增强目标物体的特征信息，所以下一步基于两者的梯度特征做融合处理；

步骤2.3.3、分别在高光区域R_k(k＝1,2…m_r)对偏振度和偏振角图像作梯度特征融合：

梯度融合处理包括以下步骤：

步骤2.3.3.1、对偏振度图像和偏振角图像进行归一化操作；

在进行偏振信息图像的融合时，先将两者做归一化的处理，通过归一化可以使I_DOP图像与I_AOP图像处于同一范围；

步骤2.3.3.2、对偏振度图像和偏振角图像计算梯度特征；

两幅图像分别在R_k区域内的梯度特征计算如下所示

其中M，N分别代表图像的宽度与高度，ΔI_x(x,y)和ΔI_y(x,y)分别代表图像f(x,y)在x，y方向上的一阶差分，计算方法如下

Δi_x(x,y)＝f(x,y)-f(x-1,y)

ΔI_y(x,y)＝f(x,y)-f(x,y-1)

步骤2.3.3.3、对偏振度图像以及偏振角图像进行加权融合，平均梯度表达了图像细微部分的变化程度，可以用其表达图像的清晰程度。因此，通过利用梯度的特征把两者进行融合处理，偏振度图像以及偏振角图像的融合公式如下

其中F_p,k表示高光区域R_k的偏振信息梯度特征融合后的图像，G_AOP和G_DOP分别为偏振角和偏振度图像对应的梯度特征，μ为权值系数，偏振角图像与偏振度图像相比含有的信息量更大，因此设置0.5<μ<1。

步骤2.3.4、对镜面反射光求和得到空间域偏振图像处理结果

其中，m_r是高光区域R_k的个数

步骤2.4、将时空域镜面反射抑制图像融合，得到最终镜面反射消除结果，可有效抑制高光区域得到无饱和像素的融合图，并提高水岸边缘相对于水面背景的对比度。将I_time与I_space中分别在高光R_k对应的子区域内进行融合，融合步骤如下

步骤2.4.1选取I_time和I_space的R_k区域对应的子图像和/>进行融合，计算协方差

和协方差矩阵

其中n为子图像和/>的像素总数，/>为对应子图像的第i个像素，/>为对应子图像的像素均值，＊表示复共轭；

步骤2.4.2计算协方差矩阵C的特征值λ₁和λ₂和特征向量v₁和v₂，得到对角矩阵D和矩阵V，V是由特征向量组成的可逆矩阵；

步骤2.4.3取最大特征值对应的特征向量v＝[v₁,v₂]^T，计算子图像的加权系数ξ：

步骤2.4.4计算融合后的子图像

步骤2.4.5重复以上步骤，直到所有的高光区域都融合完成，再将这些处理后的子图像添加到原图像中得到完整的融合图像F(x,y)

式中表示仅仅子图像替换叠加，K_f是子图像的个数，f(x,y)是原偏振图像。

步骤3、对步骤2中降噪后的图像进行处理获得河道边界线的拟合曲线方程。

河道边界的拟合曲线的获取包括以下步骤：

步骤3.1、对图像进行预处理，降噪处理；

步骤3.2、通过双通道阈值异步并行提取河道边缘特征，并对双通道提取结果进行融合得到河道边界特征像素候选点C_cand；

所述步骤3.2双通道阈值异步并行是指边缘特征提取、图像分割特征两个通道异步并行计算河道边缘，所述河道边界特征像素候选点C_cand是将边缘特征提取、图像分割特征得到的二值图像边缘像素值合并，即：

C_cand＝C_edge∪C_seg；

其中：C_edge为边缘特征提取得到的候选点集合；C_seg为图像分割特征得到的候选点集合。

所述候选点集合C_edge是通过利用高斯可变向滤波器提取任意方向的河道边缘特征，得到候选点集合C_edge；

为了提取边缘信息，使用二阶高斯可变向滤波器检测任意可能方向的河道边缘，所述的方法基于高斯函数G(x,y)的方向导数算子是可变向控制的；

其中x，y是笛卡尔坐标，σ是高斯核的方差，高斯可变向滤波器由一组基本过滤器构成，具有可分离的特性，使得边缘检测更灵活、更快速；

首先，通过计算高斯方向导数，得到基本滤波器；

河道边缘提取的二阶高斯可变向滤波器基于基本滤波器计算得到，基本滤波器通过x和y方向的方向导数计算：

其次，通过方向角将基本滤波器组合起来，得到能够计算任意方向的滤波响应，组合公式为：

其中θ代表方向角，在算法中应用两个可变向滤波器分别负责左右河道边缘特征的提取；

最后，为了突出河道边缘线特征，利用二值化方法设定阈值得到边缘提取结果：

得到候选点集合C_edge＝{(x,y)|I_edge(x,y)＝255}，其中G_thresh是阈值；

所述候选点集合C_seg首先通过高斯混合模型GMM进行图像分割，获得水面所在的区域图像I_water(u,v)；

GMM输入数据的分布是一组多元高斯分布的混合，概率密度函数为：

其中π_k是第k个高斯分量的权重，0≤π_k≤1且总和等于1，x代表数据点，K_G是分量总数，μ_k和∑_k是第k个分量多元高斯函数的均值和协方差参数；

每个分量都遵循多元正态分布，其概率密度函数形式为：

其中d是维度数，μ是d维均值向量，Σ是协方差矩阵，|Σ|是Σ的行列式；

GMM的特征向量是一个四维向量，其中包含红绿蓝(RGB)值和使非边缘噪声最小化的Sobel算子

x(k)＝[R(u,v),G(u,v),B(u,v),S(u,v)]

其中x是裁剪后的ROI图像的特征向量，k是x的索引，并且(u,v)是图像坐标；

其次，对图像分割所得河道和非河道区域确定提取水域边缘作为河道边界候选点集合C_seg；

选择I_water(u,v)最大的连接区域作为可能的河道区域，计算可以覆盖河道区域的凸包，然后找到并分离需要拟合的河道边界点区域C_seg；

ΔI(u,v)＝I(u,v)-I(u-1,v)

得到候选点集合C_seg＝{(u,v)|I_bd(u,v)＝255}，其中W_thresh是阈值。

步骤3.3、通过两阶段曲线段提取方式提取河道边界线特征，并拟合曲线方程；

如图4所示，所述两阶段曲线段提取方式首先获得河道边界的小线段，其次应用聚类算法将小线段划分为不同的组。这两个阶段的设计是为了过滤掉冗余和有噪声的信息，同时保留关键信息。第二阶段将相同方向且共线的小线段聚类在一起，可以提高检测精度。

拟合曲线方程的获取包括以下步骤：

步骤3.3.1、第一阶段对C_cand边界点提取小曲线段集合；

通过聚类算法得到小线段集合Q，小线段集合由n_s个连续的边缘像素p(x,y)点组成，当两个相邻点p_i(x_i,y_i)和p_j(x_j,y_j)满足以下条件时认为是连续的，将其聚类在一起；

dist(p_i,p_j)＝|x_i-x_j|+|y_i-y_j|<d₀

p_i(x_i,y_i)∈Q，p_j(x_j,y_j)∈Q，i≠j，

其中d₀是距离阈值；

步骤3.3.2、对小曲线段内的边缘点进行二次曲线模型拟合；

将河道形状由二次曲线进行体现，河道边界线的二次方程定义如下

x＝f(y)＝a₀+a₁y+a₂y²

当a₂＝0，代表边界线模型为直线，对应河道模型为直线的场景；

将二次多项式拟合到m个点(x_i,y_i)，i＝1,2,…m，通过求解m个方程组和3个未知数a₀、a₁、a₂得到：

/>

该模型求解过程可以写成一个矩阵方程：

AX＝B

定义模型的残差e_k用于求解最佳模型，定义如下：

求解模型的残差e_k最小化的过程，可以通过求解方程组计算：

(A^TA)X＝A^TB

考虑到矩阵A的特定形式，将矩阵代入方程通过显式执行矩阵乘法，求解方程简化为：

由上式可以看出求解的矩阵具有对称性，矩阵中元素是重复的，求解过程计算结果可以重复使用，从而可以降低计算成本，为了避免多项式的数值稳定性受到点集发散的影响，我们通过对点集归一化来提高稳定性；

步骤3.3.3、第二阶段再对小曲线段进行聚类，对小线段进行分组，将相同方向和共线的小线段分为同一组；

通过计算曲线相似性置信度，将相近且共线的小线段聚类在一起，具体通过系数差和线段之间的最短距离进行判断，计算公式如下：

Δd表示曲线l_i和l_j之间的最短距离，可通过离散采样的方法进行简化求解，Δc是系数之差，α、β表示Δd、Δc在判断相似性时的权重，λ是共线程度的阈值，代表曲线方程的系数；

步骤3.3.4、在上一步聚类的小线段结果中选择最优边界线，并确定河道左右边界；

如图5所示，位于边界同一侧的小线段具有相同的方向，而位于相反河道边界的两个小线段具有相反的方向。我们可以通过计算集群中两个相反方向的点数量来筛选出最优的左右边界线。此外，我们选择位于道河道中间部分的区域，对于左边界线，我们选择最右边的边界，而对于右边界线，我们选择最左边的边界。且左右边界必须满足如下关系，左右边界必须满足如下关系

dist(L_left,L_right)≥T_d

len(L_left))≥T_len,en(L_right))≥T_len

nums(L_left,L_right)≥T_n

nums(L_left,L_right)＝argmax{NUM(cond{f′_left(y_k)·f′_right(yk)<0})}

其中L_left和L_right是聚类为河道左右边界线区域小线段聚类结果，dist(L_left,L_right)是两条线段的距离，len(·)是线段的长度，T_d和T_len是距离和长度阈值，nums(L_left,L_right)代表同时满足左右边界水平方向对应点的切线方向相反(正负相反)的点的最大数量，cond{·}代表满足的条件关系，NUM(·)代表点的数量，argmax{·}代表满足条件的最大值；

步骤3.3.5、将上一步获得的左右边界小线段集合所在的点集合并，通过改进的RANSAC算法分别精确拟合左右边界线曲线方程，具体步骤如下：

步骤3.3.5.1、初始化河道边界线模型参数向量Param＝0，迭代次数k＝0，获取上一步河道边界特征点集P，内点集H＝NULL，最优参数向量best_Param＝NULL，最优分数best_score＝0，最大迭代次数iterations根据迭代次数上限值K设定，iterations＝INT(ξ*K)(ξ>1)，INT代表向上取整：

其中p_exp是期望迭代能够得到正确结果的概率，w是内点的占比，n_inliers和n_outliers分别是内点和外点的个数，m是河道边界线模型点的个数；

步骤3.3.5.2、从P中随机选择n个点，n＝3，并且第一个点是随机的，随后选择的点与已选择的点y值不重复，通过这些点计算Param＝[a₀,a₁,a₂]^T，并将这个n个点加入内点集H；

步骤3.3.5.3、对于任意点p_i(x_i,y_i)，p_i∈P，计算err(p_i)，如果err(p_i)<E_th，则将p_i加到内点集H上，重复执行此步骤直到所有点都判断完成，其中err(p_i)是误差：

err(p_i)＝|x_i-a₀-a₁y-a₂y_i ²|

步骤3.3.5.4、如果内点数量count(H)<T_th或Param不符合要求，H＝NULL，则转到步骤3.3.5.2，否则执行下一步；

步骤3.3.5.5、取内点集T的所有点，使用最小二乘法拟合重新计算参数Param，计算当前获得内点集得分score_k，如果score_k>best_score则best_score＝score_k,best_Param＝Param，H＝NULL，迭代次数k＝k+1，其中拟合得分计算如下：

/>

步骤3.3.5.6、如果k<iterations重复步骤3.3.5.2-3.3.5.5，直到迭代次数满足要求。

对于本领域技术人员而言，显然本发明不限于上述示范性实施例的细节，而且在不背离本发明的精神或基本特征的情况下，能够以其他的具体形式实现本发明。因此，无论从哪一点来看，均应将实施例看作是示范性的，而且是非限制性的，本发明的范围由所附权利要求而不是上述说明限定，因此旨在将落在权利要求的等同要件的含义和范围内的所有变化囊括在本发明内。不应将权利要求中的任何附图标记视为限制所涉及的权利要求。

此外，应当理解，虽然本说明书按照实施方式加以描述，但并非每个实施方式仅包含一个独立的技术方案，说明书的这种叙述方式仅仅是为清楚起见，本领域技术人员应当将说明书作为一个整体，各实施例中的技术方案也可以经适当组合，形成本领域技术人员可以理解的其他实施方式。

Claims (10)

1.一种基于偏振图像的河道边界线检测方法，其特征在于，包括以下步骤：

步骤1、通过偏振图像采集设备获取无人船前方河道场景图像数据；

步骤2、提取不同角度的偏振特征图像，对水面偏振特征图像进行处理并且进行镜面反射、高光区域消除，得到降噪后的图像；

步骤3、对步骤2中降噪后的图像进行处理获得河道边界线的拟合曲线方程。

2.根据权利要求1所述的一种基于偏振图像的河道边界线检测方法，其特征在于，所述步骤2中镜面反射、高光区域消除包括以下步骤：

步骤2.1、获取偏振光角度在θ＝0,π/4,π/2,3π/4的4个偏振方向的偏振图像，

其中S₀、S₁、S₂是斯托克斯向量参数，θ是偏振角；

步骤2.2、通过时间域镜面反射抑制，生成时间域偏振辐射图，得到时间域高光消除结果I_time；选择偏振方向与偏振角相垂直的方向生成偏振辐射图I_⊥，取垂直方向对镜面反射光抑制作用最强的偏振辐射图

在成像环境和探测设备不变的情况下，反射光数量随时间变化波动不大，利用这个特性，将取得的I_⊥组成图像序列，进行时域融合，即将时间序列图像中对应像素的最小灰度值赋给最终的图像I_time(x,y)

I_time(x，y)＝min{I_m(x，y，t₁)，I_m(x，y，t₂)，…，I_m(x，y，t_nt)}

其中t₁,t₂,…,t_nt是连续的时间序列；

步骤2.3、通过空间域镜面反射抑制，生成空间域偏振辐射图，得到空间域高光消除结果I_space；

步骤2.4、将时空域镜面反射抑制图像融合，得到最终镜面反射消除结果的融合图像。

3.根据权利要求2所述的一种基于偏振图像的河道边界线检测方法，其特征在于，所述步骤2.3中的空间域高光消除处理包括以下步骤：

步骤2.3.1、在这4个方向偏振辐射图中选择对镜面反射光抑制作用最强的偏振辐射图，通过聚类算法在其中选择具有强光的区域R_k，

其中Ig_θ(,y)是图像像素灰度值，δ_g是高光强度阈值；

步骤2.3.2、计算偏振光的线偏振度DOP和偏振角AOP，

其中S₀、S₁、S₂、S₃是斯托克斯向量参数S＝[S₀,S₁,S₂,S₃]，选用0°、45°、90°、135°四个偏振角度，用其所获得的光强信息来计算相关的偏振信息，其中，其中S₀表示光波总强度图像即I；S₁是水平方向和垂直方向线偏振光强度差；S₂是45°方向和135°两方向线偏振光强度差；S₃表示右旋与左旋圆偏振光强度差；

步骤2.3.3、分别在高光区域R_k(k＝1,2…m_r)对偏振度和偏振角图像作梯度特征融合：

步骤2.3.4、对镜面反射光求和得到空间域偏振图像处理结果

其中，m_r是高光区域R_k的个数。

4.根据权利要求3所述的一种基于偏振图像的河道边界线检测方法，其特征在于，所述步骤2.3.3中的梯度融合处理包括以下步骤：

步骤2.3.3.1、对偏振度图像和偏振角图像进行归一化操作；

步骤2.3.3.2、对偏振度图像和偏振角图像计算梯度特征；

两幅图像分别在R_k区域内的梯度特征计算如下所示

其中M，N分别代表图像的宽度与高度，ΔI_x(x,y)和ΔI_y(x,y)分别代表图像f(x,y)在x，y方向上的一阶差分，计算方法如下

ΔI_x(x，y)＝f(x，y)-f(x-1，y)

ΔI_y(x，y)＝f(x，y)-f(x，y-1)

步骤2.3.3.3、对偏振度图像以及偏振角图像进行加权融合，偏振度图像以及偏振角图像的融合公式如下

其中F_p,k表示高光区域R_k的偏振信息梯度特征融合后的图像，G_AOP和G_DOP分别为偏振角和偏振度图像对应的梯度特征，μ为权值系数，偏振角图像与偏振度图像相比含有的信息量更大，因此设置0.5<μ<1。

5.根据权利要求2所述的一种基于偏振图像的河道边界线检测方法，其特征在于，所述步骤2.4中的融合包括以下步骤：

步骤2.4.1选取I_time和I_space的R_k区域对应的子图像和/>进行融合，计算协方差

和协方差矩阵

其中n为子图像和/>的像素总数，/>为对应子图像的第i个像素，为对应子图像的像素均值，＊表示复共轭；

步骤2.4.2计算协方差矩阵C的特征值λ₁和λ₂和特征向量v₁和v₂，得到对角矩阵D和矩阵V，V是由特征向量组成的可逆矩阵；

步骤2.4.3取最大特征值对应的特征向量v＝[v₁,v₂]^T，计算子图像的加权系数ξ：

步骤2.4.4计算融合后的子图像

步骤2.4.5重复以上步骤，直到所有的高光区域都融合完成，再将这些处理后的子图像添加到原图像中得到完整的融合图像F(x,y)

式中表示仅仅子图像替换叠加，K_f是子图像的个数，f(x,y)是原偏振图像。

6.根据权利要求1所述的一种基于偏振图像的河道边界线检测方法，其特征在于，所述步骤3中的河道边界的拟合曲线的获取包括以下步骤：

步骤3.1、对图像进行预处理，降噪处理；

步骤3.2、通过双通道阈值异步并行提取河道边缘特征，并对双通道提取结果进行融合得到河道边界特征像素候选点C_cand；

步骤3.3、通过两阶段曲线段提取方式提取河道边界线特征，并拟合曲线方程。

7.根据权利要求6所述的一种基于偏振图像的河道边界线检测方法，其特征在于，所述步骤3.2双通道阈值异步并行是指边缘特征提取、图像分割特征两个通道异步并行计算河道边缘，所述河道边界特征像素候选点C_cand是将边缘特征提取、图像分割特征得到的二值图像边缘像素值合并，即：

C_cand＝C_edge∪C_seg；

其中：C_edge为边缘特征提取得到的候选点集合；C_seg为图像分割特征得到的候选点集合。

8.根据权利要求7所述的一种基于偏振图像的河道边界线检测方法，其特征在于，所述候选点集合C_edge是通过利用高斯可变向滤波器提取任意方向的河道边缘特征，得到候选点集合C_edge；

为了提取边缘信息，使用二阶高斯可变向滤波器检测任意可能方向的河道边缘，所述的方法基于高斯函数G(x,y)的方向导数算子是可变向控制的；

其中x，y是笛卡尔坐标，σ是高斯核的方差；

首先，通过计算高斯方向导数，得到基本滤波器；

河道边缘提取的二阶高斯可变向滤波器基于基本滤波器计算得到，基本滤波器通过x和y方向的方向导数计算：

其次，通过方向角将基本滤波器组合起来，得到能够计算任意方向的滤波响应，组合公式为：

其中θ代表方向角，在算法中应用两个可变向滤波器分别负责左右河道边缘特征的提取；

最后，为了突出河道边缘线特征，利用二值化方法设定阈值得到边缘提取结果：

得到候选点集合C_edge＝{(,y)|_edge(x,y)＝255}，其中G_thresh是阈值。

9.根据权利要求7所述的一种基于偏振图像的河道边界线检测方法，其特征在于，所述候选点集合C_seg首先通过高斯混合模型GMM进行图像分割，获得水面所在的区域图像I_water(u,v)；

GMM输入数据的分布是一组多元高斯分布的混合，概率密度函数为：

其中π_k是第k个高斯分量的权重，0≤π_k≤1且总和等于1，x代表数据点，K_G是分量总数，μ_k和∑_k是第k个分量多元高斯函数的均值和协方差参数；

每个分量都遵循多元正态分布，其概率密度函数形式为：

其中d是维度数，μ是d维均值向量，Σ是协方差矩阵，|Σ|是Σ的行列式；

GMM的特征向量是一个四维向量，其中包含红绿蓝(RGB)值和使非边缘噪声最小化的Sobel算子

x(k)＝[R(u，v)，G(u，v)，B(u，v)，S(u，v)]

其中x是裁剪后的ROI图像的特征向量，k是x的索引，并且(u,v)是图像坐标；

其次，对图像分割所得河道和非河道区域确定提取水域边缘作为河道边界候选点集合C_seg；

选择I_water(u,v)最大的连接区域作为可能的河道区域，计算可以覆盖河道区域的凸包，然后找到并分离需要拟合的河道边界点区域C_seg；

ΔI(u，v)＝I(u，v)-I(u-1，v)

得到候选点集合C_seg＝{(u,v)|I_bd(u,v)＝255}，其中W_thresh是阈值。

10.根据权利要求6所述的一种基于偏振图像的河道边界线检测方法，其特征在于，所述拟合曲线方程的获取包括以下步骤：

步骤3.3.1、第一阶段对C_cand边界点提取小曲线段集合；

通过聚类算法得到小线段集合Q，小线段集合由n_s个连续的边缘像素p(x,y)点组成，当两个相邻点p_i(x_i,y_i)和p_j(x_j,y_j)满足以下条件时认为是连续的，将其聚类在一起；

dist(p_i，p_j)＝|x_i-x_j|+|y_i-y_j|＜d₀

p_i(x_i，y_i)∈Q，p_j(x_j，y_j)∈Q，i≠j，

其中d₀是连续点距离阈值；

步骤3.3.2、对小曲线段内的边缘点进行二次曲线模型拟合；

将河道形状由二次曲线进行体现，河道边界线的二次方程定义如下

x＝f(y)＝a₀+a₁y+a₂y²

当a₂＝0，代表边界线模型为直线，对应河道模型为直线的场景；

将二次多项式拟合到m个点(x_i,y_i)，i＝1,2,...m，通过求解m个方程组和3个未知数a₀、a₁、a₂得到：

该模型求解过程可以写成一个矩阵方程：

AX＝B

定义模型的残差e_k用于求解最佳模型，定义如下：

求解模型的残差e_k最小化的过程，可以通过求解方程组计算：

(A^TA)X＝A^TB

考虑到矩阵A的特定形式，将矩阵代入方程通过显式执行矩阵乘法，求解方程简化为：

由上式可以看出求解的矩阵具有对称性，矩阵中元素是重复的，求解过程计算结果可以重复使用，从而可以降低计算成本，为了避免多项式的数值稳定性受到点集发散的影响，我们通过对点集归一化来提高稳定性；

步骤3.3.3、第二阶段再对小曲线段进行聚类，对小线段进行分组，将相同方向和共线的小线段分为同一组；

通过计算曲线相似性置信度，将相近且共线的小线段聚类在一起，具体通过系数差和线段之间的最短距离进行判断，计算公式如下：

Δd表示曲线l_i和l_j之间的最短距离，可通过离散采样的方法进行简化求解，Δc是系数之差，α、β表示Δd、Δc在判断相似性时的权重，λ是共线程度的阈值，代表曲线方程的系数；

步骤3.3.4、在上一步聚类的小线段结果中选择最优边界线，并确定河道左右边界；

左右边界必须满足如下关系

dist(L_left，L_right)≥T_d

len(L_left))≥T_len，len(L_right))≥T_len

nums(L_left，L_right)≥T_n

nums(L_left，L_right)＝argmax{NUM(cond{f′_left(y_k)·f′_right(y_k)＜0})}

其中L_left和L_right是聚类为河道左右边界线区域小线段聚类结果，dist(L_left,L_right)是两条线段的距离，len(·)是线段的长度，T_d和T_len是距离和长度阈值，nums(L_left,L_right)代表同时满足左右边界水平方向对应点的切线方向相反(正负相反)的点的最大数量，cond{·}代表满足的条件关系，NUM(·)代表点的数量，argmax{·}代表满足条件的最大值；

步骤3.3.5、将上一步获得的左右边界小线段集合所在的点集合并，通过改进的RANSAC算法分别精确拟合左右边界线曲线方程，具体步骤如下：

步骤3.3.5.1、初始化河道边界线模型参数向量Param＝0，迭代次数k＝0，获取上一步河道边界特征点集P，内点集H＝NULL，最优参数向量best_Param＝NULL，最优分数best_score＝0，最大迭代次数iterations根据迭代次数上限值K设定，iterations＝INT(ξ*K)(ξ>1)，INT代表向上取整：

其中p_exp是期望迭代能够得到正确结果的概率，w是内点的占比，n_inliers和n_outliers分别是内点和外点的个数，m是河道边界线模型点的个数；

步骤3.3.5.2、从P中随机选择n个点，n＝3，并且第一个点是随机的，随后选择的点与已选择的点y值不重复，通过这些点计算Param＝[a₀,a₁,a₂]^T，并将这个n个点加入内点集H；

步骤3.3.5.3、对于任意点p_i(x_i,y_i)，p_i∈P，计算err(p_i)，如果err(p_i)<E_th，则将p_i加到内点集H上，重复执行此步骤直到所有点都判断完成，其中err(p_i)是误差：

err(p_i)＝|x_i-a₀-a₁y-a₂y_i ²|

步骤3.3.5.4、如果内点数量count(H)<T_th或Param不符合要求，H＝NULL，则转到步骤3.3.5.2，否则执行下一步；

步骤3.3.5.5、取内点集T的所有点，使用最小二乘法拟合重新计算参数Param，计算当前获得内点集得分score_k，如果score_k>best_score则best_score＝score_k,best_Param＝Param，H＝NULL，迭代次数k＝k+1，其中拟合得分计算如下：

步骤3.3.5.6、如果k<iterations重复步骤3.3.5.2-3.3.5.5，直到迭代次数满足要求。