基于遥感图像3D–CNN及氮磷循环的水华预测

王 立,吴羽溪,王小艺,刘载文

(北京工商大学人工智能学院智慧环保北京实验室,北京 100048)

1 引言

随着工农业的发展,大量化工生产过程中的氮、磷等营养物质排入江河湖泊,水体富营养化程度不断加剧[1].藻类的大规模繁殖和腐烂降低了水体的溶解氧含量,导致水生动物死亡,严重破坏了水生生态系统[2].同时,藻类水华还会产生多种生物毒素,严重影响饮用水安全,影响人体健康[3].最新发布的2018年《中国生态环境公报》显示,我国城乡结合部大部分城市湖泊和中型湖泊可能面临水体富营养化[4].水污染问题已成为制约我国经济社会发展的重要因素之一,引起了国家和地方政府的重视.因此预测水华暴发的时间和程度和水体富营养化等级是保护水资源的非常重要的手段.

藻类水华的预测方法主要分为两类,基于机理的水质建模方法[5]和基于数据的水质建模方法.对于机理建模方法,模型从单层、单室、单成分、零维的简单模型发展到多层、多室、多成分、三维的复杂模型,并逐渐应用于河流湖泊的污染控制和生态系统的管理.而现有的机理建模方法,仅考虑氮或磷元素等营养物质对藻类生长过程影响,忽略水体底物中的氮磷反馈机制,没有反映氮和磷元素在藻类水华形成过程中的完整循环过程.因此本文考虑水华形成过程中氮循环、磷循环等化学过程,根据“氮–磷–藻”之间的耦合关系及底物反馈机制,对藻类水华形成过程进行生化动力学建模及分析.

对于数据建模方法,可以简单分为线性分析和非线性[6]分析两大类,线性方法目前常用的多元分析方法包括主成分分析[7]、时间序列分析[8–9]等;目前非线性模型主要包括灰色系统[10–11]和神经网络方法[12–14]等.在水华数据来源方面,传统的水质监测需要野外调查和采样分析,受人为、气候和水文条件的限制.由于覆盖范围小,数据量大,难以实现大面积水域的同时测量.遥感技术的发展为水华预测模型的发展提供了新的思路.针对遥感图像数据量大的特点,深度学习在特征提取和模型拟合方面,显示了其独特的潜力和优势[15–18].卷积神经网络(convolutional neural networks,CNN)作为一种经典的深度学习算法,克服了高维数据分析的维数灾难[19].CNN虽然对静态图像识别有很好的效果,但是对于未来时刻图像的动态预测,需要考虑图像像素在历史时刻与当前时刻的时间相关性,建立各时刻相邻帧的像素之间的相关关系,才能合理描述图像随时间的变化规律,预测图像在未来时刻的发展趋势.因此,本文对CNN加入时间维度,建立3D–CNN模型,提取遥感图像的时空分布特征.

综上,本文既研究藻类水华形成过程的生化机理,又研究该过程的数据关系,建立综合预测模型以获得更为科学准确的预测效果,本文的主要贡献有两个:1)对3D–CNN进行优化改进,并用于对遥感图像的特征提取,预测水体富营养化等级,解决海量遥感图像数据处理及时间预测问题;2)提出基于氮磷循环的藻类水华形成过程建模方法,并采用Elman神经网络(ENN)模型融合遥感图像特征,预测水华暴发时间及程度,解决提高水华预测精度的问题.

2 遥感图像特征提取及富营养化预测

2.1 基于3D–CNN的遥感图像特征提取

