基于特征区域下凸点提取的藻类荧光光谱波长选择方法
张永彬1, 朱丹丹1, 陈颖1,*, 刘喆1, 段玮靓1, 李少华2
1.燕山大学电气工程学院河北省测试计量技术及仪器重点实验室, 河北 秦皇岛 066004
2.河北先河环保科技股份有限公司, 河北 石家庄 050000
*通讯作者 e-mail: chenying@ysu.edu.cn

作者简介: 张永彬, 1992年生, 燕山大学电气工程学院硕士研究生 e-mail: 2824272678@qq.com

摘要

藻华现象的频繁发生严重影响了海洋环境和人类的生产活动, 因此对水体浮游植物的监测十分重要。 三维荧光光谱被广泛应用于水体浮游植物中藻类的群落组成分析和浓度定量分析, 然而三维荧光光谱数据中的信息冗余给藻类定性定量分析带来了一定的影响。 针对光谱信息冗余问题, 提出了特征区域积分与凸点提取相结合的三维荧光光谱波长选择方法。 以抑食金球藻、 细长聚球藻、 小球藻为研究对象, 采用Savitzky-Golay卷积平滑法对三维荧光光谱进行预处理, 解决了因外界因素造成的光谱噪声问题, 采用马氏距离法剔除三维荧光光谱数据集中的异常光谱样本, 运用浓度残差法剔除三维荧光光谱数据集中的异常浓度值样本, 然后通过偏最小二乘回归模型的内部交叉验证均方根误差衡量不同特征区域下凸点的可靠性进行波长变量的选择。 为验证波长筛选方法的有效性, 对三种藻类建立偏最小二乘回归模型, 以内部交叉验证决定系数( R2)、 内部交叉验证均方根误差(RMSECV)作为模型评价指标。 与全光谱数据建立的回归模型进行了比较, 抑食金球藻、 小球藻、 细长聚球藻的波长变量由全谱的1071个分别减少到77个、 75个、 67个, R2分别提高了0.016 4, 0.002和0.032 4, RMSECV分别降低了1.8×105, 2.0×105, 2.6×105。 与UVE方法相比, 抑食金球藻、 小球藻、 细长聚球藻的波长变量分别减少了599个、 357个、 317个, R2分别提高了0.014 5, 0.000 4和0.012 3, RMSECV分别降低了1.6×105, 7.0×104和1.6×105。 经过该方法进行波长变量选择后, 减少了冗余信息, 提高了模型预测能力。

关键词: 浮游植物; 三维荧光光谱; 特征区域; 凸点; 波长选择
中图分类号:O433.4 文献标志码:A
Wavelength Selection Method of Algal Fluorescence Spectrum Based on Convex Point Extraction From Feature Region
ZHANG Yong-bin1, ZHU Dan-dan1, CHEN Ying1,*, LIU Zhe1, DUAN Wei-liang1, LI Shao-hua2
1. Hebei Province Key Laboratory of Test/Measurement Technology and Instrument, School of Electrical Engineering, Yanshan University, Qinhuangdao 066004, China
2. Hebei Sailhero Environmental Protection Hi-tech Co., Ltd., Shijiazhuang 050000, China
*Corresponding author
Abstract

The frequent occurrence of algal bloom seriously affects the Marine environment and human production activities, so it is very important to monitor the phytoplankton in water.3D fluorescence spectroscopy has been widely used in the analysis of algae community composition and the quantitative analysis of algae concentration in water phytoplankton. However, the information redundancy in 3D fluorescence spectrum data has significantly impacted the qualitative and quantitative analysis of algae.In order to solve the problem of spectral information redundancy, a new wavelength selection method of 3D fluorescence spectrum based on the combination of feature region and convex point extraction is proposed.Taking Aureococcus anophagefferens, Chlorella Vulgaris, and Synechococcus elongatus as the research object, the Savitzky-Golay convolution smoothing method was used to preprocess the 3D fluorescence spectrum to solve the problem of spectral noise caused by external factors. The Mahalanobis distance method was used to eliminate the abnormal spectral samples in the 3D fluorescence spectrum data set.The residual concentration method was used to eliminate the abnormal concentration value samples in the 3D fluorescence spectrum data set.Then the reliability of the convex points under different characteristic regions was measured by the root mean square error of cross-validation (RMSECV) of the PLS regression model, and the wavelength variable was selected. In order to verify the effectiveness of the wavelength selection method, the PLS regression model was established for the three algae species, and the determination coefficient ( R2) and root mean square error of cross-validation (RMSECV) were used as the evaluation indexes of the model. Compared with the regression model established with the full spectrum data, the wavelength variables of Aureococcus anophagefferens, Chlorella Vulgaris, and Synechococcus elongatus respectively decreased from 1 071 to 77, 75 and 67, and R2 respectively increased by 0.016 4, 0.002 and 0.032 4. RMSECV was respectively reduced by 1.8×105, 2.0×105 and 2.6×105. Compared with the UVE method, the wavelength variables of Aureococcus anophagefferens, Chlorella Vulgaris, and Synechococcus elongatus were respectively reduced by 599, 357 and 317, and R2 was respectively increased by 0.014 5, 0.000 4 and 0.012 3, RMSECV was respectively decreased by 1.6×105, 7.0×104 and 1.6×105. After the selection of wavelength variables by the method of feature region combined with convex point extraction, the redundant information is reduced, and the model's prediction ability is improved.

