形态学特征参数是大黄鱼（Pseudosciaena crocea Richardson）良种选育、种质鉴别及特征检测的重要依据^[1-2]。我国水产科技工作者已经围绕形态学特征参数对大黄鱼不同地理群体、不同育种群体、不同家系、雌雄性别、种质特征及鉴别比较等开展了深入研究^[3-7]，有效推动了我国大黄鱼人工繁殖与苗种培育进程。然而，目前在各类大黄鱼形态学研究过程中，一方面，形态参数测定依靠传统手工尺子测量，难以实现对大样本量研究对象的快速采集，且容易损坏活鱼样本^[8]；另一方面，所采用的鱼类形态学分析方法普遍限于线性回归分析、判别分析、聚类分析等传统多元统计分析方法^[9]，在鱼类形态特征挖掘及群体识别模型构建方法方面，至今尚缺乏有效的手段。因此，研究建立一套高效的鱼类形态参数测量、形态特征挖掘和群体识别方法，对于大黄鱼良种选育、种质鉴别等具有重要意义。

本文以岱衢族大黄鱼F₂、F₃代选育群体为研究对象，通过计算机视觉（computer vision, CV）^[10]测定2类选育群体的形态参数，使用主成分分析（principal component analysis, PCA）^[11]和连续投影算法(successive projections algorithm, SPA) ^[12]对选育群体的形态参数进行分析，并获得能表达不同选育群体之间形态差异的特征变量。在此基础上，以特征变量为输入，结合稀疏表示（sparse representation, SR）^[13]方法建立岱衢族大黄鱼选育群体分类识别模型，实现不同选育群体的自动识别。

1 材料与方法 1.1 实验材料

2015年11月、2016年11月在浙江省宁波市象山港湾水产苗种有限公司育苗场（象山县黄避岙乡高泥村凤凰礁），分批采集普通网箱（5 m×10 m）养殖、投喂冰鲜鱼饵料的15月龄岱衢族大黄鱼F₂、F₃代选育群体样本各600尾，用于大黄鱼形态参数的计算机视觉测量、形态特征挖掘和选育群体识别建模研究。

1.2 实验方法 1.2.1 形态参数的计算机视觉测量 1.2.1.1 形态参数测量系统

采用自己研制的鱼类形态参数计算机视觉测量系统，对大黄鱼形态参数和质量进行测量。该系统硬件结构和实物如图 1所示，主要包括框架主体、LED光源（YX-BL25844，山东省烟台永信视觉技术有限公司）、框架主体上的透明玻璃平台、设置在透明玻璃平台四角下端的称重传感器（KNX/3kg，江苏省苏州坤宏电子有限公司）、放置在透明玻璃平台下端用于拍摄鱼体图像的工业相机（SJM300C，台湾视觉美深圳有限公司）和镜头（VM06012MP，广州智赛电子有限公司）。工业相机与计算机信号连接，采用下置方式采集透明玻璃平台上被测样本鱼体图像，镜头与被测样本距离保持固定，从而避免不同样本鱼体厚度差异所带来的尺寸测量误差。称重传感器通过称重控制器（XK3190-A27E，上海耀华称重系统有限公司）与计算机信号连接，同步传输鱼体质量数据。

用安装在计算机中的鱼类形态参数测量软件（图 2）采集被测样本的数字图像和重量数据，并利用数字图像对鱼体形态参数进行精确测量分析。该软件是在Microsoft Visual Studio 2012开发环境下采用C#语言并基于开放视觉图形处理库Emgu CV编写，其主要功能包括：1）数据采集管理。用于采集和管理被测样本的数字图像、质量和形态参数数据。2）参数设置。用于设定工业相机的图像分辨率、曝光时间、白平衡等参数；利用棋盘格标定板进行相机标定（消除图像畸变）和像素当量标定（即图像中一个像素点代表的实际物理尺寸）。3）形态参数测量。用于测量鱼体尺寸、外形轮廓曲率等参数。

本研究中鱼类形态参数计算机视觉测量系统的设计与开发细节、测量精度验证等相关内容参考余心杰等^[14]的报道。