为了提取连续的叶绿素a浓度遥感图像之间的帧间连接关系,3D–CNN在传统平面卷积神经网络的基础上增加了时间维度,本研究将多个针对叶绿素a浓度的遥感图像的连续帧图像叠加成一个立方体,并使用三维卷积核提取相邻帧之间的叶绿素a浓度的变化信息,从而捕捉到针对叶绿素a浓度的遥感图像在时间和空间维度上的分布特征.在3D–CNN提取叶绿素a浓度的遥感图像特征过程中,本研究在卷积层上进行3D卷积,从上一层获得的叶绿素a浓度特征图上的局部邻域中提取特征,然后应用加性偏压,通过激活函数传递结果.这个过程是通过将一个3D内核卷积到由多个相邻帧叠加在一起形成的立方体来实现的.通过这种构造,将卷积层中提取的叶绿素a浓度的特征映射连接到前一层中的多个相邻帧,从而捕获运动信息.在(x,y,z)位置上第i层第j层的特征图的取值计算方法如式(1)所示:

其中:f(·)代表激活函数,bij是第i层的第j个特征图的偏置项,m是遍历第i-1层特征图上的所有点,Pi是卷积核的高度,Qi是卷积核的宽度,Ri是卷积核的深度,是连接到第i-1层第m个特征图在(p,q,r)的值,p为卷积核高度Pi的求和指标,q为卷积核宽度Qi的求和指标,r为卷积核深度Ri的求和指标,是第i-1层第m个特征图中在(x+p,y+p,z+r)处的值.在卷积过程中,三维卷积核以固定步长对连续图像进行卷积运算,通过选择不同尺寸的卷积核得到各种高阶特征.同时,下采样操作也需要扩展到三维层面,总结和整合统计特征,降低数据维度,提升运算效率.

2.2 基于细菌觅食算法的网络结构优化

3D–CNN中含有大量超参数,这些参数对于神经网络的训练和预测性能有十分重要的影响.细菌觅食算法(bacterial foraging optimization,BFO)是由K.M.Passino于2002年基于Ecoli大肠杆菌在人体肠道内吞噬食物的行为提出的一种新型仿真算法[20].该算法模拟细菌群体的行为,细菌通过趋化、繁殖、驱散(迁移)3种基本操作在P维搜索空间中寻找营养物质即全局最优解.优化问题的解对应搜索空间中细菌的状态,即优化函数适应值.

本文以水体富营养化等级分类性能为适应度函数,利用细菌觅食算法优化卷积核尺寸、卷积核数量及网络层数,从而提升预测模型的分类精度.使用细菌觅食算法优化卷积神经网络结构的流程图如图1所示.

图1 细菌觅食算法优化卷积神经网络流程图Fig.1 Flow chart of convolutional neural network optimized by bacterial foraging algorithm

表1为3D–CNN模型的候选参数,包括滤波器尺寸,滤波器数量和网络层数.本研究的主要目的是最大限度地提高图像预测精度,图像预测精度随所述3D–CNN内部参数的变化而变化,建立函数如式(2):

具体优化步骤的第1步是确定滤波器大小,滤波器数量和层数的候选参数,选定的结果见表1.对于滤波器大小,待优化参数的选择有{2×2,3×3,4×4,5×5,6×6,7×7},深度均为3.根据控制变量法的原则,保持其余参数相同,仅改变滤波器大小,观察预测水体富营养化等级的准确性,如果由于特定的滤波器尺寸而导致精度急剧下降,候选优化参数将从实验中退出,即排除该组参数.滤波器数目的候选值为{4,8,12,16}.在本步骤中,所有训练过程都随着滤波器尺寸和滤波器数量的变化而进行,并记录分类精度,需要注意的是,在第1步中表现较差的候选参数不参与此步骤.网络层数主要考虑{8层、11层、14层}的3个模型,8层网络包含两块卷积层、激活层和池化层的组合、全连接层和分类层,11层与14层结构则比8层结构分别多一块和两块卷积–激活–池化的组合.采用以上各种参数组合对这些模型进行了广泛的训练和测试,观察预测模型的精度变化情况.

表1 卷积神经网络待优化参数表Table 1 Parameters of convolution neural network to be optimized

在BFO算法的优化过程中,存在着大量的初始学习参数用来运行算法.这些输入参数通常是在仿真过程选择的,表2提供了影响该算法收敛特性的输入参数的详细信息.