Keyword: Phytoplankton; 3D fluorescence spectroscopy; Feature region; Convex point extraction; Wavelength selection
引言

近年来我国工农业发展迅速, 海水富营养化加剧, 水体中浮游藻类大量繁殖, 频繁引发藻华现象, 褐潮成因藻抑食金球藻、 常见藻华中的小球藻和细长聚球藻已经成为海洋生态环境的主要污染, 影响人类健康安全和环境发展, 所以藻类种类和数量的测量对生态环境检测具有重要意义[1, 2]。 浮游植物中含有的色素具有荧光效应, 因此荧光光谱分析是一种有效的测量方法, 其中三维荧光光谱法在藻类测量中应用广泛, 如何选择光谱波长变量从而提高三维荧光光谱的信噪比是目前一个重要研究方向。

三维荧光光谱是在已设置的激发波长和发射波长范围下扫描所有波长点而得到的数据矩阵, 其数据较多, 但并不是光谱区域下的所有波长变量都含有有用信息, 只有少数光谱区域含有物质特征。 若冗余区域波长变量参与建模会影响预测精度, 增加模型的计算, 光谱预处理虽然可以提高信噪比, 但并没有去除不相关的波长变量, 不能最大程度提高模型的精度, 没有减小模型复杂度。

传统的波长选择方法是峰值法, 对固定波长点下的荧光峰进行分析, 但该方法所选择的波长变量过少, 原光谱的数据信息没有得到充分利用, 无法有效反映整体区域荧光强度变化情况[3]。 目前的波长选择方法有无信息变量消除法(uninformative variables elimination, UVE)[4]、 连续投影算法(successive projections algorithm, SPA)[5, 6]、 遗传算法(genetic algorithm, GA)[7, 8]等。 SPA通过投影方式选取线性关系最小的波长组合, 减少光谱信息中的波长变量, 但SPA所选择的波长变量的个数不能超过所用样品数量, 其受样品数量的限制较大。 UVE将随机噪声矩阵加入光谱矩阵并建立偏最小二乘(partial least-square, PLS)回归模型, 通过回归系数的稳定性进行波长变量选择。 GA是对初始群体通过自然遗传机制进行选择、 交换、 突变等算子操作产生新群体的过程, 过程中适应度值较优的变量被保留。 由于UVE的人工随机噪声矩阵和GA的随机初始种群的随机性, 导致了每次波长变量选择的结果不同, 且经过UVE波长选择后剩余的波长变量依然较多[9]

荧光区域积分结合凸点提取的波长选择方法是根据整体区域荧光强度变化选择出能够表征藻类色素的荧光特性的波长变量, 同时避免选择过程的随机性。 本文以抑食金球藻、 小球藻、 细长聚球藻为研究对象, 在对3种藻的原始光谱进行预处理后, 结合荧光区域积分和二元凸函数判定方法, 进行特征区域波长变量筛选, 并将筛选后的光谱数据作为输入建立PLS回归模型, 验证其有效性。

1 基于凸点提取的波长选择原理

特征区域下的凸点提取是通过在主要荧光区域中提取曲面峰值点及峰值周围的点来达到去除荧光光谱冗余信息的目的。 其整体过程分为两步, 首先通过荧光区域积分[10]找到荧光特征区域, 通过阈值改变荧光特征区域的大小; 然后在特征区域内投票统计所有波长点进行凸点提取后的得分, 设定阈值衡量得分可靠性。 调节特征区域和凸点投票统计后得分的阈值, 大于阈值的波长变量被认为利于模型建立, 予以保留。 其实现过程包括:

