青鳉魚的行為特征提取研究

劉翠棉，饒凱鋒，李婧，唐亮，裴琨，谷金峰，劉勇，王偉，姜杰，馬梅，王子健

1. 石家莊市環境監控中心，石家莊 050000 2. 中國科學院生態環境研究中心，環境模擬與污染控制國家重點聯合實驗室，北京 100085 3. 中國科學院生態環境研究中心，中國科學院飲用水科學與技術重點實驗室，北京 100085 4. 石家莊市環境綜合執法支隊，石家莊 050000 5. 無錫中科水質環境技術有限公司，無錫 214024 6. 中國科學院大學資源與環境學院，北京 101407

在水環境的質量監測和安全判斷中，主要有2種技術手段，一種是基于常規指標的定量分析[1-3]，一種是基于水生生物的定性分析[4-8]。其中，基于常規指標的檢測技術起步較早，研究較為廣泛，市場中的成熟產品較多，且有國家標準支撐[9]。而基于水生生物的檢測手段由于缺乏國家標準，生物個體差異明顯，具有極強的不可預測性，導致其起步較晚，不過由于簡便、快捷、直觀和無二次污染等優點，在未來的環境監測領域該技術可以應用到常規指標檢測的前端，作為常規指標檢測的啟動條件。國內常見的在線生物監測設備主要基于發光菌[4]、水溞[5]、藻類[6]和魚類[7-8]等技術，其中，魚類是水生態系統中相對于發光菌、水溞和藻類等更高級的生物，與人類對環境污染物的反應更加接近，因此，可作為水環境監測領域更理想的受試生物。

當外界環境發生改變時，魚類首先會通過行為調節機制快速適應環境的變化，這種行為改變的強弱與環境脅迫的程度有很大的關系[10]，行為改變比病理損傷或死亡發生的時間更早，準確快速識別這一改變可以在理論和應用上為水環境的在線監測提供有力的支持。為了能夠有效觀察到該變化，需采用不破壞魚類正常生活環境的非接觸式監測手段[7-8]，在自然狀態下觀測受試魚類的生理特征和運動特性。本研究中通過低壓高頻的生物傳感器采集青鳉魚在不同類型、不同濃度特征污染物脅迫下的行為電信號。

青鳉魚的行為電信號是一類非平穩非線性時間序列[11]，快速傅里葉變換(FFT)方法不能凸顯行為信號的異常變化[12]；小波變換(Wavelet)方法[13]雖然可以同時從時域和頻域上對行為信號進行解析，但小波基的選擇往往會因生物個體的差異而不具備自適應性，對于實時精確時頻分析比較困難；經驗模態分解(EMD)方法[14]能夠得到本征模態函數分量，但行為電信號被分解后各個分量不具有可解釋性。借鑒常規指標的異常檢測技術[15]，假設青鳉魚的行為電信號在某個時間間隔內是符合高斯分布的，結合圖像處理的思想，采用直方圖的方法對某個時間間隔內的行為電信號進行直方圖統計，將該直方圖與高斯模型進行比較，使得信號從高維降到低維，獲取該時間間隔內的青鳉魚的本質特征，利用該特征進行后續的相關分析。

本研究在觀察高濃度特征污染物暴露實驗下青鳉魚行為變化的基礎上，提出基于直方圖統計方法的降維算法，用低維數據作為本質特征替代固定時間間隔內的行為數據，為后續的異常行為識別提供基礎數據。

1 材料與方法(Materials and methods)

1.1 數據來源

使用的青鳉魚是符合毒理學實驗要求且采用流水繁殖的標準模式魚，利用雙層生物行為傳感器(圖1)，以20次·s-1的采樣頻率連續獲取包含青鳉魚行為信息的電信號。傳感器的上層為放置青鳉魚的暴露層，下層為空白對照層。傳感器的電極采用316L不銹鋼鍍鉻材料制成，與腔體平行，且每組電極由相對設置的2對電極組成，在腔體內形成一個低壓高頻電場。每50毫秒通過串口傳輸一次采集數據給計算機，實驗中所用工控機的主要配置為凌動D525處理器和2 G的內存、XP操作系統。實驗中使用的特征污染物為2,4,6-三氯酚。

