余心杰，吳雄飛，沈偉良1

（1.浙江大學寧波理工學院，浙江寧波315100；2.寧波市海洋與漁業(yè)研究院，浙江寧波315010）

形態(tài)學特征參數是大黃魚（Pseudosciaena crocea Richardson）良種選育、種質鑒別及特征檢測的重要依據[1-2]。我國水產科技工作者已經圍繞形態(tài)學特征參數對大黃魚不同地理群體、不同育種群體、不同家系、雌雄性別、種質特征及鑒別比較等開展了深入研究[3-7]，有效推動了我國大黃魚人工繁殖與苗種培育進程。然而，目前在各類大黃魚形態(tài)學研究過程中，一方面，形態(tài)參數測定依靠傳統(tǒng)手工尺子測量，難以實現(xiàn)對大樣本量研究對象的快速采集，且容易損壞活魚樣本[8]；另一方面，所采用的魚類形態(tài)學分析方法普遍限于線性回歸分析、判別分析、聚類分析等傳統(tǒng)多元統(tǒng)計分析方法[9]，在魚類形態(tài)特征挖掘及群體識別模型構建方法方面，至今尚缺乏有效的手段。因此，研究建立一套高效的魚類形態(tài)參數測量、形態(tài)特征挖掘和群體識別方法，對于大黃魚良種選育、種質鑒別等具有重要意義。

本文以岱衢族大黃魚F2、F3代選育群體為研究對象，通過計算機視覺（computer vision,CV）[10]測定2類選育群體的形態(tài)參數，使用主成分分析（principal component analysis,PCA）[11]和連續(xù)投影算法(successive projections algorithm,SPA)[12]對選育群體的形態(tài)參數進行分析，并獲得能表達不同選育群體之間形態(tài)差異的特征變量。在此基礎上，以特征變量為輸入，結合稀疏表示（sparse representation,SR）[13]方法建立岱衢族大黃魚選育群體分類識別模型，實現(xiàn)不同選育群體的自動識別。

1 材料與方法

1.1 實驗材料

2015年11月、2016年11月在浙江省寧波市象山港灣水產苗種有限公司育苗場（象山縣黃避岙鄉(xiāng)高泥村鳳凰礁），分批采集普通網箱（5 m×10 m）養(yǎng)殖、投喂冰鮮魚餌料的15月齡岱衢族大黃魚F2、F3代選育群體樣本各600尾，用于大黃魚形態(tài)參數的計算機視覺測量、形態(tài)特征挖掘和選育群體識別建模研究。

1.2 實驗方法

1.2.1 形態(tài)參數的計算機視覺測量

1.2.1.1 形態(tài)參數測量系統(tǒng)

采用自己研制的魚類形態(tài)參數計算機視覺測量系統(tǒng)，對大黃魚形態(tài)參數和質量進行測量。該系統(tǒng)硬件結構和實物如圖1所示，主要包括框架主體、LED光源（YX-BL25844，山東省煙臺永信視覺技術有限公司）、框架主體上的透明玻璃平臺、設置在透明玻璃平臺四角下端的稱重傳感器（KNX/3kg，江蘇省蘇州坤宏電子有限公司）、放置在透明玻璃平臺下端用于拍攝魚體圖像的工業(yè)相機（SJM300C，臺灣視覺美深圳有限公司）和鏡頭（VM06012MP，廣州智賽電子有限公司）。工業(yè)相機與計算機信號連接，采用下置方式采集透明玻璃平臺上被測樣本魚體圖像，鏡頭與被測樣本距離保持固定，從而避免不同樣本魚體厚度差異所帶來的尺寸測量誤差。稱重傳感器通過稱重控制器（XK3190-A27E，上海耀華稱重系統(tǒng)有限公司）與計算機信號連接，同步傳輸魚體質量數據。

圖1 系統(tǒng)硬件結構（A）和實物（B）圖Fig.1 System hardwarestructure(A)and entity(B)images

用安裝在計算機中的魚類形態(tài)參數測量軟件（圖2）采集被測樣本的數字圖像和重量數據，并利用數字圖像對魚體形態(tài)參數進行精確測量分析。該軟件是在Microsoft Visual Studio 2012開發(fā)環(huán)境下采用C#語言并基于開放視覺圖形處理庫Emgu CV編寫，其主要功能包括：1）數據采集管理。用于采集和管理被測樣本的數字圖像、質量和形態(tài)參數數據。2）參數設置。用于設定工業(yè)相機的圖像分辨率、曝光時間、白平衡等參數；利用棋盤格標定板進行相機標定（消除圖像畸變）和像素當量標定（即圖像中一個像素點代表的實際物理尺寸）。3）形態(tài)參數測量。用于測量魚體尺寸、外形輪廓曲率等參數。