1.2.1.2 形态参数采集与处理

利用上述鱼类形态学参数的计算机视觉测量系统对大黄鱼样本的形态参数进行快速活体测量。先将活鱼样本分批次（每组20尾）放入质量浓度为15 mg/L的丁香酚水溶液中，麻醉1~2 min至鱼体侧翻不运动，逐个取出样本放到透明玻璃平台上，通过鱼类形态参数测量软件，采集鱼体数字图像和质量数据，并在计算机中保存。完成全部样本数字图像和质量数据采集之后，在软件中采用手工鼠标标记方式，对鱼体形态参数进行批量测定。

1.2.2 形态特征挖掘

在数据挖掘领域，特征提取（feature extraction, FE）和特征选择(feature selection, FS)是从原始特征中找出本质特征的有效手段，其中特征提取是指由原始特征经过线性或非线性变换得到较少数量但具有表达能力的新特征，特征选择是指直接从原始特征中选取最相关的特征子集。目前，在鱼类形态学研究中，主要通过对样本鱼体形态性状参数之间的相关性关系和通径分析，寻找出形态特征（显著相关的形态性状参数或参数比例）作为生长预测、种群划分等的依据，本质上还是属于特征选择范畴。为了探索建立鱼类形态特征挖掘方法，本文分别从基于主成分分析（PCA）的特征提取和基于连续投影算法（SPA）的特征选择2个角度，开展岱衢族大黄鱼不同选育群体的形态特征挖掘研究。

1.2.2.1 基于PCA的形态特征提取

PCA是将多个相关的特征变量转化为少数几个相互独立的新特征变量的有效分析方法。假设X为采集到的一个n×m的大黄鱼样本数据矩阵，其中的每一列对应一个特征变量，每一行对应一个样本，则X∈R^n×m。先将X标准化为：

$ \mathit{\boldsymbol{X}}\prime = [\mathit{\boldsymbol{X}} - {\left( {1\;\;1{\rm{ }}\;\; \ldots \;{\rm{ }}1} \right)^{\rm{T}}}M]{\rm{diag}}(\frac{1}{{{s_1}}}, \frac{1}{{{s_2}}}, \ldots , \frac{1}{{{s_m}}}). $

(1)

其中：M=[m₁ m₂ … m_m]，为X的均值；[s₁ s₂，…，s_m]为X的标准差。

X′的PCA模型可表示为：

$ \mathit{\boldsymbol{X}}\prime = {t_1}p_1^{\rm{T}} + {t_2}p_2^{\rm{T}} + {\rm{ }} \ldots {\rm{ }} + {t_m}p_m^{\rm{T}}. $

(2)

式中：t₁，t₂，…，t_m∈Rⁿ，为得分向量，也叫主元；p₁，p₂，…，p_m∈R^m，为载荷向量。如果可以用前k（k＜m，一般取累计贡献率大于85%对应的主元个数）个主元来代替数据中的主要变化，那么X′的PCA模型可表示为：

$ \mathit{\boldsymbol{X}}\prime = {t_1}p_1^{\rm{T}} + {t_2}p_2^{\rm{T}} + {\rm{ }} \ldots {\rm{ }} + {t_k}p_k^{\rm{T}} + E = {\mathit{\boldsymbol{X}}_k}\prime {\rm{ }} + E. $

(3)

其中E为误差矩阵。X′可以近似表示为：

$ {\mathit{\boldsymbol{X}}_k}\prime = {t_1}p_1^{\rm{T}} + {t_2}p_2^{\rm{T}} + {\rm{ }} \ldots {\rm{ }} + {t_k}p_k^{\rm{T}}. $

(4)

本文采用主成分分析对大黄鱼样本数据矩阵X进行变换，用变换后的前k（k＜24）个主元代替原始24个特征变量信息，实现鱼体形态特征提取。

1.2.2.2 基于SPA的形态特征选择