在線監測設備要求系統具有較低的誤報率和漏報率，為了減少誤報的干擾，實驗檢測的污染物濃度≥1 TU。所以，在本文的在線系統中，采集軟件開啟穩定運行4 h±5 min后進行特征污染物暴露實驗，針對2,4,6-三氯酚設計其暴露濃度梯度分別為1、2、5和10 TU。其中，1 TU為該種污染物在48 h流水暴露的情況下，青鳉魚的半致死濃度。

圖1 雙層低壓高頻生物傳感器Fig. 1 Double-layer low-voltage high-frequency biological behavior sensor

1.2 分析方法

1.2.1 行為信號的復雜性

由于生物行為具有不可預測性，青鳉魚在生物傳感器中的行為信號差異明顯，不同個體的青鳉魚的原始行為信號是完全不同的。在2個生物行為傳感器中分別各放置一條青鳉魚，截取同一分鐘內2個傳感器的原始信號，如圖2所示，2個原始信號完全不同，說明個體差異明顯。即使是同一條青鳉魚，在不同時空范圍內運動所產生的原始電信號也完全不同。圖3(a)和(b)分別顯示的是同一條青鳉魚在相鄰的2 min內的原始電信號，圖3(c)和(d)是同一條青鳉魚在2個不同傳感器中獲取的1 min的原始數據，這也證明了生物行為的無規律性和復雜性。

工控條件下環境噪聲或實驗室條件下人為造成的傳感器震動等都會造成青鳉魚行為上的顯著變化。由于青鳉魚被限定在傳感器上層活動，其移動范圍有限，無法有效避開噪聲源，這便會在有限空間內引起行為信號發生異常改變，這一改變會對預警設備產生一定的影響，如果算法不能有效識別并濾除由此產生的異常信號便會引發設備誤報警。截取人為敲動傳感器前后各30 min且經過趨勢算法處理過的行為信號(圖4)，可知，在敲動傳感器時，行為信號有個較為明顯的躍升，與環境脅迫閾值模型[16]極其相似，很難判斷是否為污染導致的行為變化。

低壓高頻傳感器在將監測數據傳輸至采集控制軟件的過程中，要經過模數轉換的采集卡，通常采集卡會包含一些電路噪聲，主要包括：內部的導電微粒不連續地造成的低頻噪聲、半導體PN結兩端勢壘區電壓的變化引起累積在此區域的電荷數量改變而產生的散粒噪聲、長期使用過程中導電體內部電子的無規則運動產生的高頻熱噪聲等。如果這些噪聲產生的信號掩蓋了青鳉魚產生的電信號，將很難通過算法來識別出真正的行為信號。由于實際設備中加持在傳感器上的是交流電，工頻交流電也會給采集到的信號附加上50 Hz的噪聲信號。

另外，青鳉魚的正常生物鐘現象也會給預警算法帶來一定的困難，進入和退出生物鐘時的行為信號分別與重金屬和有機特征污染物的暴露特性相對應。圖5顯示了7 d空白對照實驗中的行為信號，橫坐標為時間尺度(單位：min)，縱坐標為青鳉魚的行為信號強度。由圖5可知，青鳉魚每天基本在固定的時間點會進入和退出生物鐘，進入生物鐘時行為信號下降到一定的程度，處于生物鐘的過程中，青鳉魚也不是完全處于睡眠狀態，有時行為信號也會有一定的波動。而退出生物鐘時，行為信號會恢復到進入生物鐘之前的強度，但也有例外的情況，恢復后的行為信號也許強于或弱于前一天的信號。長時間運行后行為信號總體上會是一個慢慢下降的過程，如果在規定的時間內不進行運維，整個信號會慢慢降低到魚死亡時對應的水平。

