CN106707287A

CN106707287A - 基于扩展卡尔曼滤波结合最近邻聚类算法的鱼群数量估计方法

Info

Publication number: CN106707287A
Application number: CN201611206595.1A
Authority: CN
Inventors: 韩军; 荆丹翔; 王杰英; 杜鹏飞; 章旻昊; 任佳
Original assignee: Zhejiang University ZJU
Current assignee: Zhejiang University ZJU
Priority date: 2016-12-23
Filing date: 2016-12-23
Publication date: 2017-05-24

Abstract

本发明公开了一种基于扩展卡尔曼滤波结合最近邻聚类算法的鱼群数量估计方法，其步骤：(1)固定双频识别声呐；(2)通过走航的方式进行水下探测，并获取声学数据；(3)数据处理，并利用扩展卡尔曼滤波结合最近邻聚类算法统计鱼体个数；(4)计算扫过水体体积，求取鱼群密度；(5)根据已知水域的蓄水量，估算整个水域中的鱼体数量。本发明利用双频识别声呐发射水下声波，遇到障碍物返回声波信号的原理，通过接收回波信号，并进行信号处理、获取水下鱼群信息，从而进行鱼群数量的估计，具有简单易行、高效而又不损坏鱼类资源等优点，且比传统的利用目标强度进行积分的方法具有更高的精度，为渔业资源评估提供了一种新途径。

Description

基于扩展卡尔曼滤波结合最近邻聚类算法的鱼群数量估计方法

技术领域

本发明属于渔业资源评估技术领域，尤其涉及一种基于扩展卡尔曼滤波结合最近邻聚类算法的鱼群数量估计方法。

背景技术

渔业资源评估是现代渔业开发过程中的重要环节，其中鱼群数量统计是渔业资源评估最基本的要求。传统方法主要依靠采样捕捞，这对鱼类资源本身具有损害；或者采用计量鱼探仪，利用回波积分法或回波计数法进行量测，这只能粗略估算鱼群数量，误差较大。现代社会对渔业资源的质量和产量提高、有效保护海洋生态系统，实现海洋资源的可持续发展提出了更高要求。例如应运而生的海洋牧场，就是传统捕捞和养殖方式向增殖和管理利用方式转变的一种现代渔业形式。在这样的情况下，如何快速准确，又不需要大量人力、物力消耗地高精度鱼群数量估计成为目前急需解决的问题。

发明内容

本发明的目的是针对上述不足，提供一种基于扩展卡尔曼滤波结合最近邻聚类算法的鱼群数量估计方法。

为了达到上述目的，本发明所采用的技术方案如下：基于扩展卡尔曼滤波结合最近邻聚类算法的鱼群数量估计方法，包括以下步骤：

(1)将双频识别声呐固定在调查船的船舷外侧或船底，并将其朝下浸没于水中，将姿态传感器固定在调查船上；

(2)将双频识别声呐通过网线与上位机相连，上位机实时获取走航时的声学数据，同时将姿态传感器接入上位机，上位机获取当前时刻调查船的姿态信息；

(3)对步骤(2)获得的声学数据进行图像构建及图像预处理，得到处理后的声学图像，通过目标检测算法识别图像中的鱼体目标，再通过扩展卡尔曼滤波结合最近邻聚类算法进行目标关联与跟踪，并统计声学图像中的鱼体数量；

(4)将GPS接入上位机，上位机记录下探测航迹，并结合双频识别声呐的探测范围，得到双频识别声呐扫过水体的体积，将步骤(3)获得的目标数量除以双频识别声呐扫过水体的体积，得到鱼群密度；

(5)根据已知水域的蓄水量，乘上步骤(4)中的鱼群密度，得到整个水域中鱼体数量。

进一步的，所述声学数据由多帧数据组成，每一帧数据代表一幅矩形声学图像。

进一步的，所述双频识别声呐布置于水面附近，深度不超过1米。

进一步的，所述双频识别声呐波束探测方向与水面夹角选取0°至90°之间任意角度。

进一步的，所述调查船的航速不大于6节。

进一步的，所述步骤(3)具体为：

(3.1)根据笛卡尔坐标系到极坐标系的坐标转换关系，将权利要求2中的矩形声学图像转换成扇形图，并通过插值将扇形图填充完整；