SPA的目标是选取彼此之间共性较少的变量组合，且尽可能地保留原始数据信息，从而达到特征选择的目的。SPA的工作原理是通过迭代的方法，从一个变量开始，每次循环都计算其在未选入变量上的投影，将投影向量最大的变量引入到变量组合中，每个新选入的变量，都与前一个的线性关系最小。对于由样本数n和特征变量数m组成的大黄鱼样本数据矩阵X∈R^n×m，k（k＜m－1）为需要选择的特征变量个数，SPA的特征选择步骤如下：

1）利用式（1），对样本数据矩阵X进行标准化处理，得到X′；

2）第一次迭代（p=1）开始前，在X′中任选一列向量x′_j，记为x′_K(0)，即K(0)=j，j∈1, …, m；

3）将未选入的列向量的位置计为s，s ={j, 1≤j≤ m, j $ \notin $ {K(0), …, K(p－1)}}；

4）计算剩余列向量x′_j(j∈s)与当前所选向量x′_K(p－1)的投影：

$ {P_{x{\prime _j}}} = x{\prime _j} - [x\prime _j^{\rm{T}}x{\prime _{K(p - 1)}}]x{\prime _{K(p - 1)}}{[x\prime _{K(p - 1)}^{\rm{T}}x{\prime _{K(p - 1)}}]^{ - 1}}; $

(5)

5）提取投影向量最大特征变量的列向量位置：

$ K(p) = {\rm{arg}}[{\rm{max}}(\parallel {P_{x{\prime _j}}}\parallel )], j \in s; $

(6)

6）令$ {{x'}_j} = {P_{{{x'}_j}}}, j \in s;p = p + 1$。如果p＜k，返回到第3）步循环计算。

对于每一个初始K(0)，循环一次后利用多元线性回归（multiple linear regression, MLR）进行交互验证分析，最小均方根误差（root mean square error, RMSE）对应的K(p)即为特征选择结果。本文采用SPA对大黄鱼样本数据矩阵X中的24个特征变量进行选择，最终选取k个特征变量组合{K(p)，p=0，…，k －1}来表达大黄鱼样本形态差异。

1.2.3 稀疏表示方法

对大黄鱼不同选育群体的形态特征提取和选择之后，不同选育群体识别的基本问题就是基于形态特征数据，根据L类训练样本，确定测试样本属于哪一类。目前常用的模式识别方法是基于反向传播（back propagation, BP）神经网络或支持向量机（support vector machine, SVM）等建立相应的分类模型，将测试样本的特征数据输入该模型，得到相应的识别结果^[15]。与上述方法需通过多个参数优化选择过程来建立一个分类模型不同，稀疏表示法直接通过测试样本在每个训练样本上的投影系数进行识别分类^[16]，识别过程简单，易实现。

假设大黄鱼选育群体种类有L类；每一类的训练样本数分别为N₁，N₂，…，N_i，…，N_L；每个训练样本的特征变量为m维，用列向量f来表示，则f∈R^m；同一类训练样本数据位于一个线性子空间中，那么新测试样本数据可看成是同类样本数据的线性组合。

将充足的训练样本给第i个类，然后将第i类中的N_i个训练样本数据组成矩阵A_i（称之为训练样本矩阵）的列，具体表示为：

$ {\mathit{\boldsymbol{A}}_i} = [{f_{i, 1}}, {f_{i, 2}}, \ldots , {f_{i, {N_i}}}]{\rm{ }} \in {R^{m \times {N_i}}}. $

(7)

任意来自于相同类的测试样本y∈R^m，近似满足于以下由训练样本组成的线性组合：

$ y = {a_{i, 1}}{f_{i, 1}} + {a_{i, 2}}{f_{i, 2}} + {\rm{ }} \ldots {\rm{ }} + {a_{i, {N_i}}}{f_{i, {N_i}}}. $

(8)

随后，为所有类的训练样本集定义一个新的训练样本矩阵A作为所有L类训练样本的集合：

$ \mathit{\boldsymbol{A}} = [{\mathit{\boldsymbol{A}}_1}, {\mathit{\boldsymbol{A}}_2}, \ldots , {\mathit{\boldsymbol{A}}_L}] = [{f_{1, 1}}, {f_{1, 2}}, \ldots , {f_{L, {N_L}}}]. $