圖2 2條青鳉魚在2個生物傳感器中的原始信號Fig. 2 Original signals from two different biosensors with two different medakas

圖3 同一條青鳉魚在同一生物傳感器中相鄰2 min內的原始信號((a)和(b))；同一青鳉魚在不同生物傳感器中的原始信號((c)和(d))Fig. 3 The original signal of the same medaka fish within 2 minutes of the same biosensor ((a) and (b)); the original signal of the same medaka fish in different biosensors ((c) and (d))

圖4 噪聲對青鳉魚的行為影響Fig. 4 Effect of noise on the behavior of medaka

上面描述的是已知情況下的行為信號特性，實際使用過程中可能還存在一些未被發現的情況。而且僅從已知的情況來看，行為信號的曲線也是很復雜的，在進行預警算法研發的過程中，需要找到一個或幾個特征能夠表征青鳉魚的真實狀態屬性，使用這些特征進行算法分析。

1.2.2 常用信號處理方法

在信號處理領域，EMD特征分析、FFT法和小波變換等方法都具有很強的信號分解能力，但在實際應用過程中，也都存在各自的局限性。例如，EMD法能夠將行為信號的趨勢提取出來，但在噪聲信號比較強的情況下，分解出來的趨勢不具有解釋性，而且識別的效果也不是很理想；同樣FFT法雖然能夠得到行為信號在頻域上的一些特性，但對實時系統而言，信號的趨勢及一些細節特征未能很好地提取出來，處理的效果略差于EMD法；在小波變換的方法中，需要進一步分析才能對青鳉魚的行為信號處理有很好的效果，但青鳉魚的個體差異導致了每次更換標準模式魚后小波的基函數選擇比較困難。圖6(a)、(b)和(c)分別顯示了EMD法、FFT法和小波變換法這3種方法的行為分析結果，X軸為時間間隔，Y軸為行為強度。

圖5 青鳉魚生物鐘的行為信號Fig. 5 The behavioral signals of biological clock of medaka

由圖6可知，1 TU的2,4,6-三氯酚作用于青鳉魚后，通過分析EMD法的第6個分量可以在15 min左右進行有效預警，但該分量不具備穩定性，更換一批青鳉魚或是更換2,4,6-三氯酚的濃度，第6個分量都不會有這一相同的結果，或許有的情況下EMD算法的結果非常有效，但也有不管使用哪個分量都會失效的情況。圖6(c)中采用的是db1小波，數據處理完后，數據量減少1/2，從圖中無法判斷對什么位置進行預警效果最好，只能再通過諸如序貫貝葉斯等算法進一步分析才能比較有效地給出答案，不過小波變換的方法在相同參數下，結果比EMD法要穩定一些。另外，FFT方法處理后，即便再使用其他算法進行深度分析，也很難給出預警結果。對于復雜的生物個體而言這些經典的信號處理算法往往不能直接使用，所以，需要對原始的行為信號進行特征提取，得到表征生物個體特性的信號，再使用FFT法、EMD法和小波變換法等會取得好的效果。

1.2.3 問題的產生與分析

任何一種信號處理方法，都需要找到一個能夠表征青鳉魚的行為特征的數據才能進行有效的分析。顯然，原始的青鳉魚的電信號數據是不能直接用于處理的，如何消除個體差異、噪聲等的影響是值得思考的問題。筆者所在團隊的前期實驗研究表明，不管是低濃度或高濃度特征污染物，還是單一或混合污染物，青鳉魚在暴露實驗開始后，污染物都會產生短暫的興奮作用，而達到一個臨界點后產生抑制作用，期間還會出現死亡或者掙扎死亡的情況。因此，在實驗的過程中，筆者觀察并記錄了特征污染物作用下，不同濃度、不同種類的污染物對青鳉魚的影響。從空白的行為信號中，隨機挑取1 min的數據進行畫圖，用同樣的方法獲取中毒點附近約10 min中的1 min數據畫圖，觀察圖形發現了一個顯著的問題：按照前文所述的實驗條件，1 min的數據包含1 200個數據點(采集頻率為20次·s-1)，暴露前和暴露后圖形具有極大的差別，如圖7所示，圖7(a)是暴露前某1分鐘的正常行為數據，圖7(b)是暴露后觀察到青鳉魚的狀態發生明顯改變的1 min數據，可以看出在振幅和交替的頻率方面2個是完全不同的。