Step1: 识别光谱数据中的凸点。 三维荧光光谱可以看作是一个关于发射波长(x轴)、 激发波长(y轴)、 荧光强度(z轴)的二元函数z=f(x, y), 采用Savitzky-Golay多项式拟合微分法求偏导, 计算x的偏导数时固定y轴, 计算y的偏导数时固定x轴。 根据二元凸函数[11]判定定理, 三维光谱的任一点(x, y)为凸点需要满足以下条件

2f(x, y)x2> 0, 2f(x, y)y2> 0(1)

2f(x, y)x22f(x, y)y2-2f(x, y)xy20(2)

对所有样本三维荧光光谱的每个波长点的凸点进行统计, 将得到的统计值与样本总数的比值作为个波长点凸点得分。

Step2: 将整个荧光区域划分为若干个单位区域, 由于三维荧光光谱数据是离散数据, 其单位区域的体积积分为

Φi=EXEMI(λEX, λEM)ΔλEXΔλEM(3)

式(3)中, I(λ EX, λ EM)为在激发波长λ EX、 发射波长λ EM处的荧光强度, Δ λ EX为激发波长间隔(取5 nm), Δ λ EM为发射波长间隔(取5 nm)。

设最大单位区域体积积分的 1ηqq=10, 9, 8, 7, 6, 5, 4, 3)作为阈值, 大于该阈值的所有单位区域体积积分的并集作为特征区域。

Step3: 求光谱数据的凸点得分与特征区域的交集, 得到特征区域下的光谱数据凸点得分结果。

Step4: 设定特征区域下的光谱数据凸点得分阈值η t(0.3< η t< 1), 小于设定阈值, 则认为该凸点是不利于建模的波长变量或者是由于噪声引起的波长变量。

Step5: 将选择的波长变量建立偏最小二乘回归模型, 采用留一交叉验证法计算内部交叉验证均方根误差(RMSECV)。

Step6: 遍历所有η q值, 重复Step2— Step5过程, 当RMSECV值最小时所选择的波长变量为最终的波长变量选择结果

2 实验部分

所用的抑食金球藻和小球藻由秦皇岛海洋环境监测中心站提供, 细长聚球藻购自中国科学院淡水藻种库。 采用FS920三维荧光光谱仪对三藻种进行光谱测量, 测量条件: 激发波长范围为400~650 nm, 发射波长范围为630~730 nm, 激发步长和发射步长均为5 nm, 狭缝宽度为 5 nm, 设置发射波长滞后激发波长10 nm, 用1 cm石英比色池每天取样一次进行光谱测量。

实验中所选藻种的培养周期均为14 d, 为采集到更多丰度下的光谱信息, 每种藻培养3份, 每份保持不同起始丰度。 经过多周期测量, 最终获得藻种光谱数据共306个。 不同色素有不同的荧光特性, 藻类中色素的种类决定荧光峰的位置, 色素的含量决定荧光峰的强弱, 图1给出了三种藻类的等高线谱。 图1(a)中440/685, 450/685和470/685 nm三处的荧光峰分别由抑食金球藻的叶绿素a(Chla)、 叶绿素c(Chlc)和叶绿素b(Chlb)造成; 图1(b)中440/685, 470/685, 480/685 nm三处的荧光峰分别由小球藻的Chla、 Chlb和类胡萝卜素造成, 420/685 nm处的荧光峰为Chla的次吸收带引起的次荧光峰; 图1(c)中420/680和440/680 nm分别为细长聚球藻中Chla的次荧光锋和主荧光峰, 620/655和620/680 nm处的荧光峰均由藻蓝蛋白(PC)引起, 但PC对光的吸收和释放特性主要主要表现在620/655 nm处的强荧光峰。 Chla存在于三种藻中, 但在细长聚球藻中Chla的荧光发射波长为680 nm, 可能是由于细长聚球藻中的Chla与蛋白质结合, 改变了荧光的响应特性, 从而引起发射光波长的偏移。