(9)

则y的线性表示可以根据所有 L类训练样本写成另一种形式：

$ y = \mathit{\boldsymbol{A}}x. $

(10)

式中：列向量$ x = {[0, 0, \ldots , {a_{i, 1}}, {a_{i, 2}}, \ldots , {a_{i, j}}, \ldots , {a_{i, {N_i}}}, 0, 0, \ldots , 0]^{\rm{T}}}$是一个系数向量；元素a_{i, j}表示y在A中第i类的第j个训练样本上的投影系数。

若y来自第i类样本，则x在其所属类别的A_i上的投影系数不为0，而在其他训练样本类别的A_l (l≠ i)上的投影系数为0。也就是说，y的解是稀疏的。

稀疏表示和压缩感知的理论研究表明，如果解x是稀疏的，则可通过解决如下的L-1最小化范数问题来求解式（10）：

$ \hat x = {\rm{ arg}}\;\;{\rm{ min}}\parallel x{\parallel _1}{\rm{subject}}\;\;{\rm{ to}}\;\mathit{\boldsymbol{A}}x = y. $

(11)

式中：${\hat x} $为x的近似解；||·||₁表示向量中各元素的绝对值之和。

在理想情况下，${\hat x} $的非零系数中最大的系数值如果对应第i个类，且较大的系数基本都同第i个类有关，就可将y归到那个类中。但是，在实际应用中噪声等问题可能会导致${\hat x} $中有一部分与其他类有关的小的非零系数项出现。因此，对于每个类i，可定义一个函数δ_i，用来在稀疏解${\hat x} $中选取与第i类有关的系数。δ_i(${\hat x} $)是一个新的向量，其中除了与第i类有关的系数外，令其他的系数值都为零，重新生成给定的测试样本y的一个近似表示，将其记为：

$ {{\hat y}_i} = \mathit{\boldsymbol{A}}{\delta _i}(\hat x{\rm{ }}). $

(12)

显然，${{\hat y}_i} $与y距离越小，则${{\hat y}_i} $属于第i类的可能性越高。因此，可通过计算${{\hat y}_i} $与y的最小残差值，即求解以下方程

$ \mathop {{\rm{min}}}\limits_i {r_i}(y){\rm{ }} = {\rm{ }}\mathop {{\rm{min}}}\limits_i \parallel y - \mathit{\boldsymbol{A}}{\delta _i}(\hat x){\parallel _2} $

(13)

得到残差最小的i对应的类别，即为y的类别标志。

1.2.4 基于形态特征和稀疏表示的大黄鱼选育群体识别方法

通过上述方法，建立完整的基于形态特征和稀疏表示的岱衢族大黄鱼F₂、F₃代选育群体识别方法的步骤如下。

1）形态特征挖掘：分别采用PCA和SPA方法，对所采集的岱衢族大黄鱼F₂、F₃代选育群体的24个形态参数进行特征提取和选择，得到k个特征变量。

2）输入：在每一种选育群体中随机选取300个样本作为训练样本，剩余的300个样本作为测试样本。利用所有训练样本的k个特征变量，构成如下训练样本矩阵：

$ \mathit{\boldsymbol{A}} = [{\mathit{\boldsymbol{A}}_{{{\rm{F}}_{\rm{2}}}}}, {\mathit{\boldsymbol{A}}_{{{\rm{F}}_{\rm{3}}}}}\left] {{\rm{ }} = } \right[{f_{{{\rm{F}}_{\rm{2}}}, 1}}, \ldots , {f_{{{\rm{F}}_{\rm{2}}}, 300}}, {f_{{{\rm{F}}_{\rm{3}}}, 1}}, \ldots , {f_{i, j}}, \ldots , {f_{{{\rm{F}}_{\rm{3}}}, 300}}]. $

(14)