由圖7可知，只要能夠找到一種辦法用極低的維數表達出這個差別，便可以掌握青鳉魚的行為信號規律。直觀地看，能夠表征這個差別的只有振幅和交替的頻率(產生振幅變化的時間間隔比未暴露前短)，用2個變量來描述振幅和交替頻率的變化較為困難，通常的做法是將1 200個點進行降維，若讓降維模型不產生過擬合或欠擬合，需要的樣本點數須是維數的15倍左右[17]，即實時緩沖區需要存儲18 000 min(15 min×1 200)的行為信號，按一天24 h計算，需要12.5 d后算法才能訓練好模型，并給出第一個算法結果，顯然這個過程在實際應用中是不可接受的。

圖6 常用信號處理算法的結果注：FFT法表示快速傅里葉變換法，該法處理后強度都集中在低頻部分；小波變換后信號復雜度沒有發生變化；EMD法表示經驗模態分解法，該法處理后信號趨勢已經變得明顯；14:34表示2,4,6-三氯酚的加藥時間點，14:49表示相對于閾值的預警點。Fig. 6 Results of common signal processing algorithmsNote: The FFT processing result is concentrated in the low frequency part; the signal does not change after wavelet transform; the signal trend of EMD processing has become obvious; 14:34 is administration time for 2,4,6-trichlorophenol; 14:49 is time for early-warning based on threshold.

圖7 暴露實驗前(a)、后(b)青鳉魚的行為特征Fig. 7 Behavior feature of medaka before (a) and after (b) exposure experiment

1.2.4 基于直方圖的特征提取

圖像處理領域中比較常用的一種處理方法是直方圖，直方圖可以檢驗數據分布的類型，分析數據是否服從正態分布，判斷數據有無異常，同時還可以直觀地判斷分布中心是否偏離理論值，判斷數據分布范圍是否滿足一定的要求。工業應用的很多場景數據都是近似符合高斯模型的，如圖8所示，在暴露實驗前后高斯模型圖的中心位置、高度和面積等都會有所不同。

圖8 高斯模型Fig. 8 Gaussian model

在青鳉魚行為特征提取的研究中，假設利用低壓高頻傳感器采集到的青鳉魚的行為信號在1 min內也是近似符合高斯分布的。那么將圖7中的2種情況，用直方圖對1 200個點進行統計，由于行為信號強度在[0 1]區間內，直方圖統計過程中以0.01為間隔，畫出的直方圖如圖9所示，圖9(a)是圖7(a)對應的直方圖，而圖9(b)是圖7(b)對應的直方圖。由圖9(a)和圖9(b)可知，正常狀態下的直方圖的高度、非零部分的寬度和面積都不同，且高度與振幅的改變有很大的關系，振幅的波動越劇烈，振幅在[0 1]區間的分布越均勻，導致直方圖的高度就越小。同時，非零寬度和面積也都隨著暴露實驗時間的改變而改變。因此，可以將直方圖所對應的這100維數據作為這1 min內的行為特征，觀察圖9的直方圖還可以對這個行為特征進行優化，使用的時候可以僅用非零寬度所對應的那部分數據作為最終的行為特征，這樣維數會更低，效果會更好。

基于直方圖統計的算法解決了分析過程中遇到的降維問題，可以實時地將青鳉魚的行為信號從1 200維降到100維，且不涉及模型的訓練，避免了因模型訓練而帶來的過擬合和欠擬合現象。該方法的另外一個優點是不用花很長的時間訓練模型，也不需要獲取經驗或未來數據的均值和方差，直接對每分鐘的數據進行處理即可。獲得該100維數據后可以使用無監督的降維方法，比如主成分分析(PCA)法等對其再次降維，保留直方圖中特征值最大的特征向量，用累積貢獻率來截取最終的行為特征，為后續的識別算法提供基礎數據。