(3.2)将(3.1)得到的扇形图像进行线性拉伸处理；

(3.3)利用图像差分原理，将拉伸后的图像中的背景去除，留下图像中的有效目标，即水中的鱼体；

(3.4)利用扩展卡尔曼滤波结合最近邻聚类算法对图像中的有效目标进行跟踪处理，统计目标个数；

(3.5)将拉伸后的图像进行实时显示，并用不同颜色将每一个目标的轨迹显示出来。

进一步的，所述步骤(3.4)具体为：

(3.4.1)首先定义观测量为，其中r是声纳图像中有效目标相对于扇形对应的圆心的观测距离，α是对应的观测角度，系统的量测方程：

其中δr为距离值量测误差，δα为角度量测误差，(x,y)是有效目标在声纳图像中的位置，第i个有效目标在k时刻的系统状态向量为其中(x_k，i，y_k，i)是第i个有效目标在k时刻的位置，是第i个有效目标在k时刻的速度；

(3.4.2)计算第i个有效目标在前一个时刻状态量X_k-1(i)的预测值其中i＝1,2,…,m_k-1，m_k-1是k-1时刻的状态个数，预测方程为：

其中F是状态转移矩阵，Δt是采样时间间隔；

(3.4.3)分别计算预测值与当前时刻观测量Z_k(j)的权值w(i,j)，其中j＝1,2,…,n_k，n_k是当前时刻观测量的个数，权值就是与Z_k(j)的距离值；

(3.4.4)找到最小权值w_min＝min{w(i,j)}，以及对应的i和j，判断w_min是否小于设定的门限值ε，如果w_min<ε，将预测量与观测量Z_k(j)作为一对匹配值，最后在{w(i,j)}中去掉这对组合；

(3.4.5)重复步骤(3.4.2)、(3.4.4)m_k-1次；

(3.4.6)如果m_k-1个状态量中还有未匹配的状态，即表示此预测量对应的目标消亡了，目标总个数加一；如果n_k个观测量中还存在未匹配的量，即表示出现了新的目标，此目标的初始位置即为这个观测值，初始速度设为其他所有目标的平均速度。

(3.4.7)利用扩展卡尔曼滤波算法求得当前时刻的状态量X_k(i)，i＝1,2,…,m_k，其中m_k表示当前时刻状态量的个数；当所有有效目标经过N次状态更新后，就完成了有效目标数量的统计，其中N是所述声学数据的总帧数。

进一步的，所述利用扩展卡尔曼滤波算法求得当前时刻的状态量X_k(i)具体为：

(3.4.7.1)计算预测协方差矩阵：

P_k|k-1(i)＝FP_k-1(i)F^T+Γ_k(i)Q_k-1(i)Γ_k(i)^T

其中P_k-1(i)是前一个时刻的协方差矩阵，Q_k-1(i)是前一个时刻的系统过程噪声，Γ_k(i)是噪声协方差矩阵，上标T代表矩阵转置；

(3.4.7.2)计算卡尔曼增益矩阵：

其中R_k(i)是系统测量噪声。

(3.4.7.3)利用卡尔曼滤波方程，计算当前时刻的状态量X_k(i)：

h(k)是量测转移矩阵，Z_k(i)是当前时刻目标i对应的观测量。(3.4.7.4)更新当前时刻的协方差矩阵：

其中I表示单位阵。

本发明的有益效果如下：

1)本发明将双频识别声呐固定在船舷外侧或船底的水下，通过走航的方式采集水下声学数据，利用上位机软件将声学数据构建成声学图像，通过目标检测算法识别鱼体目标，再利用多目标跟踪算法进行目标关联并统计声学图像中的鱼体数量，结合GPS记录的走航轨迹计算双频识别声呐扫过的水体体积，获得鱼群密度，最后根据已知水域蓄水量估算出整片水体内的鱼群数量。

2)本发明中的上位机软件通过读取声学数据，构建声学图像，完成图像线性拉伸、目标提取以及利用扩展卡尔曼滤波结合最近邻聚类算法对图像中的有效目标进行跟踪处理并统计目标个数的任务，实现对双频识别声呐扫过水域内鱼体个数的自动计数。