图1 三藻种等高线光谱图
(a): 抑食金球藻; (b): 小球藻; (c): 细长聚球藻
Fig.1 Contour spectra of 3 phytoplankton
(a): Aureococcus anophagefferens; (b): Chlorella vulgaris; (c): Synechococcus elongatus

3 结果与讨论
3.1 光谱数据去噪

由于一些因素的影响, 通过光谱仪获得的荧光光谱既有有用信号, 同时也有噪声信号。 本文通过Savitzky-Golay卷积平滑法去除噪声信号, 平滑窗口宽度为5, 多项式次数为2。 图2和图3为三种藻类样本激发光谱平滑前后对比图, 可看出各藻种光谱的大体形状没有发生变化, 消除了光谱上的毛刺, 将误差信息重新分配到整个光谱中, 使光谱更加平滑, 降低了噪声干扰。

图2 三藻种平滑去噪前的激发光谱图
(a): 抑食金球藻; (b): 小球藻; (c): 细长聚球藻
Fig.2 Excitation spectra of 3 phytoplankton before smooth denoising
(a): Aureococcus anophagefferens; (b): Chlorella vulgaris; (c): Synechococcus elongatus

图3 三藻种平滑去噪后的激发光谱图
(a): 抑食金球藻; (b): 小球藻; (c): 细长聚球藻
Fig.3 Excitation spectra of 3 phytoplankton after smooth denoising
(a): Aureococcus anophagefferens; (b): Chlorella vulgaris; (c): Synechococcus elongatus

3.2 异常样本数据剔除

在光谱采集过程中, 由于测量人员、 测量环境、 测量仪器等因素的影响, 导致所测量的样本中存在一些光谱异常样本和化学值异常样本, 这些异常数据会影响模型的精度, 需要在建模之前剔除[12]

光谱异常样本通过马氏距离法剔除。 计算出各个样本到平均光谱的马氏距离, 通过平均值和标准差来设定阈值, 调节阈值权重系数来改变阈值范围, 剔除不利于建模的样本。 光谱异常样本剔除结果如图4所示, 本文中抑食金球藻、 小球藻、 细长聚球藻的阈值权重系数分别取0.5, 0和0.5, 剔除样本数分别为17个、 19个、 13个。

图4 三藻种马氏距离分布图
(a): 抑食金球藻; (b): 小球藻; (c): 细长聚球藻
Fig.4 Mahalanobis distance distribution of 3 phytoplankton
(a): Aureococcus anophagefferens; (b): Chlorella vulgaris; (c): Synechococcus elongatus

剩余的样本中可能存在化学值异常样本, 本实验中化学值指藻种浓度, 可用浓度残差法剔除。 通过留一交叉验证法依次获得数据集各样本的浓度残差组成浓度残差向量, 设置浓度残差的F统计性检验的显著性水平阈值, 浓度残差大于该阈值的样本为化学值异常样本。 浓度残差剔除异常样本的结果如图5, 抑食金球藻、 小球藻、 细长聚球藻的显著性水平均为0.5, 剔除的浓度异常样本分别为32个、 34个、 34个, 最终可用于建模的样本数分别为53个、 55个、 49个。

图5 三藻种浓度残差分布图
(a): 抑食金球藻; (b): 小球藻; (c): 细长聚球藻
Fig.5 Concentration residual distribution of 3 phytoplankton
(a): Aureococcus anophagefferens; (b): Chlorella vulgaris; (c): Synechococcus elongatus

3.3 基于荧光区域积分和凸点提取的波长选择方法

3.3.1 全波长凸点提取

对各藻类样本数据集下的全谱波长点依次进行凸点统计, 如果满足凸点条件, 则将该波长点保留并将凸点数加1, 以每个波长点统计的凸点数与样本总数的比值作为该点的凸

点统计得分, 图6为各藻类凸点得分统计图, 抑食金球藻、 小球藻、 细长聚球藻的全谱波长点均为1 071个, 被判定为凸点的波长点个数分别为369个、 422个和501个。

图6 三藻种凸点统计图
(a): 抑食金球藻; (b): 小球藻; (c): 细长聚球藻
Fig.6 Convex pointsstatistics of 3 phytoplankton
(a): Aureococcus anophagefferens; (b): Chlorella vulgaris; (c): Synechococcus elongatus

3.3.2 初始特征区域的确定