表2 细菌觅食算法的初始参数Table 2 Initial parameters of bacterial foraging algorithm

2.3 水体富营养化等级预测

基于3D–CNN的遥感图像特征提取,可对水体富营养化等级进行定性预测.水体富营养化是藻类水华形成的重要原因,其程度可以作为水质评价的重要参考.首先根据国际经济合作与发展组织(OECD)规定的水体富营养化分级标准,将湖泊营养程度贫营养(＜3 μg/L)、中营养(3～11 μg/L)、富营养(11～78 μg/L)、重富营养(＞78 μg/L)4种状态,如表3所示.

表3 富营养化分级标准(OECD)Table 3 Eutrophication classification standard(OECD)

针对太湖的实际情况,由于太湖常年处于富营养化状态,叶绿素a浓度的波动范围无法覆盖整个分级标准,而是集中分布在19～25 μg/L之间,若根据OECD的分级标准,则处于贫营养状态和重富营养状态的样本数量为0,CNN网络训练过程中会出现严重的样本分布不均问题,由此会影响到预测模型精度.

为了便于网络的训练,本研究提出根据太湖实际叶绿素a浓度范围,将OECD中富营养状态等级由轻至重进一步细化为4个等级,分别是富营养I级(11～20 μg/L)、富营养II级(21～22 μg/L)、富营养III级(23～24 μg/L)、富营养IV级(25～78 μg/L)共4种状态,见表4.

表4 富营养状态细化分级标准(太湖)Table 4 Detailed classification standard of eutrophication(Taihu Lake)

基于3D–CNN对水体富营养化等级预测过程的结构图如图2所示,通过对t-2,t-1及t时刻的叶绿素遥感图像进行连续的卷积和池化操作,预测未来时刻的水体富营养化等级k(t+i).

图2 基于3D–CNN的水体富营养化等级预测Fig.2 Prediction of water eutrophication level based on 3D–CNN

3 水华形成过程分析及水华暴发预测

3.1 氮–磷–藻的生化动力学建模及分析

随着工业化进程的加快和人类活动的影响,化工生产中的过量氮和磷等营养物质进入湖泊,是形成藻类水华的最根本原因.

以太湖的藻类水华为例,其水华形成过程主要考虑的化学过程包括氮循环、磷循环以及溶解氧的平衡等,如图3所示.

图3 水华形成的生化过程Fig.3 Biochemical process of water bloom formation

天然水体中的藻类通过光合作用对氮、磷元素吸收的化学表达式为

其中(CH2O)106(NH3)16(H3PO4)为藻类的化学分子式.

同时,藻类又通过呼吸作用把氮、磷元素释放于水中的化学表达式为

因此,根据“氮–磷–藻”3者的耦合关系及底物反馈机制,建立藻类水华形成过程的三维生化动力学时变参数模型如下所示:

其中:ca为t时刻叶绿素a的浓度,N,P分别为t时刻总氮、总磷的浓度;G(t)为时变的藻类生长率,

T(t)为水温随时间变化的t时刻的数值,gmax为藻类最大生长率;N0,P0分别为总氮、总磷在初始时刻的浓度;T1为藻类的环境损失率;T2为氮的耗损率,gN为藻类对氮的吸收率,dN为藻类代谢分解对氮的贡献率;T3为磷的耗损率,gP为藻类对磷的吸收率,dP为藻类代谢分解对磷的贡献率,g(N),g(P)为总氮、总磷增长函数的Monod方程,

其中KN和KP分别代表藻类生长对氮和磷的半饱和常数.

针对模型中3个微分方程的众多常值参数,本文采用四阶龙格库塔数值算法与遗传算法相结合,实现对模型中常值参数的智能优化率定[20],流程如图4所示.

图4 常值参数优化率定流程图Fig.4 Flow chart of constant parameter optimization calibration