3)本发明提供的方法简单易行、快速高效，不需要消耗大量人力物力，并且无损于水域中的渔业资源，相对于传统的捕捞估计或采用目标强度进行数量估计的方法具有更高精度，为渔业资源评估提供了一种新途径。

附图说明

图1是本发明方法的流程图；

图2是本发明中双频识别声呐布置示意图；

图3是本发明中设计的数据处理软件流程图；

图4是本发明中坐标转换示意图；

图5是本发明中多目标跟踪算法流程图；

图6是本发明中体积计算示意图；

图2中：GPS模块1、姿态传感器2、调查船3、双频识别声呐4、笔记本电脑5。

具体实施方式

下面结合具体实施例及附图对本发明做进一步详细说明，但本发明不仅局限于次。

本例中采用的双频识别声呐主要是有3个透镜和一个声呐阵列组成的多波束系统，它可以向水下发射频率为1.8MHz或1.1MHz的超声波，波束探测的范围最小是5米，最大是40米，接收数据的速率最大为每秒钟接收20帧图像，探测视场范围在水平方向为29°，垂直方向为14°，空气中重量约7千克，功率约30W。图1是利用双频识别声呐进行鱼群数量估计的实现流程图，主要过程阐述如下：

第一步，将双频识别声呐固定在调查船上。如图2所示，调查船以4节的均速沿着x轴运动，双频识别声呐固定在船舷外侧，它的正方向沿着y轴布置，并且波束发射方向与水面夹角设为60°，双频识别声呐浸没于水中并在水面附近，此例中距离水面约0.3米。

第二步，调查船在水面航行，双频识别声呐发射高频探测波束，接收回波信号，通过网线将双频识别声呐与上位机相连，上位机实时获取走航时的声学数据，同时将姿态传感器接入上位机，上位机获取当前时刻调查船的姿态信息；姿态传感器可以采用OCTANS光纤罗经(iXsea)，但不限于此，它用于调整水面晃动引起的声学数据波动。

第三步，对上一步中获得的声学数据进行图像构建及图像预处理，通过目标检测算法识别图像中的鱼体目标，再通过扩展卡尔曼滤波结合最近邻聚类算法进行目标关联与跟踪，并统计声学图像中的鱼体数量。

由于双频识别声呐传回的声学数据由多帧数据组成，每一帧数据代表一幅声学图像，为了实现鱼群数量统计，此例中利用Visual Studio 2015中的Microsoft FoundationClasses(MFC)模块设计了数据处理软件，软件具体实现流程如图3所示。

(1)读入一帧数据，由于每一帧声学数据是由一个96×512的二维数组构成，数组元素为0～255的灰度值。

(2)由于双频识别声呐的实际接收图像是扇形，因此需要将矩形声学图像转换成扇形图，如图4所示，矩形ABCD是原始数据，扇环A′B′C′D′是转换后的图形，以扇环对应的圆心O′建立坐标系，假设原图中一点P转化到扇环中的对应点是P′(x,y)，根据笛卡尔坐标系到极坐标系的转换公式：

其中ρ是O′P′的长度，θ是以O′为原点，x轴为极轴所对应的极角。经过坐标转换后，二维数组中的每个像素点被分配到了扇环的对应位置，由于这是一个非线性变换，因此扇环中有部分区域未被填充，尤其是扇环的顶部位置，因此，利用如下的公式进行差值处理：

其中v是需要插值的位置(x,y)对应的像素值，(x_i,y_i)是围绕在(x,y)四周的点，对应像素值为v_i，i的取值为1～4.

(3)将上步中得到的扇形图像进行线性拉伸处理，即灰度图像的增强处理，公式如下：

g(x,y)＝T[f(x,y)]＝af(x,y)+b

其中f(x,y)是原图像的灰度值，T(·)是映射函数，g(x,y)是增强后的图像，a是图像增强系数，b是偏移系数。a的取值为255除以图像中最大灰度值的商，b取值为0.

(4)利用图像差分原理，即设定图像阈值T，当像素值大于T时，该像素点判断为目标，否则是背景。将拉伸后的图像中的背景去除，留下图像中的有效目标，即水中的鱼体。

(5)利用扩展卡尔曼滤波结合最近邻聚类算法(NN‐EKF算法)对图像中的有效目标进行跟踪处理，统计目标个数；该算法具体实现过程如图5所示：