式中：A_F2、A_F3表示岱衢族大黄鱼F₂、F₃代选育群体的训练样本矩阵；f_{i, j}为第i种选育群体的第j个训练样本的k个特征变量组成的向量列。矩阵A为F₂、F₃代选育群体的全体训练样本矩阵，其大小为k× 600。选取某个选育群体中的测试样本y。

3）L-1范数最小化求解：利用式(11)，得到系数向量${\hat x} $。这里的计算误差阈值设为0.01，迭代次数为150次。

4）计算残差值：利用式(13)，计算F₂、F₃代选育群体2个类的${{\hat y}_i} $与y的最小残差值。

5）输出：具有最小残差的类作为判定测试样本y的类别。

上述PCA算法、SPA算法、稀疏表示算法在Matlab 7.9下编程实现，其中稀疏表示算法中的L-1范数最小化方法采用美国斯坦福大学研究人员编写的Matlab软件包^[17]实现；计算机硬件环境为Intel(R) Core(TM) i5-2300 CPU、主频2.80 GHz、内存4 096 MB的PC机，软件环境为Windows 7.0操作系统。

2 结果与讨论 2.1 形态测量结果

利用鱼体形态参数的计算机视觉测量系统采集大黄鱼活鱼数字图像和体质量的速度约150尾/h，满足大样本量鱼类快速、活体采集要求。实验共采集岱衢族大黄鱼F₂、F₃代选育群体的数字图像和体质量样本1 200个（F₂、F₃代各600尾），测定了体质量、体长、体宽、头部曲率、全长、全高、吻长、尾柄长、尾柄高、头长等10个参数。2个群体样本都为15月龄，样本的体质量、全长等规格相近（表 1）。计算相关形态参数比值得到的24个形态参数描述见表 2。

2类选育群体样本的外形如图 3所示。根据样本外形观测，F₃代选育群体与F₂代选育群体外形具有一定差异，这是由于在岱衢族大黄鱼选育过程中进行了外形淘汰选育，F₃代群体样本的总体体型偏长，尤其是尾柄形态较细长。

2.2 形态特征挖掘 2.2.1 主成分分析特征提取

对岱衢族大黄鱼F₂、F₃代选育群体样本的24个形态参数进行主成分分析，得到前10个主成分（PC₁~PC₁₀）的累计方差贡献率（表 3）和2个群体的24个形态参数的主成分载荷量分布（图 4）。

由表 3和图 4可知：第1主成分PC1贡献率为58.38%，对其载荷量较大的是全长/尾柄长（TL/ CPL）、体长/尾柄长（BL/CPL）、尾柄长/尾柄高（CPL/CPH）、体高/体质量（BH/BM）；第2主成分PC₂贡献率为21.87%，对其载荷量较大的是全长/头长（TL/HL）、体长/头长（BL/HL）；第3主成分PC₃贡献率为8.94%，对其载荷量较大的是全长/体长（TL/ BL）、体长/吻长（BL/KL）。前3个主成分的累计贡献率达到89.19%，说明PCA提取得到的这3个主成分因子已经能够较好地表达2类选育群体的形态参数信息。基于PCA的载荷分析表明，上述对第1、2、3主成分取值影响大的8个形态参数，可以被认为是选育群体样本之间形态差异的重要因子。根据2类选育群体样本的第1、2、3主成分绘制散点图(图 5)，从中可以看出，虽然F₂、F₃代选育群体的部分样本有重叠交错，但是这2类群体样本整体上还是可以区分的。

总之，由于岱衢族大黄鱼F₂、F₃代选育群体形态参数具有差异，可以利用PCA提取的主成分或者利用PCA载荷分析得到的8个特征形态参数，对不同选育群体样本进行识别。

2.2.2 连续投影算法特征选择

利用SPA从岱衢族大黄鱼F₂、F₃代选育群体样本的24个形态参数中选择共线性最小的有效特征形态参数。经SPA选择的不同变量个数的最小均方根误差（RMSE）分布图如图 6所示。从中可以看出，提取出7个特征形态参数时RMSE达到低点，此后RMSE趋于平稳，因此，可以从原始24个形态参数中选择这7个特征形态参数，它们分别是全长/体长（TL/BL）、全长/头长（TL/HL）、全长/尾柄长（TL/ CPL）、体长/头长（BL/HL）、尾柄长/尾柄高（CPL/ CPH）、全长/尾柄高（TL/CPH）、体长/尾柄高（BL/ CPH）。