计算各藻类全谱范围下所有单位荧光区域的体积积分, 初始特征区域阈值取最大单位区域体积积分的1/10, 大于该阈值的单位体积积分的并集作为初始特征区域, 图7为各藻类提取到的初始候选特征区域结果图, 抑食金球藻、 小球藻、 细长聚球藻在初始特征区域下的波长点个数分别为340个、 409个、 525个, 由图6和图7可知虽然主要的色素峰值点既被判定为凸点, 又包含在初始特征区域内, 但在初始候选特征区域外, 有些波长点既不是色素的特征峰又无其他明显的光谱特征信息, 却多次被判定为凸点。

图7 初始候选特征区域
(a): 抑食金球藻; (b): 小球藻; (c): 细长聚球藻
Fig.7 Initial candidate feature region
(a): Aureococcus anophagefferens; (b): Chlorella vulgaris; (c): Synechococcus elongatus

3.3.3 特征区域和凸点提取阈值的确定

图8为初始特征区域下的凸点得分, 抑食金球藻、 小球藻、 细长聚球藻在该区域下被判定为凸点的波长个数分别为131个、 142个、 259个, 已有效去除不含特征信息的凸点, 但此时由噪声引起的光谱凸点波长变量包含在内, 同时由于特征区域阈值选取问题, 初始特征区域并不精确, 区域内可能存在与藻种信息相关性不大的波长点, 依旧有光谱冗余。 图9为阈值nqnt取不同值时RMSECV的变化情况, RMSECV值越小, 模型的预测能力越好, 合适的nqnt值可以提高模型的预测能力。 图9(a)中nq取1/8, 1/7, 1/6和1/5时RMSECV值相等, 图9(b)中nq取1/4和1/3时RMSECV值相等, 是由于在这几个特征区域阈值下, 光谱的积分区域发生了变化, 但波长变量的数量并没有改变, 这时模型的预测能力由凸点的阈值nt决定。 当抑食金球藻、 小球藻、 细长聚球藻的候选特征区域阈值分别为1/10, 1/6和1/8, 凸点得分阈值分别为0.4, 0.5和0.9时RMSECV最佳。 图10为最佳候选特征区域阈值和最佳凸点得分阈值下各藻种筛选波长的最终结果, 抑食金球藻、 小球藻、 细长聚球藻剩余波长点数分别为77, 75和67。

图8 初始候选特征区域下的凸点统计
(a): 抑食金球藻; (b): 小球藻; (c): 细长聚球藻
Fig.8 Statistics of convex points under the initial candidate feature region
(a): Aureococcus anophagefferens; (b): Chlorella vulgaris; (c): Synechococcus elongatus

图9 η qη t取不同值时RMSECV的变化曲线
(a): 抑食金球藻; (b): 小球藻; (c): 细长聚球藻
Fig.9 Change curve of RMSECV with different η q and η t
(a): Aureococcus anophagefferens; (b): Chlorella vulgaris; (c): Synechococcus elongatus

图10 特征区域结合凸点提取波长选择结果
(a): 抑食金球藻; (b): 小球藻; (c): 细长聚球藻
Fig.10 Wavelength selection results of feature region combined with convex point extraction
(a): Aureococcus anophagefferens; (b): Chlorella vulgaris; (c): Synechococcus elongatus

3.4 建模方法及模型评价

分别用全谱波长数据、 特征区域结合凸点提取后波长数据、 UVE提取后波长数据建立PLS回归模型, 采用留一交叉验证法进行内部交叉验证, 以内部交叉验证决定系数(R2)、 内部交叉验证均方根误差(RMSECV)作为评价指标评估波长选择的效果, 其中R2反映模型验证的稳定性, RMSECV反映模型的预测能力, R2越接近1, RMSECV越接近0, 模型的预测性能及稳定性越好。 图11为经过UVE提取的波长变量结果, 表1为各PLS回归模型的R2、 RMSECV值。 抑食金球藻经该方法波长筛选后的R2相对全谱提高了0.016 4, RMSECV降低了1.8× 105, 波长点由全谱的1 071个减少到77个, 经UVE波长筛选后的R2相对全谱提高了0.001 9, RMSECV降低了3.0× 104, 波长点由全谱的1 071个减少到676个; 小球藻经该方法波长筛选后的R2相对全谱提高了0.002, RMSECV降低了2.0× 105, 波长点由全谱的1 071个减少到75个, 经UVE波长筛选后的R2相对全谱提高了0.001 6, RMSECV降低了1.3× 105, 波长点由全谱的1 071个减少到432个; 细长聚球藻经该方法波长筛选后的R2相对全谱提高了0.032 4, RMSECV降低了2.6× 105, 波长点由全谱的1 071个减少到67个, 经UVE波长筛选后的R2相对全谱提高了0.010 1, RMSECV降低了1.0× 105, 波长点由全谱的1 071个减少到384个。 由以上分析可知, 三藻种在经过特征区域结合凸点提取后R2均得到提高, 其中小球藻的R2提高效果并不显著, 但在没有降低模型预测能力的基础上波长变量减小到原来的7%, 该波长提取方法是有效的。 相比UVE, 特征区域结合凸点提取所得到的模型精度更高, 所选择的波长变量更少, 且均分布在含有荧光特征的区域, 保留了所含色素的特征峰。