圖2 魚類形態(tài)參數測量軟件Fig.2 Software for measuring fish morphological parameters

本研究中魚類形態(tài)參數計算機視覺測量系統(tǒng)的設計與開發(fā)細節(jié)、測量精度驗證等相關內容參考余心杰等[14]的報道。

1.2.1.2 形態(tài)參數采集與處理

利用上述魚類形態(tài)學參數的計算機視覺測量系統(tǒng)對大黃魚樣本的形態(tài)參數進行快速活體測量。先將活魚樣本分批次（每組20尾）放入質量濃度為15 mg/L的丁香酚水溶液中，麻醉1～2 min至魚體側翻不運動，逐個取出樣本放到透明玻璃平臺上，通過魚類形態(tài)參數測量軟件，采集魚體數字圖像和質量數據，并在計算機中保存。完成全部樣本數字圖像和質量數據采集之后，在軟件中采用手工鼠標標記方式，對魚體形態(tài)參數進行批量測定。

1.2.2 形態(tài)特征挖掘

在數據挖掘領域，特征提?。╢eature extraction,FE）和特征選擇(feature selection,FS)是從原始特征中找出本質特征的有效手段，其中特征提取是指由原始特征經過線性或非線性變換得到較少數量但具有表達能力的新特征，特征選擇是指直接從原始特征中選取最相關的特征子集。目前，在魚類形態(tài)學研究中，主要通過對樣本魚體形態(tài)性狀參數之間的相關性關系和通徑分析，尋找出形態(tài)特征（顯著相關的形態(tài)性狀參數或參數比例）作為生長預測、種群劃分等的依據，本質上還是屬于特征選擇范疇。為了探索建立魚類形態(tài)特征挖掘方法，本文分別從基于主成分分析（PCA）的特征提取和基于連續(xù)投影算法（SPA）的特征選擇2個角度，開展岱衢族大黃魚不同選育群體的形態(tài)特征挖掘研究。

1.2.2.1 基于PCA的形態(tài)特征提取

PCA是將多個相關的特征變量轉化為少數幾個相互獨立的新特征變量的有效分析方法。假設X為采集到的一個n×m的大黃魚樣本數據矩陣，其中的每一列對應一個特征變量，每一行對應一個樣本，則X∈Rn×m。先將X標準化為：

其中：M=[m1m2… mm]，為X的均值；[s1s2，…，sm]為X的標準差。

X′的PCA模型可表示為：

式中：t1，t2，…，tm∈Rn，為得分向量，也叫主元；p1，p2，…，pm∈Rm，為載荷向量。如果可以用前k（k＜m，一般取累計貢獻率大于85%對應的主元個數）個主元來代替數據中的主要變化，那么X′的PCA模型可表示為：

其中E為誤差矩陣。X′可以近似表示為：

本文采用主成分分析對大黃魚樣本數據矩陣X進行變換，用變換后的前k（k＜24）個主元代替原始24個特征變量信息，實現(xiàn)魚體形態(tài)特征提取。

1.2.2.2 基于SPA的形態(tài)特征選擇

SPA的目標是選取彼此之間共性較少的變量組合，且盡可能地保留原始數據信息，從而達到特征選擇的目的。SPA的工作原理是通過迭代的方法，從一個變量開始，每次循環(huán)都計算其在未選入變量上的投影，將投影向量最大的變量引入到變量組合中，每個新選入的變量，都與前一個的線性關系最小。對于由樣本數n和特征變量數m組成的大黃魚樣本數據矩陣X∈Rn×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?{K(0),…,K(p－1)}}；

4）計算剩余列向量x′j(j∈s)與當前所選向量x′K(p－1)的投影：

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

6）令 x′j=Px′j，j∈s；p=p+1。如果 p＜k，返回到第3）步循環(huán)計算。

對于每一個初始K(0)，循環(huán)一次后利用多元線性回歸（multiple linear regression,MLR）進行交互驗證分析，最小均方根誤差（root mean square error,RMSE）對應的K(p)即為特征選擇結果。本文采用SPA對大黃魚樣本數據矩陣X中的24個特征變量進行選擇，最終選取k個特征變量組合{K(p)，p=0，…，k－1}來表達大黃魚樣本形態(tài)差異。

1.2.3 稀疏表示方法