通过比较PCA载荷分析得到的8个特征形态参数和SPA选择的7个特征形态参数之后发现，2种特征挖掘方法都筛选出了全长/体长（TL/BL）、全长/头长（TL/HL）、全长/尾柄长（TL/CPL）、体长/头长（BL/HL）、尾柄长/尾柄高（CPL/CPH）这5个特征形态参数。从形态参数比值含义来看，5个特征形态参数主要反映了大黄鱼头部、尾鳍、尾柄部位长度在鱼体长中的占比，集中表现在全长、体长、头长、尾柄长4个指标，这与直接观察所得到的F₃代群体较F₂代群体样本总体体型偏长、尾柄形态细长的形态差异的结论一致。

2.3 岱衢族大黄鱼F₂、F₃代选育群体识别

分别利用PCA提取的前10个主成分（PCA变换主元特征）、PCA载荷分析得到的8个特征形态参数组合（PCA选择特征）及SPA选择的7个特征形态参数组合（SPA选择特征），按照1.2.4节稀疏表示识别步骤，对岱衢族大黄鱼F₂代（DQF₂）、F₃代（DQF₃）选育群体样本进行识别。不同特征对选育群体的识别结果如表 4所示。可以看出，3组特征对岱衢族大黄鱼F₂代（DQF₂）、F₃代（DQF₃）选育群体样本进行识别的平均准确率都在75%以上，且岱衢族大黄鱼F₂代、F₃代群体之间存在较显著的形态差异，能够进行形态区分。其中：利用PCA提取的前10个主成分对2类选育群体样本进行识别，得到最高的平均识别准确率为88.3%。这是由于主成分分析能够在丢失最少信息的条件下，将原有众多变量浓缩成较少的能够反映原有变量绝大部分信息的综合指标，因此能够全面反映群体间的综合形态差异，得到较好的识别效果。利用PCA载荷分析选择的8个特征形态参数组合和利用SPA选择的7个特征形态参数组合对2类选育群体样本的识别准确率分别为79.0%、80.5%，两者识别效果没有明显差距，但是相比PCA提取的前10个主成分，其识别效果不够理想。究其原因，可能是岱衢族大黄鱼F₂、F₃代选育群体代系之间虽然存在形态差异，但是形态多样性较低，难以直接使用较少的形态参数因子来概括群体之间的形态差异。因此，在岱衢族大黄鱼选育群体代系之间的形态差异描述中要充分利用数据挖掘手段，采用特征提取方法获得综合性形态特征指标来对其形态差异进行描述，从而克服在传统鱼类形态学研究中传统特征选择方法所获得的单一性形态特征指标不全面等缺点。

3 结论

本文采用计算机视觉技术测定大黄鱼形态参数，通过主成分分析和连续投影算法提取形态特征变量，并结合稀疏表示方法实现了对岱衢族大黄鱼F₂、F₃代选育群体的识别建模，为鱼类形态学研究提供了新思路，为建立岱衢族大黄鱼外形指标及开展外形评价提供了有效手段。主成分分析和连续投影算法对形态参数进行特征提取和选择的结果表明，全长/体长、全长/头长、全长/尾柄长、体长/头长、尾柄长/尾柄高这5个形态特征变量是岱衢族大黄鱼F₂、F₃代群体形态差异的主要指标。稀疏表示方法对岱衢族大黄鱼F₂、F₃代选育群体的识别结果表明，岱衢族大黄鱼F₂、F₃代选育群体代系之间的形态多样性较低，难以直接使用较少的形态参数因子来概括群体之间的形态差异；而通过PCA提取主成分得到的特征变量，能够全面反映群体间的综合形态差异，对岱衢族大黄鱼F₂、F₃代选育群体的识别效果最好。