青鳉鱼的行为特征提取研究

刘翠棉，饶凯锋，李婧，唐亮，裴琨，谷金峰，刘勇，王伟，姜杰，马梅，王子健

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.