對大黃魚不同選育群體的形態(tài)特征提取和選擇之后，不同選育群體識別的基本問題就是基于形態(tài)特征數據，根據L類訓練樣本，確定測試樣本屬于哪一類。目前常用的模式識別方法是基于反向傳播（back propagation,BP）神經網絡或支持向量機（support vector machine,SVM）等建立相應的分類模型，將測試樣本的特征數據輸入該模型，得到相應的識別結果[15]。與上述方法需通過多個參數優(yōu)化選擇過程來建立一個分類模型不同，稀疏表示法直接通過測試樣本在每個訓練樣本上的投影系數進行識別分類[16]，識別過程簡單，易實現(xiàn)。

假設大黃魚選育群體種類有L類；每一類的訓練樣本數分別為N1，N2，…，Ni，…，NL；每個訓練樣本的特征變量為m維，用列向量f來表示，則f∈Rm；同一類訓練樣本數據位于一個線性子空間中，那么新測試樣本數據可看成是同類樣本數據的線性組合。

將充足的訓練樣本給第i個類，然后將第i類中的Ni個訓練樣本數據組成矩陣Ai（稱之為訓練樣本矩陣）的列，具體表示為：

任意來自于相同類的測試樣本y∈Rm，近似滿足于以下由訓練樣本組成的線性組合：

隨后，為所有類的訓練樣本集定義一個新的訓練樣本矩陣A作為所有L類訓練樣本的集合：

則y的線性表示可以根據所有L類訓練樣本寫成另一種形式：

式 中 ：列 向 量 x=[0,0,…,ai,1,ai,2,…,ai,j,…,ai,Ni,0,0,…,0]T是一個系數向量；元素ai,j表示y在A中第i類的第j個訓練樣本上的投影系數。

若y來自第i類樣本，則x在其所屬類別的Ai上的投影系數不為0，而在其他訓練樣本類別的Al(l≠i)上的投影系數為0。也就是說，y的解是稀疏的。

稀疏表示和壓縮感知的理論研究表明，如果解x是稀疏的，則可通過解決如下的L-1最小化范數問題來求解式（10）：

式中：x?為x的近似解；‖‖·1表示向量中各元素的絕對值之和。

在理想情況下，x?的非零系數中最大的系數值如果對應第i個類，且較大的系數基本都同第i個類有關，就可將y歸到那個類中。但是，在實際應用中噪聲等問題可能會導致x?中有一部分與其他類有關的小的非零系數項出現(xiàn)。因此，對于每個類i，可定義一個函數δi，用來在稀疏解x?中選取與第i類有關的系數。δi(x?)是一個新的向量，其中除了與第i類有關的系數外，令其他的系數值都為零，重新生成給定的測試樣本y的一個近似表示，將其記為：

顯然，y?i與y距離越小，則 y?i屬于第i類的可能性越高。因此，可通過計算y?i與y的最小殘差值，即求解以下方程得到殘差最小的i對應的類別，即為y的類別標志。

1.2.4 基于形態(tài)特征和稀疏表示的大黃魚選育群體識別方法

通過上述方法，建立完整的基于形態(tài)特征和稀疏表示的岱衢族大黃魚F2、F3代選育群體識別方法的步驟如下。

1）形態(tài)特征挖掘：分別采用PCA和SPA方法，對所采集的岱衢族大黃魚F2、F3代選育群體的24個形態(tài)參數進行特征提取和選擇，得到k個特征變量。

2）輸入：在每一種選育群體中隨機選取300個樣本作為訓練樣本，剩余的300個樣本作為測試樣本。利用所有訓練樣本的k個特征變量，構成如下訓練樣本矩陣：

式中：AF2、AF3表示岱衢族大黃魚F2、F3代選育群體的訓練樣本矩陣；fi,j為第i種選育群體的第j個訓練樣本的k個特征變量組成的向量列。矩陣A為F2、F3代選育群體的全體訓練樣本矩陣，其大小為k×600。選取某個選育群體中的測試樣本y。

3）L-1范數最小化求解：利用式(11)，得到系數向量x?。這里的計算誤差閾值設為0.01，迭代次數為150次。

4）計算殘差值：利用式(13)，計算 F2、F3代選育群體2個類的y?i與y的最小殘差值。

5）輸出：具有最小殘差的類作為判定測試樣本y的類別。

上述PCA算法、SPA算法、稀疏表示算法在Matlab 7.9下編程實現(xiàn)，其中稀疏表示算法中的L-1范數最小化方法采用美國斯坦福大學研究人員編寫的Matlab軟件包[17]實現(xiàn)；計算機硬件環(huán)境為Intel(R)Core(TM)i5-2300 CPU、主頻2.80 GHz、內存4 096 MB的PC機，軟件環(huán)境為Windows7.0操作系統(tǒng)。