2 結果(Results)

為了驗證該特征提取算法的有效性，實驗過程中使用不同梯度濃度的2,4,6-三氯酚作為特征污染物，對青鳉魚進行暴露實驗。詳細實驗記錄列于表1，暴露通常在軟件開啟后4 h±5 min的范圍內進行，整個實驗持續24 h。

2,4,6-三氯酚是環境中主要的有機污染物之一，也是工業生產中的重要原料，它對皮膚及粘膜具有強烈的腐蝕作用，對各種細胞有直接損害。因此，將其作為特征污染物，對青鳉魚進行暴露實驗。圖10(a)是采集到的原始信號，其中，紅色虛線的地方是暴露實驗的開始時間，精確到分鐘，從紅色虛線位置開始對青鳉魚的行為活動觀察了1 h，觀察的結果與圖10(a)的趨勢是吻合的。圖10(b)～圖10(e)分別截取了第100、242、305、1 300分鐘的1 min數據，按照上文闡述的算法，對其進行直方圖統計，第100分鐘時為未開始毒性暴露實驗前的正常行為，可以看出整體比較符合高斯模型，而第242分鐘是剛剛開始毒性暴露實驗的最初1 min數據的直方圖，與第100分鐘的數據相比，略有變化，但不明顯。第305分鐘是筆者認為的最早且最佳的預警點，此刻預警不會產生誤報也不會漏報，由直方圖可知，高斯模型的高度下降很明顯，非零部分的寬度也變得很大，顯著區別于正常和剛經毒性暴露時刻的行為。由圖10(a)可知，第1 300分鐘時魚已經死亡，對應的圖10(e)表現為高度值很大，非零寬度值很小。為了模擬實時在線分析的狀態，把原始信號按1次·min-1的頻率進行統計，在統計過程中由于后續信號的值未知，因此是真實的在線分析，對得到的直方圖使用在線PCA方法進行進一步降維，抽取其中一維數據作圖于圖10(f)。暴露實驗開始前以及未產生顯著變化的行為信號都被壓縮在一個平穩的狀態，而行為變化比較劇烈的位置被凸顯出來，圖中15:37所標記的圓形點是相對于閾值的預警點，從暴露實驗開始至該點約持續了63 min，當然若略降低閾值條件，預警時間還可以再提前一點。綜合來看，基于直方圖統計的方法，比單純FFT法、EMD法和小波變換法等方法性能要優越得多，也可以將提取到的特征應用在EMD法、小波變換法等方法中，通過分析這些特征，EMD法、小波變換法等方法的預警效果也不會差。不過這種特征上的差別，如果利用聚類法、一類支持向量機(OneClass SVM)法等方法效果會更好一些。

2 TU的2,4,6-三氯酚對魚的影響會更大一些，行為曲線的變化也會更明顯，如圖11所示。在暴露實驗開始后一段時間內，青鳉魚對毒性的反應更加明顯，因此，中毒后會比中毒前正常信號的振幅要大，且比1 TU的更加明顯。對于選定點的直方圖也與圖10很相似。由圖11(f)可知，在暴露開始后的第23分鐘，行為變化就已經遠遠超出了正常信號的水平，選擇此刻進行預警，同樣不會產生誤報和漏報的現象。從環境脅迫閾值模型的角度來描述，可以解釋為該時刻已經開始毒性累積，且到達了一定程度，在后續的一段時間魚會不停地掙扎直到最終死亡。

圖9 青鳉魚行為信號的直方圖Fig. 9 Behavioral histogram of medaka

表1 1 TU的2,4,6-三氯酚(TCP)實驗數據Table 1 1 TU of 2,4,6-trichlorophenol (TCP) exposure experiment data