由于“氮–磷–藻”三维生化动力学系统属于高维非线性动力学系统,其动力学行为较复杂,故首先采用分岔理论确定系统的分岔集,求出各个范围的分岔点后再采用中心流形方法对系统降维,具体研究临界局部稳定性和分岔类型,确定水华暴发的临界条件.

水华暴发不同程度的临界条件如表5所示.

表5 水华暴发程度及其临界条件Table 5 Outbreak degree and critical conditions of water bloom

其中,G0,G1分别为三维生化动力学系统第1个平衡点和第2个平衡点的雅克比矩阵特征值为0时对应的G(t)值.

根据对模型的非线性动力学分析,当时变参数G(t)在(0,G0)范围内,系统稳定,无水华产生;当G(t)处于(G0,G1)时,系统振荡稳定,此时当叶绿素a峰值小于阈值时,产生轻度水华,当叶绿素a峰值大于阈值时,产生中度水华;当G(t)大于G1时,系统不稳定,产生重度水华.

3.2 基于ENN融合的水华暴发预测

“氮–磷–藻”三维生化动力学时变参数模型能够解释藻类水华形成过程的生化机理,但无法反映水域的整体分布情况,而遥感图像能够反映整体水域变化的时空分布信息,因此,本文综合基于机理建模与数据建模方法,将生化动力学模型与遥感图像3D–CNN特征结合,并考虑水质监测站点的叶绿素a浓度历史数据,采用ENN融合3者信息,预测水华暴发的时间和程度.其中ENN预测模型的输入数据共由3部分组成:由3D–CNN提取的遥感图片特征、叶绿素a浓度计算值以及历史时刻站点采样值.该预测模型在传统的反向传播模型基础上新增一个承接层作为一步延时算子,使系统具有适应时变特性的能力,增强网络的全局稳定性.另外由于Elman神经网络是动态反馈型网络,能够内部反馈、存储和利用过去时刻输出信息,实现动态系统的映射并直接反应系统的动态特性.本研究通过训练Elman神经网络输出叶绿素a浓度的预测值,基于“氮–磷–藻”三维生化动力学时变参数模型,利用分岔理论和中心流形方法确定水华爆发的时间及程度.

藻类水华预测模型的结构图如图5所示.其中Im(t-2,t-1,t)代表采样站点叶绿素a浓度的历史时刻测量值,Ics(t)代表当前时刻叶绿素a浓度的机理模型计算值,以上二者与3D–CNN提取的特征共同输入至ENN模型,得到叶绿素a浓度的精确预测值.

图5 预测模型结构图Fig.5 Structure diagram of prediction model

此外,本文对预测模型引入反馈校正机制,以减小模型失配、干扰等因素对控制系统的影响,利用模型预测输出值与实际输出值的差值来修正模型输出的预测值,误差的计算方法如式(6)所示:

其中:Vtrue表示叶绿素a浓度的真值,Vpred表示叶绿素a浓度的预测值,i表示预测的时间长度,采用i维的加权反馈向量h=[h1h2··· hi]对预测模型进行修正,修正方法如式(7)所示:

预测模型的反馈机制示意图如图6所示.

图6 反馈机制示意图Fig.6 Schematic diagram of feedback mechanism

综上所述,本研究总体流程图如图7所示.

图7 研究思路流程图Fig.7 Flow chart of research ideas

4 实验验证

4.1 数据集

太湖叶绿素a浓度遥感图像数据的来源是在晴天(云覆盖率70%以下)的情况下,针对太湖地区MODIS卫星250米分辨率遥感图像,进行叶绿素a浓度反演,获取太湖叶绿素a浓度的时空分布情况.本研究选取2009–2010年间的152幅RGB图像,每3幅序列图像作为3D–CNN的一组输入,共获得150组数据,其中前115组数据作为训练样本,后35组作为测试样本.图像分辨率为273×276像素.

4.2 遥感图像特征提取及富营养化预测结果