(a)首先定义观测量为：

其中r是声纳图像中有效目标相对于扇形对应的圆心的观测距离，α是对应的观测角度，系统的量测方程：

其中δr为距离值量测误差，δα为角度量测误差，(x,y)是有效目标在声纳图像中的位置，第i个有效目标在k时刻的系统状态向量为：

其中(x_k，i，y_k，i)是第i个有效目标在k时刻的位置，是第i个有效目标在k时刻的速度；

(b)计算第i个有效目标在前一个时刻状态量X_k-1(i)的预测值其中i＝1,2,…,m_k-1，m_k-1是k-1时刻的状态个数，预测方程为：

其中F是状态转移矩阵，Δt是采样时间间隔；

(c)分别计算预测值与当前时刻观测量Z_k(j)的权值w(i,j)，其中j＝1,2,…,n_k，n_k是当前时刻观测量的个数，权值就是与Z_k(j)的距离值；

(d)找到最小权值w_min＝min{w(i,j)}，以及对应的i和j，判断w_min是否小于设定的门限值ε，如果w_min<ε，将预测量与观测量Z_k(j)作为一对匹配值，最后在{w(i,j)}中去掉这对组合；

(e)重复步骤(c)、(d)m_k-1次；

(f)如果m_k-1个状态量中还有未匹配的状态，即表示此预测量对应的目标消亡了，目标总个数加一；如果n_k个观测量中还存在未匹配的量，即表示出现了新的目标，此目标的初始位置即为这个观测值，初始速度设为其他所有目标的平均速度。

(g)利用扩展卡尔曼滤波算法求得当前时刻的状态量X_k(i)，i＝1,2,…,m_k，其中m_k表示当前时刻状态量的个数。当所有有效目标经过N次状态更新后，就完成了有效目标数量的统计，其中N是所述声学数据的总帧数。

其中利用扩展卡尔曼滤波算法求得当前时刻的状态量X_k(i)具体过程为：

计算预测协方差矩阵：

P_k|k-1(i)＝FP_k-1(i)F^T+Γ_k(i)Q_k-1(i)Γ_k(i)^T

计算卡尔曼增益矩阵：

其中R_k(i)是系统测量噪声。

利用卡尔曼滤波方程，计算当前时刻的状态量X_k(i)：

h(k)是量测转移矩阵，Z_k(i)是当前时刻目标i对应的观测量。

更新当前时刻的协方差矩阵：

其中I表示单位阵。

(6)将拉伸后的图像进行实时显示，并用不同颜色将每一个目标的轨迹显示出来。

第四步，将GPS接入上位机，上位机记录下探测航迹，并结合双频识别声呐的探测范围，得到双频识别声呐扫过水体的体积，如图6所示，双频识别声呐探测范围是一个截面为扇形的五面体，探测的视角为α，此例中α为14°，探测的宽度为β，此例中β是29°，探测长度为h_i，根据水深可以自动调节；假设双频识别声呐从t₀时刻到t₀+Δt时刻移动了L_i的距离，则扫过的水体体积近似为：

公式中的L_i可以通过接入双频识别声呐的GPS模块计算得到。当航迹遍布调查水域的大部分区域时，利用步骤三中统计的鱼体个数除以总扫描体积V＝∑V_i，就得到鱼体密度ρ_fish，单位为尾/立方米。

利用上述原理，对上海的青草沙水库进行一次走航调查，走航总时长约为8小时，获得10GB声学数据，走航距离52.3公里，双频识别声呐扫过总体积为2.65×10⁵m³，自动计数总个数为3.2×10⁵尾，最后估算得到整个水库鱼体数量为3.6×10⁸尾。

Claims

1.一种基于扩展卡尔曼滤波结合最近邻聚类算法的鱼群数量估计方法，其特征在于，包括以下步骤：

(1)将双频识别声呐固定在调查船的船舷外侧或船底，并将其朝下浸没于水中，将姿态传感器固定在调查船上。

(2)将双频识别声呐通过网线与上位机相连，上位机实时获取走航时的声学数据，同时将姿态传感器接入上位机，上位机获取当前时刻调查船的姿态信息。