注：1 TU 2,4,6-三氯酚對應的濃度為2.3 mg·L-1。

Note: The corresponding concentration of 2,4,6-trichlorophenol is 2.3 mg·L-1.

圖10 1 TU的2,4,6-三氯酚暴露實驗注：14:34表示2,4,6-三氯酚的加藥時間點，15:37表示相對于閾值的預警點。Fig. 10 1 TU of 2,4,6-trichlorophenol exposure experimentNote: 14:34 is administration time for 2,4,6-trichlorophenol; 15:37 is time for early-warning based on threshold.

濃度≥5 TU的2,4,6-三氯酚對青鳉魚的作用效果更加明顯，從暴露實驗開始到觀察到魚的行為發生改變，僅僅是幾分鐘的時間，圖12(a)也證實了這一點，直方圖與1 TU或2 TU的很類似，在線PCA分析的結果表明，在暴露實驗開始后的第12分鐘(如果濃度更高，時間會更短)魚的行為就已經發生了劇烈的變化，其實對于這種濃度的行為信號，即使不作任何特征提取，直接使用經典的信號處理方法，也幾乎都能準確預警。只是若不做特征處理，暴露實驗開始前一段時間有可能會發生誤報，而暴露實驗開始后預警的時間可能會在30 min左右，利用特征的話，可以將時間大大縮短。

圖11 2 TU的2,4,6-三氯酚暴露實驗注：14:33表示2,4,6-三氯酚的加藥時間點，14:56表示相對于閾值的預警點。Fig. 11 2 TU of 2,4,6-trichlorophenol exposure experimentNote: 14:33 is administration time for 2,4,6-trichlorophenol; 14:56 is time for early-warning based on threshold.

3 討論(Discussion)

低壓高頻的行為傳感器，獲取的是魚在電場加持的環境下由于行為變化而產生的電信號，而信號受到多種因素尤其是生物個體差異的影響導致其屬于非平穩、非線性的范疇，實際設備中采用8個通道、每個通道放3條魚的形式來消除生物個體帶來的信號差異性。如果不對信號進行去噪、特征提取，實際使用過程中將會遇到困難，通常的信號處理算法可以對某一種或某幾種情況適用，而不能解決所有情況下的預警準確性問題。

長期大量的實驗研究表明，對于魚類敏感的有機污染物而言，不管是低濃度或高濃度，還是單一或混合有機污染物，青鳉魚在毒物暴露后其行為都符合環境脅迫閾值模型。主要表現在暴露后，青鳉魚會受到短暫的興奮作用，達到一個行為變化臨界點后改為受到抑制作用，在高濃度甚至一些低濃度污染物暴露實驗中還可能會出現掙扎死亡的情況。而重金屬類的污染物對魚的作用機制不同，因此，基于電信號的技術中，行為曲線是一個逐漸累積至緩慢下降的過程，下降的時間和曲率與污染物濃度也存在劑量-響應關系。青鳉魚的行為變化與水體有機污染物之間存在良好的劑量-響應關系，對青鳉魚的行為電信號處理過程中，首要的任務是如何凸顯異常信號，并減少正常信號的波動特征。通過提取能夠表征青鳉魚本質的行為來準確發現其異常變化，從而達到準確及時判斷是否為水質污染所導致的效果。異常發生后，可以通過模式識別的相關算法對特征進行處理，最終確定是否有污染發生。將該模型和算法整合到生物綜合毒性的在線連續監測設備中，能夠進一步優化突發性污染事故的生物綜合毒性在線連續監測技術和設備。

圖12 5 TU的2,4,6-三氯酚暴露實驗注：14:31表示2,4,6-三氯酚的加藥時間點，14:43表示相對于閾值的預警點。Fig. 12 5 TU of 2,4,6-trichlorophenol exposure experimentNote: 14:31 is administration time for 2,4,6-trichlorophenol; 14:43 is time for early-warning based on threshold.