利用细菌觅食算法优化后的CNN结构及具体参数如表6所示.优化后的CNN结构包含两次卷积和池化过程,由全连接层将特征进行整合,输出富营养化预测等级,本文同时给出3D–CNN和2D–CNN的结果作为对比.

表6 2D–CNN和3D–CNN的优化结果Table 6 Optimization results of 2D–CNN and 3D–CNN

为了清晰起见,本文展示了典型的遥感图像和各富营养化水平的叶绿素a浓度反演图.2009年8月26日的遥感影像对应第I级富营养化状态,7月1日、8月17日、7月25日的影像分别代表II,III,IV级.遥感图像清晰地捕捉到太湖流域的轮廓,从图8可以明显看出,随着遥感图像中绿色部分的增加,水体富营养化程度也随之增加.在图8的第3列中,灰度的明暗表示叶绿素a浓度的具体反演情况,显示藻类水华的空间分布状态.

图8 各富营养化等级对应的典型遥感图片Fig.8 Typical remote sensing images corresponding to each eutrophication level

2D–CNN网络将3个时刻的叶绿素反演遥感图像垂直拼接作为一组输入,通过训练网络预测下一时刻的水体富营养化等级,本质上实现的仍是平面卷积.而3D–CNN则将3张遥感反演图按照时间顺序排列,将平面卷积拓展为立体卷积操作,提取包含遥感图像的帧间联系的时空特征.高效的特征提取是对图像进行分析和进一步处理的前提与基础,因此,3D–CNN能够实现更加精准而全面的预测.图9和图10分别展示了基于2D–CNN和3D–CNN建模过程的特征示意图,灰度图的明暗程度代表叶绿素浓度的具体取值,即藻类水华的严重程度,经过卷积和池化操作后特征图的数量由卷积核数量决定.从两幅图可以看出,浅层网络用于提取灰度图的纹理、细节特征;深层网络则用于提取灰度图的轮廓、形状特征.随着卷积和池化过程的交替进行,图片的像素逐渐降低,提取的特征越来越具有代表性.

图9 基于2D–CNN提取的太湖特征图Fig.9 Feature maps of Taihu Lake based on 2D–CNN

图10 基于3D–CNN提取的太湖特征图Fig.10 Feature maps of Taihu Lake based on 3D–CNN

对于水体富营养化等级的预测结果,将上述2D–CNN和3D–CNN分别用于太湖的遥感反演叶绿素a浓度图,并从10次运行结果中选取最佳值.

在35组试验数据中,属于I类水体的数据有9组,II,III,IV类数据分别有16,6,4组.表7列举出了针对各富营养化等级的预测正确样本数及准确率,并最终给出了2D–CNN和3D–CNN的最佳总体准确率,分别为57.14%和71.43%.

表7 水体富营养化等级预测结果Table 7 Prediction results of water eutrophication level

表8详细对比了2D–CNN与3D–CNN方法的具体性能.本研究训练集样本数为115组,测试集样本数为35组,分别针对训练集与测试集从预测准确率与训练时间两方面进行评估.结果表明,3D–CNN方法获得了更加出色的预测准确率,在10次运行结果中,使用该方法的最佳预测精度和平均预测精度均优于使用2D–CNN得到的预测结果.

表8 2D–CNN与3D–CNN的详细性能比较Table 8 Detailed performance comparison between 2D–CNN and 3D–CNN

4.3 水华形成过程分析及水华暴发预测结果

根据水华形成过程的“氮–磷–藻”3维生化动力学时变参数模型,采用本文方法对其常值参数率定结果如表9所示.

表9 常值参数率定结果Table 9 Calibration results of constant parameters

水温为模型的时变参数,采样数据如图11所示.