2 結果與討論

2.1 形態(tài)測量結果

利用魚體形態(tài)參數的計算機視覺測量系統(tǒng)采集大黃魚活魚數字圖像和體質量的速度約150尾/h，滿足大樣本量魚類快速、活體采集要求。實驗共采集岱衢族大黃魚F2、F3代選育群體的數字圖像和體質量樣本1 200個（F2、F3代各600尾），測定了體質量、體長、體寬、頭部曲率、全長、全高、吻長、尾柄長、尾柄高、頭長等10個參數。2個群體樣本都為15月齡，樣本的體質量、全長等規(guī)格相近（表1）。計算相關形態(tài)參數比值得到的24個形態(tài)參數描述見表2。

表1 岱衢族大黃魚F2、F3代選育群體樣本規(guī)格Table 1 Population sample size of Daiqu large yellow croaker F2 and F3

表2 形態(tài)參數描述Table2 Description of morphological parameters

2類選育群體樣本的外形如圖3所示。根據樣本外形觀測，F(xiàn)3代選育群體與F2代選育群體外形具有一定差異，這是由于在岱衢族大黃魚選育過程中進行了外形淘汰選育，F(xiàn)3代群體樣本的總體體型偏長，尤其是尾柄形態(tài)較細長。

2.2 形態(tài)特征挖掘

2.2.1 主成分分析特征提取

對岱衢族大黃魚F2、F3代選育群體樣本的24個形態(tài)參數進行主成分分析，得到前10個主成分（PC1～PC10）的累計方差貢獻率（表3）和2個群體的24個形態(tài)參數的主成分載荷量分布（圖4）。

由表3和圖4可知：第1主成分PC1貢獻率為58.38%，對其載荷量較大的是全長/尾柄長（TL/CPL）、體長/尾柄長（BL/CPL）、尾柄長/尾柄高（CPL/CPH）、體高/體質量（BH/BM）；第2主成分PC2貢獻率為21.87%，對其載荷量較大的是全長/頭長（TL/HL）、體長/頭長（BL/HL）；第3主成分PC3貢獻率為8.94%，對其載荷量較大的是全長/體長（TL/BL）、體長/吻長（BL/KL）。前3個主成分的累計貢獻率達到89.19%，說明PCA提取得到的這3個主成分因子已經能夠較好地表達2類選育群體的形態(tài)參數信息?；赑CA的載荷分析表明，上述對第1、2、3主成分取值影響大的8個形態(tài)參數，可以被認為是選育群體樣本之間形態(tài)差異的重要因子。根據2類選育群體樣本的第1、2、3主成分繪制散點圖(圖5)，從中可以看出，雖然F2、F3代選育群體的部分樣本有重疊交錯，但是這2類群體樣本整體上還是可以區(qū)分的。

圖3 岱衢族大黃魚F2代選育個體(A)和F3代選育個體(B)外形Fig.3 Body shape of F2(A)and F3(B)generations of Daiqu largeyellow croaker

表3 前10個主成分累計方差貢獻率Table 3 Cumulative variance contribution rate of the first 10 principal components

圖4 主成分載荷量Fig.4 Principal component loads

圖5 F2、F3代選育群體的第1、2、3主成分散點分布Fig.5 Distribution of the first,second,and third principal dispersal pointsof F2 and F3 generations

總之，由于岱衢族大黃魚F2、F3代選育群體形態(tài)參數具有差異，可以利用PCA提取的主成分或者利用PCA載荷分析得到的8個特征形態(tài)參數，對不同選育群體樣本進行識別。

2.2.2 連續(xù)投影算法特征選擇

圖6 SPA選擇的不同變量數的最小均方根誤差（RMSE）分布Fig.6 Minimum root mean squareerror(RMSE)distributions of different variablesselected by SPA

利用SPA從岱衢族大黃魚F2、F3代選育群體樣本的24個形態(tài)參數中選擇共線性最小的有效特征形態(tài)參數。經SPA選擇的不同變量個數的最小均方根誤差（RMSE）分布圖如圖6所示。從中可以看出，提取出7個特征形態(tài)參數時RMSE達到低點，此后RMSE趨于平穩(wěn)，因此，可以從原始24個形態(tài)參數中選擇這7個特征形態(tài)參數，它們分別是全長/體長（TL/BL）、全長/頭長（TL/HL）、全長/尾柄長（TL/CPL）、體長/頭長（BL/HL）、尾柄長/尾柄高（CPL/CPH）、全長/尾柄高（TL/CPH）、體長/尾柄高（BL/CPH）。