(3)对步骤(2)获得的声学数据进行图像构建及图像预处理，得到处理后的声学图像，通过目标检测算法识别图像中的鱼体目标，再通过扩展卡尔曼滤波结合最近邻聚类算法进行目标关联与跟踪，并统计声学图像中的鱼体数量。

(4)将GPS接入上位机，上位机记录下探测航迹，并结合双频识别声呐的探测范围，得到双频识别声呐扫过水体的体积，将步骤(3)获得的目标数量除以双频识别声呐扫过水体的体积，得到鱼群密度。

2.根据权利要求1所述的基于扩展卡尔曼滤波结合最近邻聚类算法的鱼群数量估计方法，其特征在于，所述声学数据由多帧数据组成，每一帧数据代表一幅矩形声学图像。

3.根据权利要求1所述的基于扩展卡尔曼滤波结合最近邻聚类算法的鱼群数量估计方法，其特征在于，所述双频识别声呐布置于水面附近，深度不超过1米。

4.根据权利要求1所述的基于扩展卡尔曼滤波结合最近邻聚类算法的鱼群数量估计方法，其特征在于，所述双频识别声呐波束探测方向与水面夹角选取0°至90°之间任意角度。

5.根据权利要求1所述的基于扩展卡尔曼滤波结合最近邻聚类算法的鱼群数量估计方法，其特征在于，所述调查船的航速不大于6节。

6.根据权利要求1所述的基于扩展卡尔曼滤波结合最近邻聚类算法的鱼群数量估计方法，其特征在于，所述步骤(3)具体为：

(3.2)将(3.1)得到的扇形图像进行线性拉伸处理；

(3.4)利用扩展卡尔曼滤波结合最近邻聚类算法对图像中的有效目标进行关联、跟踪处理，统计目标个数；

7.根据权利要求1所述的基于扩展卡尔曼滤波结合最近邻聚类算法的鱼群数量估计方法，其特征在于，所述步骤(3.4)具体为：

(3.4.1)首先定义观测量为其中r是声纳图像中有效目标相对于扇形对应的圆心的观测距离，α是对应的观测角度，系统的量测方程：

H : \{\begin{matrix} r = \sqrt{x^{2} + y^{2}} + δ r \\ α = t a n \frac{x}{y} + δ α \end{matrix}

其中δr为距离值量测误差，δα为角度量测误差，(x,y)是有效目标在声纳图像中的位置，第i个有效目标在k时刻的系统状态向量为其中(x_k,i，y_k,i)是第i个有效目标在k时刻的位置，是第i个有效目标在k时刻的速度；

{\hat{X}}_{k | k - 1} (i) = {FX}_{k - 1} (i)

其中F是状态转移矩阵，Δt是采样时间间隔；

(3.4.5)重复步骤(3.4.3)、(3.4.4)m_k-1次；

8.根据权利要求7所述的利用双频识别声呐进行鱼群数量估计的方法，其特征在于，所述利用扩展卡尔曼滤波算法求得当前时刻的状态量X_k(i)具体为：

(3.4.7.1)计算预测协方差矩阵：

P_k|k-1(i)＝FP_k-1(i)F^T+Γ_k(i)Q_k-1(i)Γ_k(i)^T

(3.4.7.2)计算卡尔曼增益矩阵：

K_{k} (i) = S_{k} (i) {(\frac{\partial H}{\partial X} S_{k} (i) + R_{k} (i))}^{- 1}

S_{k} (i) = P_{k | k - 1} (i) {(\frac{\partial H}{\partial X})}^{T}

其中R_k(i)是系统测量噪声。

(3.4.7.3)利用卡尔曼滤波方程，计算当前时刻的状态量X_k(i)：

X_{k} (i) = {\hat{X}}_{k | k - 1} (i) + K_{k} (i) (Z_{k} (i) - H {h (k) {\hat{X}}_{k | k - 1} (i)})

h(k)是量测转移矩阵，Z_k(i)是当前时刻目标i对应的观测量。

(3.4.7.4)更新当前时刻的协方差矩阵：

P_{k} (i) = (I - K_{k} (i) \frac{\partial H}{\partial X}) P_{k | k - 1} (i) + K_{k} (i) R_{k} (i) K_{k} {(i)}^{T}

其中I表示单位阵。