对于另一个代理变量APP的下载种类数量(新媒体使用广度)来说，它对于新生代女性农民工工作匹配的影响也存在一个适当的范围。APP的下载种类数量也必须控制在一定的区间范围之内才能对新生代女性农民工的工作匹配产生正向的影响。考虑到新生代女性农民工的群体特点，本文认为这里还有其特殊的原因。在实地调研访谈的过程中可以发现，新生代女性农民工对征婚交友类网站 APP 的使用普遍增多。这些征婚交友类网站APP的出现在一定程度上满足了新生代女性农民工对于婚恋交友的需求，但是这种APP下载种类数量的增多对其工作匹配并不会产生积极的影响作用。具体结果如表5所示：

代入常值参数和时变参数后,计算得到总氮和总磷的时间历程,以及“氮–磷–藻”三维生化动力学系统的相轨迹如图12–13所示.

图12 总氮、总磷的时间历程Fig.12 Time history of TN and TP

图13 “氮–磷–藻”三维生化动力学系统相轨迹Fig.13 “N-P-alga”three dimensional biochemical kinetic system phase trajectory

由于藻类生长率G(t)是根据水温等影响因素随时间变化的,导致“氮–磷–藻”三维生化动力学系统的第2个平衡点也随时间变化.通过对式(3)求解平衡点及稳定性分析可知,在第31个采样点处,系统处于震荡稳定状态,此时需再通过判断叶绿素a浓度的预测值是否超过阈值确定水华暴发的程度是轻度还是中度.

对于叶绿素a浓度的预测结果,图14展示出了实测值分别与遥感图像3D–CNN特征预测值,以及基于“氮–磷–藻”三维生化动力学时变参数模型、遥感图像3D–CNN特征和叶绿素a浓度水质监测历史数据的ENN融合预测值的比较曲线.

图14 叶绿素a浓度预测曲线Fig.14 Prediction curve of chlorophyll a concentration

本研究针对富营养化等级与叶绿素浓度实行不同的评价标准,对预测模型的性能进行评定.对于富营养化等级,精度μ的计算公式如式(8)所示:

对于叶绿素浓度的预测情况,本研究采用标准误差指标衡量,计算方式如式(9)所示：

其中:Vtrue代表叶绿素浓度的真值,Vpred表示叶绿素浓度的预测值.

计算表明,3D–CNN预测值和ENN融合预测值的精度分别达到88.04%和93.97%,RMSE分别为1.0810和0.3864.实践证明,ENN融合预测相比3D–CNN预测的精度更高.

同时,根据图14可知,叶绿素a浓度在第31个采样点处达到峰值,与三维生化动力学系统分析结果相吻合,而该峰值并未达到叶绿素a浓度的阈值20 μg/L,因此可知水华暴发的程度是轻度水华.

5 结论

本文以太湖流域的遥感图像及“氮–磷–藻”生化动力学系统为研究对象,分析藻类水华形成过程,提出基于3D–CNN的遥感图像特征提取及“氮–磷–藻”三维生化动力学时变模型的水体富营养化等级及水华暴发预测方法.通过太湖流域监测数据的实例验证,得到以下结论.

1) 针对遥感图像特征提取,本文采用的3D–CNN方法可以同时提取太湖流域整体叶绿素a浓度的时间和空间分布特征,克服了传统方法缺乏对流域整体时空分布信息分析,以及难以对遥感图像海量数据分析的问题,为水体富营养化等级及水华暴发预测提高了信息的完备性.

2) 针对水华形成过程分析,本文考虑了藻类水华形成过程中重要的氮循环和磷循环的化学过程,所建立的“氮–磷–藻”三维生化动力学时变参数模型能够体现三者的耦合关系和反馈机制,同时提出相应的模型参数智能优化率定及非线性动力学分析方法,为水华暴发预测提供了过程的可解释性.

3) 针对水华暴发预测,本文将遥感图像的3D–CNN特征、“氮–磷–藻”三维生化动力学时变参数模型与水质监测历史数据三者的信息通过ENN融合,使得水华预测模型同时具备数据驱动和机理驱动模型的信息完整且机理准确优势,从而提高了对水华暴发时间和程度预测的精度.