通過比較PCA載荷分析得到的8個特征形態(tài)參數和SPA選擇的7個特征形態(tài)參數之后發(fā)現(xiàn)，2種特征挖掘方法都篩選出了全長/體長（TL/BL）、全長/頭長（TL/HL）、全長/尾柄長（TL/CPL）、體長/頭長（BL/HL）、尾柄長/尾柄高（CPL/CPH）這5個特征形態(tài)參數。從形態(tài)參數比值含義來看，5個特征形態(tài)參數主要反映了大黃魚頭部、尾鰭、尾柄部位長度在魚體長中的占比，集中表現(xiàn)在全長、體長、頭長、尾柄長4個指標，這與直接觀察所得到的F3代群體較F2代群體樣本總體體型偏長、尾柄形態(tài)細長的形態(tài)差異的結論一致。

2.3 岱衢族大黃魚F2、F3代選育群體識別

分別利用PCA提取的前10個主成分（PCA變換主元特征）、PCA載荷分析得到的8個特征形態(tài)參數組合（PCA選擇特征）及SPA選擇的7個特征形態(tài)參數組合（SPA選擇特征），按照1.2.4節(jié)稀疏表示識別步驟，對岱衢族大黃魚F2代（DQF2）、F3代（DQF3）選育群體樣本進行識別。不同特征對選育群體的識別結果如表4所示?？梢钥闯觯?組特征對岱衢族大黃魚F2代（DQF2）、F3代（DQF3）選育群體樣本進行識別的平均準確率都在75%以上，且岱衢族大黃魚F2代、F3代群體之間存在較顯著的形態(tài)差異，能夠進行形態(tài)區(qū)分。其中：利用PCA提取的前10個主成分對2類選育群體樣本進行識別，得到最高的平均識別準確率為88.3%。這是由于主成分分析能夠在丟失最少信息的條件下，將原有眾多變量濃縮成較少的能夠反映原有變量絕大部分信息的綜合指標，因此能夠全面反映群體間的綜合形態(tài)差異，得到較好的識別效果。利用PCA載荷分析選擇的8個特征形態(tài)參數組合和利用SPA選擇的7個特征形態(tài)參數組合對2類選育群體樣本的識別準確率分別為79.0%、80.5%，兩者識別效果沒有明顯差距，但是相比PCA提取的前10個主成分，其識別效果不夠理想。究其原因，可能是岱衢族大黃魚F2、F3代選育群體代系之間雖然存在形態(tài)差異，但是形態(tài)多樣性較低，難以直接使用較少的形態(tài)參數因子來概括群體之間的形態(tài)差異。因此，在岱衢族大黃魚選育群體代系之間的形態(tài)差異描述中要充分利用數據挖掘手段，采用特征提取方法獲得綜合性形態(tài)特征指標來對其形態(tài)差異進行描述，從而克服在傳統(tǒng)魚類形態(tài)學研究中傳統(tǒng)特征選擇方法所獲得的單一性形態(tài)特征指標不全面等缺點。

表4 不同特征對岱衢族大黃魚F2代（DQF2）、F3代（DQF3）選育群體樣本的識別結果Table4 Identification resultsof F2 generation(DQF2)and F3 generation(DQF3)sampleswith different characteristics

3 結論

本文采用計算機視覺技術測定大黃魚形態(tài)參數，通過主成分分析和連續(xù)投影算法提取形態(tài)特征變量，并結合稀疏表示方法實現(xiàn)了對岱衢族大黃魚F2、F3代選育群體的識別建模，為魚類形態(tài)學研究提供了新思路，為建立岱衢族大黃魚外形指標及開展外形評價提供了有效手段。主成分分析和連續(xù)投影算法對形態(tài)參數進行特征提取和選擇的結果表明，全長/體長、全長/頭長、全長/尾柄長、體長/頭長、尾柄長/尾柄高這5個形態(tài)特征變量是岱衢族大黃魚F2、F3代群體形態(tài)差異的主要指標。稀疏表示方法對岱衢族大黃魚F2、F3代選育群體的識別結果表明，岱衢族大黃魚F2、F3代選育群體代系之間的形態(tài)多樣性較低，難以直接使用較少的形態(tài)參數因子來概括群體之間的形態(tài)差異；而通過PCA提取主成分得到的特征變量，能夠全面反映群體間的綜合形態(tài)差異，對岱衢族大黃魚F2、F3代選育群體的識別效果最好。