图11 UVE波长选择结果
(a): 抑食金球藻; (b): 小球藻; (c): 细长聚球藻
Fig.11 Wavelength selection results of UVE
(a): Aureococcus anophagefferens; (b): Chlorella vulgaris; (c): Synechococcus elongatus

表1 三种方法的R2和RMSECV Table 1 R2 and RMSECV of three methods
4 结论

利用实验得到的三种藻类不同浓度下的三维荧光光谱数据, 提出了特征区域结合凸点提取的波长变量选择方法。 根据藻类光谱找到包含荧光峰的初始特征区域, 利用二元凸函数判定定理统计初始特征区域下的凸点, 设定特征区域和凸点统计结果的不同阈值, 通过计算PLS回归模型的RMSECV值衡量特征区域下波长变量凸点提取的可靠性从而选择波长去留。 结果表明, 用该方法得到的波长变量建立的PLS模型优于用全谱波长变量建立的PLS模型, 且选择的波长变量在色素荧光峰附近, 能够很好地表征藻类色素的荧光特性, 为三维荧光光谱分析提供了依据。

参考文献
[1] Nunes S, Latasa M, Delgado M, et al. Deep-Sea Research Part Ⅰ, 2019, 151(5): 103059. [本文引用:1]
[2] Nankabirwa A, DeCrop W, Vand er Meeren T, et al. Ecological Indicators, 2019, 107: 105563. [本文引用:1]
[3] PU Shu-juan, HE Pei-xiang, ZHANG Yue, et al(蒲淑娟, 贺培翔, 张悦, ). Petrochemical Industry Application(石油化工应用), 2021, 40(5): 102. [本文引用:1]
[4] CHEN Yuan-yuan, WANG Zhi-bin, WANG Zhao-ba(陈媛媛, 王志斌, 王召巴). Spectroscopy and Spectral Analysis(光谱学与光谱分析), 2017, 37(1): 299. [本文引用:1]
[5] Zhang Chu, Liu Fei, Kong Wenwen, et al. Sensors, 2015, 15(7): 16576. [本文引用:1]
[6] XIONG Zhi-xin, MA Pu-fan, LIANG Long, et al(熊智新, 马璞璠, 梁龙, ). Transactions of China Pulp and Paper(中国造纸学报), 2019, 34(4): 46. [本文引用:1]
[7] HUA Chen-zhi, ZHAO Ling, SONG Jian-jun(花晨芝, 赵凌, 宋建军). Journal of Sichuan Normal University(四川师范大学学报), 2019, 42(6): 825. [本文引用:1]
[8] MAI Shu-kui, YANG Yang, ZHAO Xiao-bo, et al(买书魁, 杨洋, 赵小波, ). Food Science and Technology(食品科技), 2019, 44(2): 301. [本文引用:1]
[9] Chang Haitao, Zhu Lianqing, Lou Xiaoping, et al, Sensors, 2016, 16(6): 827. [本文引用:1]
[10] CHEN Shi-yu, LI Yan, LI Ai-min(陈诗雨, 李燕, 李爱民). Environmental Science & Technology(环境科学与技术), 2015, 38(5): 64. [本文引用:1]
[11] ZHANG Xue-mao(张学茂)). Journal of Langfang Teachers College(廊坊师范学院学报), 2011, 11(4): 18. [本文引用:1]
[12] SHI Lu-zhen, ZHANG Jing-chuan, WANG Yan-qun, et al(石鲁珍, 张景川, 王彦群, ). Journal of Chinese Agricultural Mechanization(中国农机化学报), 2016, 37(6): 99. [本文引用:1]