CN105912790A - 基于深度回归模型的遥感水质监测方法 - Google Patents

Abstract

本发明公开了一种基于深度回归模型的遥感水质监测方法，对待处理遥感图像预处理，得到选定水体的遥感数据；对监测水体区域采样点进行采样，得到采样点的实际水质参数和对应的遥感数据；将所有采样点的实际水质参数和遥感数据作为训练样本，利用多层自动编码机逐层进行训练，得到每层自动编码机的输出，直至最后一层自动编码机训练完毕，得到多层自动编码机的输出；最后将该多层自动编码机的输出利用支持向量机回归算法反演水质参数，得到水质回归模型；利用水质回归模型对待处理遥感图像进行反演，得到待处理遥感图像中水体需要监测项目的浓度值。本方法对各种水体适应性好、监测精度高、适合长期连续监测。

Description

基于深度回归模型的遥感水质监测方法

技术领域

本发明涉及水质监测，具体指一种基于深度回归模型的遥感水质监测方法，属于水质监测分析技术领域。

背景技术

水是自然环境中最重要的影响因子之一，是人类生产、生活不可缺少的重要资源之一，水质的好坏直接关系到人类的生存与发展。水质监测是水质评价与水污染防治的主要依据，随着水体污染问题的日益严重，水质污染问题成为了社会经济可持续发展必须解决的重大问题，尤其是内陆水体，其水质好坏直接影响到国民生产和人们的日常生活。因此，准确、快捷的水质监测显得尤为重要。

常规监测手段通常是在水域内定点定剖面，常年累月进行监测，采样进行物理和化学分析才能得到准确的水质信息，受人力、物力和气候、水文条件的限制，采集的数据量不可能太多，对于整个水域而言，这些监测点数据只具有局部和典型的代表意义，但是对于大多数区域的水质的状况并不能反映出来，且成本高、速度慢、难于长时间跟踪监测，不能给出这些表面水质指标在空间和时间上的分布状况，不能及时的全面反映水域生态环境的变化。随着遥感技术的发展，尤其是遥感行业与光谱分辨率的提高，为应用遥感技术进行水体污染监测和研究开辟了新途径。利用卫星遥感信息进行大面积范围内水质空间分布及动态的定量分析，能够在一定程度上弥补水面采样观测时空间隔大且费时的缺陷，还能发现常规方法难以揭示的污染物排放源、迁移扩散方向、影响范围以及与清洁水混合解释的特征，有利于查明污染物的来龙去脉，为科学地布设地面监测点提供依据，卫星遥感还以其特有的周期性优势成为持久监测区域乃至全球尺度上湖泊，水库的有效手段。利用遥感技术进行水质监测主要是由于被污染水体具有独特的有别于清洁水体的光谱特征，这些光谱特征在于其特定波长的光的吸收或反射，而且这些光谱特征能够为遥感器所捕获并在遥感图像中体现出来。由于遥感图像具有直观明了、宏观性强的特点，能清楚地反映出区域或整个流域污染现状和空间分布特征。利用多时相遥感数据还可以对同一流域水体污染历史和污染趋势做出研究和预测，为水资源保护和规划以及可持续发展提供动态基础数据和科学决策依据。

在水质预测研究初期，主要采用一些简单的时间序列预测方法，包括移动平均法、指数平滑法、马尔可夫预测法和回归预测法等。这些方法在分析未来的结果与过去、现在各种因素之间的关系时，效果比较好，且数据处理简单，但需要时间稳定延续，并不适合长期预测，并且数据需要服从特殊的统计分布，只能按相同规律构建线性模型。

但是，实际应用中水质时间序列呈现非线性模型，主要预测方法有：灰色预测方法、人工神经网络、混沌理论、支持向量机算法等。灰色系统是邓聚龙教授于1982年提出的，在分析一些贫信息系统时具有较好的适应性，主要缺点是数据灰度越高，预测精度越低，不适合长期预测。人工神经网络在高精确的非线性映射中，具有较强的自适应能力，运用BP神经网络模型构建了河流地表水中有机物和其主要影响因素间的映射关系，其预测结果客观，具有较强泛化能力，但BP神经网络在样本不足时，计算精度不高，同时神经网络常因结构过大而容易出现“过拟合”。支持向量机是一种以结构风险最小化准则为基础的机器学习方法，在解决非线性、小样本及高维数问题中展现出很多特有的优势，近年支持向量机在水质预测中，取得了较大的成果，但核函数和参数选取困难，同时模型对不同水质数据特点的适应性差。

2006年，加拿大多伦多大学Hinton教授建立了深度学习模型，它通过模拟建立人脑的分层模型结构，能对原始数据逐级提取出由底层到高层更加抽象的特征。由于深度学习具有优异的特征学习能力，得到的特征对数据有更本质的刻划。近年来，不仅在理论验证方面得到了成功，而且在实际应用方面，也得到了巨大且迅速的推广和发展。在理论模型方面，主要分成了三种类型结构，其中包括生成性深度结构、区分性深度结构和混合型结构。在实际应用方面，主要在语言识别和信号处理、物体识别、自然语言处理和多任务学习方面得到了成功的应用。目前，深度学习的模型主要有：受限玻尔兹曼机(RestrictedBoltzmann Machine)及其深层模型一一深度置信网络(Deep Belief Network)、深度玻尔兹曼机(Deep Boltzmann Machine)；自动编码机(Autoencoder)、降噪自动编码机(Denoising Autoencoder)及其深层模型一一堆栈自动编码机(Stacked Autoencoder)；Pooling Unit以及卷积神经网络(Convolutional Neural Networks)等。

发明内容

针对现有技术存在的上述不足，本发明的目的在于提供一种基于深度回归模型的遥感水质监测方法，本方法对各种水体适应性好、监测精度高、适合长期连续监测。

本发明的技术方案是这样实现的：

基于深度回归模型的遥感水质监测方法，步骤如下，

1)将遥感图像数据导入计算机，提取出包含需要监测水体的矩形范围作为待处理遥感图像；

2)利用掩模矩阵对待处理遥感图像每个波段进行掩膜，并从每个波段中去掉光谱值为0的背景信息，再做归一化处理，得到选定水体的遥感数据；

3)在监测水体区域选取不同的点作为采样点，对采样点进行需要监测项目的水质参数实地采集，得到采样点的实际水质参数，并确定所有采样点对应的遥感数据；采样点的实际水质参数包含监测项浓度较大，中等及较小的数据分布；

4)将所有采样点的实际水质参数和采样点的遥感数据作为训练样本，利用多层自动编码机逐层进行训练，以逐层提取遥感水质图像的特征，从而得到每层自动编码机的输出，并作为下一层自动编码机的输入；如此递推，直至最后一层自动编码机训练完毕，得到多层自动编码机的输出；最后将该多层自动编码机的输出利用支持向量机回归算法反演水质参数，得到水质回归模型；

5)利用所构建的水质回归模型对待处理遥感图像进行反演，得到待处理遥感图像中水体需要监测项目的浓度值；根据其处于不同的浓度区间赋予不同的颜色，从而直观显示整个水体水质参数的浓度情况，实现水体水质遥感监测。

第4)步多层自动编码机中每层自动编码机按如下方法训练，

假设单层自动编码机输入层x包含d个输入神经元，隐层y含有h个神经元；重构层z含有与输入层相同的神经元数目；网络中各个神经元的激活函数设定为sigmoid函数：

$f\left(x\right)=sigmoid\left(x\right)=\frac{1}{1+e^{-x}}---\left(1\right)$

在训练过程中，模型首先将输入向量x(x∈R^d)映射到隐层中，以得到隐层特征y(y∈R^h)，然后将隐层特征y(y∈R^h)映射到重构层，得到重构层特征z(z∈R^d)；这两个映射环节分别用如下两个公式表示：

y＝sigmoid(W_yx+b_y) (2)

z＝sigmoid(W_zy+b_z) (3)

其中W_y，W_z分别是输入层-隐层以及隐层-重构层的连接权值矩阵；b_y，b_z是隐层与重构层每个神经元的偏置；为了减少模型中需要训练的参数，给出如下约束：

W_y＝W_z＝W (4)

对于这样的自动编码机，训练就意味着学习模型中包含的三组参数W，b_y，b_z；

模型的训练则是一个对该目标函数的最优化过程，即调节W,b_y,b_z来最小化重构的输出与输入之间的差异；故确定这样的目标函数：

$\begin{array}{cc}idealreconstruction=\arg & \underset{W,b_y,b_z}{\min }\[\cos t(x,z)\]\end{array}---\left(5\right)$

这里的z在x给定的情况下取决于模型的参数W,b_y,b_z；cost(x,z)是重建与输入的差异的度量，定义目标函数为：

$\cos t(x,z)=-\frac{1}{m}\underset{i=1}{\overset{m}{\Σ}}\underset{k=1}{\overset{d}{\Σ}}\[x_{ik}log\left(z_{ik}\right)+(1-x_{ik})log(1-z_{ik})\]---\left(6\right)$

式中，d代表输入向量的大小，m代表每一份小批量数据的大小；x_ik或z_ik表示小批量数据中的第i个输入或重建的第k个元素；里层的求和是在输入向量的维度上进行的，而外层的求和则是对于小批量数据的每一个样本进行；

在确定好了优化目标函数以后，通过随机梯度下降法来训练网络侧参数，参数更新的规则定义如下，其中η代表学习率：

$W=W-\mathrm{\η}\frac{\∂\cos t(x,z)}{\∂W}---\left(7\right)$

$b_y=b_y-\mathrm{\η}\frac{\∂\cos t(x,z)}{\∂b_y}---\left(8\right)$

$b_z=b_z-\mathrm{\η}\frac{\∂\cos t(x,z)}{\∂b_z}---\left(9\right)$

推导式(7)～(9)的具体形式；

首先，用标量形式重写重构的输出：

${\mathrm{net}}_{ip}^y=\underset{q=1}{\overset{d}{\Σ}}{x}_{iq}{W}_{qp}+b_{yp}---\left(10\right)$

${\mathrm{net}}_{ik}^z=\underset{p=1}{\overset{h}{\Σ}}{W}_{kp}f\left({\mathrm{net}}_{ip}^y\right)+b_{zk}---\left(11\right)$

$z_{ik}=f\left({\mathrm{net}}_{ik}^z\right)=f\left(\underset{p=1}{\overset{h}{\Σ}}{W}_{kp}f\right(\underset{q=1}{\overset{d}{\Σ}}{x}_{iq}{W}_{qp}+b_{yp})+b_{zk})---\left(12\right)$

式中(或)表示对于小批量数据中的第i个样本，第p(或k)个隐层(或输出层)神经元的净输入；

式(1)的一阶和二阶导数如下：

f'(x)＝f(x)[1-f(x)] (13)

f”(x)＝f(x)[1-f(x)][1-2f(x)] (14)

根据式(10)～(12)，就可以计算重构的输出z对于W，b_y，b_z三个参数的偏导数；使用标量形式来表示这些偏导数：

$\frac{\∂z_{ik}}{\∂W_{rs}}=\left\{\begin{array}{cc}{f}^{\′}\left({\mathrm{net}}_{ik}^z\right)\[W_{ks}{f}^{\′}\left({\mathrm{net}}_{is}^y\right)x_{ir}+f\left({\mathrm{net}}_{is}^y\right)\]& k=r\\ f^{\′}\left({\mathrm{net}}_{ik}^z\right)\[W_{ks}{f}^{\′}\left({\mathrm{net}}_{is}^y\right)x_{ir}\]& k\≠r\end{array}\right.---\left(15\right)$

$\frac{\∂z_{ik}}{\∂b_{yr}}=f^{\′}\left({\mathrm{net}}_{ik}^z\right)W_{kr}{f}^{\′}\left({\mathrm{net}}_{ir}^y\right)---\left(16\right)$

$\frac{\∂z_{ik}}{\∂b_{zr}}=f^{\′}\left({\mathrm{net}}_{ik}^z\right)---\left(17\right)$

式中W_rs表示连接第r个输入节点与第s个隐层节点的连接权值；b_yr、b_zr代表第r个隐层、输出层节点的偏置；

综合式(10)～(17)，代价函数式(6)对于W，b_y，b_z三个参数的偏导数表达式为：

$\left\{\begin{array}{c}\frac{\∂\cos t}{\∂W_{rs}}=-\frac{1}{m}\underset{i=1}{\overset{m}{\Σ}}\{\underset{k=1}{\overset{d}{\Σ}}\[\frac{x_{ik}-z_{ik}}{z_{ik}(1-z_{ik})}{f}^{\′}\left({\mathrm{net}}_{ik}^z\right)W_{ks}{f}^{\′}\left({\mathrm{net}}_{is}^y\right)x_{ir}\]+f^{\′}\left({\mathrm{net}}_{ik}^z\right)f^{\′}\left({\mathrm{net}}_{is}^y\right)\}\\ \frac{\∂\cos t}{\∂b_{yr}}=-\frac{1}{m}\underset{i=1}{\overset{m}{\Σ}}\underset{k=1}{\overset{d}{\Σ}}\frac{x_{ik}-z_{ik}}{z_{ik}(1-z_{ik})}{f}^{\′}\left({\mathrm{net}}_{ik}^z\right)W_{kr}{f}^{\′}\left({\mathrm{net}}_{ir}^y\right)\\ \frac{\∂\cos t}{\∂b_{zr}}=-\frac{1}{m}\underset{i=1}{\overset{m}{\Σ}}\underset{k=1}{\overset{d}{\Σ}}\frac{x_{ik}-z_{ik}}{z_{ik}(1-z_{ik})}{f}^{\′}\left({\mathrm{net}}_{ik}^z\right)\end{array}\right.---\left(18\right)$

将这些表达式代入式(7)～(9)，即可得到权值的更新公式；

在训练完模型后，把隐层-重构层挪去，隐层的输出y即作为下一层自动编码机的输入，或者作为支持向量机回归算法的输入。

第4)步支持向量回归算法具体过程为，设训练样本集S为{(x_i,y_i),i＝1,2，...,n}，其中x_i∈R^d为输入向量，y_i为输出向量；对于非线性支持向量回归，通过核函数k(x_i,x_j)将输入样本空间非线性映射到高维特征空间进行线性回归，非线性回归估计函数f(x_i)可表示为：

式(19)中表示向量w与的内积，w为权值向量，为非线性映射函数，b为阈值，b∈R；

为使支持向量回归保持较好的稀疏性，引入不敏感损失函数ε最小化经验风险，得到损失函数为

c(x,y,f)＝|y-f(x)|_ε＝max{0,|y-f(x)|-ε} (20)

若点x的观察值y和回归预测函数值f(x)之间的差别小于ε时，损失为0；

为求解支持向量回归函数中的w与b，在满足不敏感损失函数ε前提下将最小化；引入松弛变量ξ,ξ^*≥0，进一步将式(19)转化为如下优化问题：

$\underset{\mathrm{\ω},\mathrm{\ξ}}{\min }\left\{\frac{1}{2}\right|\left|\mathrm{\ω}\right|{|}^2+C\underset{i=1}{\overset{n}{\Σ}}(\mathrm{\ξ}+{\mathrm{\ξ}}^{*})\}---\left(21\right)$

其中约束条件为：

式(22)中，C为惩罚因子，用于平衡模型复杂项和训练误差项的权重参数，C越大表示对训练误差大于ε的样本惩罚越大；不敏感损失函数ε，规定了回归函数的误差要求；ξ,ξ^*都为松弛变量。

将Lagrange函数引入式(21)，用核函数K(x_i,x_j)将高维空间内积运算通过原空间运算来实现，有

α_i和α_i ^*为拉格朗日算子，(α_i-α_i ^*)≠0的训练样本为支持向量，得到回归函数：

式(24)中n为支持向量个数，核函数采用高斯径向基(RBF)函数：

K(x_i,x_j)＝exp{-|x_i-x_j|²/(2σ²)} (25)。

相比现有技术，本发明具有如下有益效果：

(1)提出了一种基于深度回归的遥感水质监测方法。遥感图像数据呈现复杂的非线性分布，与传统的浅层回归模型相比，运用深度学习方法可以挖掘出隐藏其中的本质规律，消除波段的相关性，降低冗余度，为建立效果良好的水质反演模型提供了有力的保证。

(2)建立了一种适用于小样本的非线性水质反演模型，即利用深度学习得到少量水质遥感数据的本质特征，并对本质特征数据点和对应的水质参数实地监测数据利用支持向量回归算法进行学习，得到适合整个水体的水质反演模型。适用于对地面样本数据少的水质遥感图像进行水质评价和动态监测。

附图说明

图1-自动编码机结构示意图。

图2-堆栈自动编码回归模型示意图。

图3-基于深度回归模型的遥感水质检测流程图。

图4-渭河流域示意图。

图5-渭河流域的遥感影像。

图6-3组测试样本反演结果。

图7-13组监测数据的反演结果。

图8-整个水体氨氮参数的浓度分布图。

具体实施方式

为能全面有效地监测某水域的水质质量，本发明提出了一种深度回归模型遥感水质监测方法，根据堆栈自动编码机的基本原理，构建了基于堆栈自动编码回归(StackedAutoencoder Regression,SAER)的遥感水质监测模型，并用该模型对真实的遥感图像进行了水质参数反演，得到了水质质量分布情况图。

下面结合附图对本发明基于深度回归模型的遥感水质监测方法各步骤进行详细说明，具体可结合图3。

1)将遥感图像数据导入计算机，提取出包含需要监测水体范围作为待处理遥感图像；

2)利用掩模矩阵对待处理遥感图像每个波段进行掩膜，并从每个波段中去掉光谱值为0的背景信息，再做归一化处理，得到选定水体的遥感数据；

3)在监测水体区域选取不同的点作为采样点，对采样点进行需要监测项目的水质参数实地采集，得到采样点的实际水质参数，并确定所有采样点对应的遥感数据；

4)将所有采样点的实际水质参数和采样点的遥感数据作为训练样本，利用多层自动编码机逐层进行训练，以逐层提取遥感水质图像的特征，从而得到每层自动编码机的输出，并作为下一层自动编码机的输入；如此递推，直至最后一层自动编码机训练完毕，得到多层自动编码机的输出；最后将该多层自动编码机的输出利用支持向量机回归算法反演水质参数，得到水质回归模型；

5)利用所构建的水质回归模型对待处理遥感图像进行反演，得到待处理遥感图像中水体需要监测项目的浓度值；根据其处于不同的浓度区间赋予不同的颜色，从而直观显示整个水体水质参数的浓度情况，实现水体水质遥感监测。

在步骤4)中，本发明用到了自动编码机、堆栈自动编码回归和支持向量机回归算法。下面分别对其进行详细说明。

一、自动编码机

自动编码机是一种数据驱动的、非监督地学习数据特征的神经网络模型，其结构如图1所示。可以把它看成是一个输出节点数与输入节点数相等的多层神经网络。在输入-隐层中，模型将计算数据(图中用x表示)在各个隐层节点上的值(图中用y表示)。这些值构成了数据的特征，也称作“表达”。在隐层一输出层中，则根据这一组特征的值重新恢复出输入数据(图中用z表示)。在重建z的过程中，只用了特征y而没有用输入x，因此，所有z中的信息都只由y得来；正因如此，y中是否保留了足够多的输入信息即可通过z的重建质量来评价。因此，提取特征的过程就可以被这么设定：训练这样一个神经网络，使得它重建得到的输出与输入尽可能地一致。这样的模型使得可以根据数据本身来无监督地训练网络。

假设自动编码机包含d个输入神经元，隐藏层y则含有h个神经元。重建层z含有与输入层相同的神经元数目。网络中各个神经元的激活函数设定为sigmoid函数:

$f\left(x\right)=sigmoid\left(x\right)=\frac{1}{1+e^{-x}}---\left(1\right)$

在训练过程中，模型首先将输入向量x(x∈R^d)映射到隐层中，以得到隐层特征y(y∈R^h)这一步所对应的网络如图1中框线内部分所示。称这部分网络为“编码器”。然后y被下一层网络(称为“解码器”)映射到输出层，称这一层的输出为“重建的输出”z(z∈R^d)。这两个步骤可以用如下两个公式表示:

y＝sigmoid(W_yx+b_y) (2)

z＝sigmoid(W_zy+b_z) (3)

其中W_y，W_z分别是输入一隐层，隐层一输出层的连接权值矩阵；层与输出层每个神经元的偏置。为了减少模型中需要训练的参数，b_y，b_z是隐给出如下约束:

W_y＝W_z＝W (4)

此时称这样的自动编码机具有“约束的权值”。因此，对于这样的自动编码机，训练就意味着学习模型中包含的三组参数W，b_y，b_z。

模型的训练则是一个对该目标函数的最优化过程。即调节W,b_y,b_z来最小化重建的输出与输入之间的差异。所以可以确定这样的目标函数:

$\begin{array}{cc}idealreconstruction=\arg & \underset{W,b_y,b_z}{\min }\[\cos t(x,z)\]\end{array}---\left(5\right)$

这里的z在x给定的情况下取决于模型的参数W,b_y,b_z。cost(x,z)是重建与输入的差异的度量，本文使用这两者之间的互熵来度量这种差异。在实际的实现中，权值的更新是基于小批量数据的。即，累积一小批数据的差异后再更新权值。这样的做法有助于提高算法在处理大规模数据时的学习速率。即定义的目标函数为：

$\cos t(x,z)=-\frac{1}{m}\underset{i=1}{\overset{m}{\Σ}}\underset{k=1}{\overset{d}{\Σ}}\[x_{ik}log\left(z_{ik}\right)+(1-x_{ik})log(1-z_{ik})\]---\left(6\right)$

式中，d代表输入向量的大小，m代表每一份小批量数据的大小。x_ik(或z_ik)表示小批量数据中的第i个输入(或重建)的第k个元素。里层的求和是在输入向量的维度上进行的；而外层的求和则是对于小批量数据的每一个样本进行。

在确定好了优化目标及罚函数以后，可以通过随机梯度下降法来训练网络侧参数。参数更新的规则定义如下(其中η代表学习率):

$W=W-\mathrm{\η}\frac{\∂\cos t(x,z)}{\∂W}---\left(7\right)$

$b_y=b_y-\mathrm{\η}\frac{\∂\cos t(x,z)}{\∂b_y}---\left(8\right)$

$b_z=b_z-\mathrm{\η}\frac{\∂\cos t(x,z)}{\∂b_z}---\left(9\right)$

现在推导式(7)～(9)的具体形式。这之中最重要的部分即是对罚函数对于三个参数W，b_y，b_z与偏导数形式的推导。

首先，用标量形式重写重建的输出:

${\mathrm{net}}_{ip}^y=\underset{q=1}{\overset{d}{\Σ}}{x}_{iq}{W}_{qp}+b_{yp}---\left(10\right)$

${\mathrm{net}}_{ik}^z=\underset{p=1}{\overset{h}{\Σ}}{W}_{kp}f\left({\mathrm{net}}_{ip}^y\right)+b_{zk}---\left(11\right)$

$z_{ik}=f\left({\mathrm{net}}_{ik}^z\right)=f\left(\underset{p=1}{\overset{h}{\Σ}}{W}_{kp}f\right(\underset{q=1}{\overset{d}{\Σ}}{x}_{iq}{W}_{qp}+b_{yp})+b_{zk})---\left(12\right)$

式中(或)表示对于小批量数据中的第i个样本，第p(或k)个隐层(或输出层)神经元的净输入。

式(1)的一阶和二阶导数如下:

f'(x)＝f(x)[1-f(x)] (13)

f”(x)＝f(x)[1-f(x)][1-2f(x)] (14)

根据式(10)～(12)，就可以计算重建的输出z对于W，b_y，b_z三个参数的偏导数。为了简化标记，依然使用标量形式来表示这些偏导数:

$\frac{\∂z_{ik}}{\∂W_{rs}}=\left\{\begin{array}{cc}{f}^{\′}\left({\mathrm{net}}_{ik}^z\right)\[W_{ks}{f}^{\′}\left({\mathrm{net}}_{is}^y\right)x_{ir}+f\left({\mathrm{net}}_{is}^y\right)\]& k=r\\ f^{\′}\left({\mathrm{net}}_{ik}^z\right)\[W_{ks}{f}^{\′}\left({\mathrm{net}}_{is}^y\right)x_{ir}\]& k\≠r\end{array}\right.---\left(15\right)$

$\frac{\∂z_{ik}}{\∂b_{yr}}=f^{\′}\left({\mathrm{net}}_{ik}^z\right)W_{kr}{f}^{\′}\left({\mathrm{net}}_{ir}^y\right)---\left(16\right)$

$\frac{\∂z_{ik}}{\∂b_{zr}}=f^{\′}\left({\mathrm{net}}_{ik}^z\right)---\left(17\right)$

式中W_rs表示连接第r个输入节点与第s个隐层节点的连接权值。b_yr(或b_zr)代表第r个隐层(输出层)节点的偏置。

综合式(10)～(17}，代价函数(式(6))对于W，b_y，b_z三个参数的偏导数表达式为:

$\left\{\begin{array}{c}\frac{\∂\cos t}{\∂W_{rs}}=-\frac{1}{m}\underset{i=1}{\overset{m}{\Σ}}\{\underset{k=1}{\overset{d}{\Σ}}\[\frac{x_{ik}-z_{ik}}{z_{ik}(1-z_{ik})}{f}^{\′}\left({\mathrm{net}}_{ik}^z\right)W_{ks}{f}^{\′}\left({\mathrm{net}}_{is}^y\right)x_{ir}\]+f^{\′}\left({\mathrm{net}}_{ik}^z\right)f\left({\mathrm{net}}_{is}^y\right)\}\\ \frac{\∂\cos t}{\∂b_{yr}}=-\frac{1}{m}\underset{i=1}{\overset{m}{\Σ}}\underset{k=1}{\overset{d}{\Σ}}\frac{x_{ik}-z_{ik}}{z_{ik}(1-z_{ik})}{f}^{\′}\left({\mathrm{net}}_{ik}^z\right)W_{kr}{f}^{\′}\left({\mathrm{net}}_{ir}^y\right)\\ \frac{\∂\cos t}{\∂b_{zr}}=-\frac{1}{m}\underset{i=1}{\overset{m}{\Σ}}\underset{k=1}{\overset{d}{\Σ}}\frac{x_{ik}-z_{ik}}{z_{ik}(1-z_{ik})}{f}^{\′}\left({\mathrm{net}}_{ik}^z\right)\end{array}\right.---\left(18\right)$

将这些表达式代入式(7)～(9}，即可得到权值的更新公式。

在训练完模型后，把隐层一输出层挪去，隐层的输出y即可作为下一层某个分类器的输入，或者输出给下一级的自动编码机。

自动编码机这一模型的优势即在于这种基于重建的训练。值得注意的是，在重建的过程中，只利用了隐藏层的特征y，而y是在编码机中由输入得到的。如果模型可以由y成功地重建输入，那就意味着y这一组特征包含了来自输入的足够充分的信息。因此，学习所得到的那个由权值和偏置所定义的非线性变换，即可被看作一个很好的非线性特征提取器。

二、堆栈自动编码回归

2.1堆栈自动编码回归模型

堆栈自动编码回归是指利用普通的自动编码机层叠后，再在顶层加上一个单层的支持向量机回归得到的整个深层神经网络回归模型。这样含有多层自动编码机的网络，可以逐层地提取遥感水质图像的特征，得到深层更抽象更具可分性的特征，最后利用支持向量回归反演水质参数。

如图2示出了一个含有3个堆叠起来的自动编码机与一个支持向量回归所构成的堆栈自动编码回归模型。该网络一共含有5层，将初始的输入数据转化为了输出层水质参数的浓度值。第一层自动编码机将第0层的输入映射成为第1层上的特征。训练这层网络的方法与前面所介绍的训练单层的自动编码机方法相同，训练后只保留了编码机部分的网络。当训练完第一层自动编码机后，下一层的自动编码机将第一层的特征映射为第二层上的特征。训练方法也与单层自动编码机一样，只是这时的输入改为前一层己训练网络的输出。如此递推，直至将每一层自动编码机训练完毕。

如此逐层地堆叠自动编码机可以在每一层上都最小化信息的损失。同时，提取出的特征逐层地变得更抽象、具有更好的不变性。这也是使用堆栈自动编码机用于提取深层特征的出发点所在。

2.2支持向量回归

Vapnik等人于1963年提出了支持向量机(Support Vector Machines,SVM)原理。SVM以VC维理论(统计学)和结构风险最小化为基础，凭借有限的样本信息在训练样本的学习精度和正确识别任意样本的能力，即模型复杂度和学习能力之间寻求最好的推广能力。SVM原理是在统计学习理论的基础上提出的一种机器学习方法，在许多领域都得到了推广和应用。

支持向量回归(Support Vector Regression,SVR)是在SVM原理基础上建立的一种回归技术，在时间序列的预测、非线性建模与预测及优化控制等方面，具有泛化能力强和预测准确度高等优点。设训练样本集S为{(x_i,y_i),i＝1,2，...,n}，其中x_i∈R^d为输入向量，y_i为输出向量。对于非线性支持向量回归，通过核函数将输入样本空间非线性映射到高维特征空间进行线性回归，非线性回归估计函数f(x_i)可表示为：

式(19)中表示向量w与的内积，w为权值向量，为非线性映射函数，b为阈值，b∈R。

为使支持向量回归保持较好的稀疏性，引入不敏感损失函数ε最小化经验风险，得到损失函数为

c(x,y,f)＝|y-f(x)|_ε＝max{0,|y-f(x)|-ε} (20)

若点x的观察值y和回归预测函数值f(x)之间的差别小于ε时，损失为0。

为求解支持向量回归函数中的w与b，在满足不敏感损失函数ε前提下将最小化。引入松弛变量ξ,ξ^*≥0，进一步将式(19)转化为如下优化问题：

$\underset{\mathrm{\ω},\mathrm{\ξ}}{\min }\left\{\frac{1}{2}\right|\left|\mathrm{\ω}\right|{|}^2+C\underset{i=1}{\overset{n}{\Σ}}(\mathrm{\ξ}+{\mathrm{\ξ}}^{*})\}---\left(21\right)$

其中约束条件为：

式(22)中，C为惩罚因子，用于平衡模型复杂项和训练误差项的权重参数，C越大表示对训练误差大于ε的样本惩罚越大。不敏感损失函数ε，规定了回归函数的误差要求。

将Lagrange函数引入式(21)，用核函数K(x_i,x_j)将高维空间内积运算通过原空间运算来实现，有

α_i和α_i ^*为拉格朗日算子，(α_i-α_i ^*)≠0的训练样本为支持向量，得到回归函数：

式(24)中n为支持向量个数，核函数采用高斯径向基(RBF)函数：

K(x_i,x_j)＝exp{-|x_i-x_j|²/(2σ²)} (25)

以下结合具体实施例对本发明所涉及监测方法进行详细说明。

渭河流域是陕西省的人口密集地带和重要工业基地，流域内集中了全省61％的人口和上千家工矿企业，其工业总产值占全省82％，城市生活污水和工业污水多以排入支流后排入渭河。渭河沿岸工业企业集中，城镇人口密集，自然是工业废水和城市生活污水的主要产区，来自两岸诸多点源的废污水，将渭河作为唯一的排泄通道，自上而下排量沿程不断增加，水污染超标程度不断增加。根据环境监测部门监测显示，1998年，全省废污水年排放量为7.03亿m²，其中工业废水3.45亿m²，城镇生活污水3.58亿m²，而渭河干流两岸纳入的废水污水量达到5.52亿m²，占全省污水排放量的78.5％。在渭河全段尤其以洋河入渭河至零河入渭河口段，即西安市区所在区段，年排放量达1.97亿m²，占渭河全程年排放量的40.3％其次为林家村至蔡家坡(即宝鸡市所在区段)，年废污水量为1.25亿m²，占渭河全程年排放量的25.8％，咸阳公路桥至洋河入渭河口段(即咸阳市所在区段)，年废污水排放量为1.05亿m²，占渭河全程年排放量的21.2％；渭南市区所在河段，年废污水排放量为0.6亿m²，占渭河全程排放量的12.2。另外据2005年统计，陕西省境内渭河干支流排污口增至290个，这52条大小河流，年废污水排放量达8.07亿m²，占全省废污水排放量的81.6％，与1998年相比，年排放量增加2.54亿m²，7年中平均每年以0.36亿m²的速度递增。渭河流域水土流失严重，大面积的水土流失对渭河水质造成严重污染。农药、化肥的使用，大部分残留随着农田退水和地表径流进入水体，造成污染。一些工矿企业的工业废弃物和垃圾随意堆积在河滩或直接倾入渭河中，造成了水体污染。多年调查统计数据表明:渭河是陕西省多条大河流中接纳废污水量最多的河流，是污染最严重的河流。图4为渭河流域示意图。

利用本发明对2006年渭河陕西段水质遥感图像数据进行水质监测。2006年渭河陕西段水质遥感图像有4个波段，共3451×3654个数据。选定[1260:2100，400:2300]范围以包含所有监测点，包括耿镇和成阳等地段，共841×1901个数据，其中包含60944个水体数据，其遥感影像如图5所示。同时获取了2000-2006年间渭河陕西段部分监测点共13组水质遥感数据和对应的实地监测数据，数据分布在渭河陕西段9个监测断面，水质参数指标包括各监测点的高锰酸盐指数(COD_mn)、氨氮(NH₃-N)和重铬酸盐指数(COD_cr)等，由于渭河陕西段的污染情况主要为有机物污染，因此在水质参数指标上首选耗氧有机物指标及氧平衡指标，本例选择氨氮(NH₃-N)进行数据分析和实验。

具体步骤为：

1)将遥感图像数据读入计算机，并提取出包含水体的矩形范围：计算机读入2006年渭河陕西段水质遥感图像，选定[1260:2100，400:2300]的矩形范围，得到包含所有监测点的遥感图像。

2)通过ENVI软件平台，人工解译生成掩模矩阵。利用掩模矩阵对遥感图像每个波段进行掩膜，并从每个波段中去掉光谱值为0的背景信息，再做归一化处理，得到选定水体60944行4列的遥感数据矩阵B＝{x₁，x₂，…，x_i，…，x₆₀₉₄₄}^T，i∈[1，60944]，x_i是第i个数据点4个波段的光谱值构成的列向量。

3)采用获取的13组水质遥感数据和对应的实地监测数据，从中选取10组数据作为训练样本，剩下的3组作为测试样本。考虑到数据的大小不一致，为了使训练样本和测试样本的数据分布较为均匀，尽量使训练样本中的浓度数据能够包含较大，中等及较小的数据，而测试样本也符合这种分布。

4)通过训练样本对本发明的深度回归模型进行建模，在建模时，网络层数设置为4，其中第一层为输入层，中间二层为隐含层，最C后一层为输出层。两隐含层的节点数分别为6和3，输出层的惩罚参数和核参数σ分别设置为50和5。

5)根据构建的模型，对测试样本进行反演，其反演结果如图6所示，图7为13组数据的反演结果，根据图5和图6的反演结果可知，该模型能够达到较理想的效果。其中：+：代表原数据，即真实值；○：代表预测数据。

6)利用所构建的模型对渭河流域陕西段进行反演，将整个水体水质参数的浓度值根据其处于不同的区间赋予不同的颜色，浓度值在0～4范围内的地段以绿色显示；浓度值在4～8范围内的地段以红色显示；浓度值在8～12范围内的地段以蓝色显示；浓度值超过12的地段以紫色显示；从而直观显示整个水体水质参数的浓度情况，实现水体的水质评价和监测。图8为整个水体氨氮参数的浓度分布图。

7)为了更好的查看氨氮参数的浓度分布，可以将各个浓度区间内的氨氮分别通过图进行显示。

根据实验结果可知，咸阳段的氨氮浓度值最高，这是由于沿河两岸城镇居民数量居多，大量生活污水直接排入地表水体；农业面源也是主要污染源，大量使用化肥农药导致地表水体中氨氮浓度的升高；拥有多家工业企业(如咸阳造纸厂、兴平造纸厂)排放的污水，也造成了氨氮浓度较大。与此同时，河流上游地区取水灌溉，使河流的流量降低，河流的同化、自净能力减弱。

最后需要说明的是，本发明的上述实例仅仅是为说明本发明所作的举例，而并非是对本发明的实施方式的限定。尽管申请人参照较佳实施例对本发明进行了详细说明，对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其他不同形式的变化和变动。这里无法对所有的实施方式予以穷举。凡是属于本发明的技术方案所引申出的显而易见的变化或变动仍处于本发明的保护范围之列。

Claims (4)

1.基于深度回归模型的遥感水质监测方法，其特征在于：步骤如下，

1)将遥感图像数据导入计算机，提取出包含需要监测水体的矩形范围作为待处理遥感图像；

2)利用掩模矩阵对待处理遥感图像每个波段进行掩膜，并从每个波段中去掉光谱值为0的背景信息，再做归一化处理，得到选定水体的遥感数据；

3)在监测水体区域选取不同的点作为采样点，对采样点进行需要监测项目的水质参数实地采集，得到采样点的实际水质参数，并确定所有采样点对应的遥感数据；

4)将所有采样点的实际水质参数和采样点的遥感数据作为训练样本，利用多层自动编码机逐层进行训练，以逐层提取遥感水质图像的特征，从而得到每层自动编码机的输出，并作为下一层自动编码机的输入；如此递推，直至最后一层自动编码机训练完毕，得到多层自动编码机的输出；最后将该多层自动编码机的输出利用支持向量机回归算法反演水质参数，得到水质回归模型；

5)利用所构建的水质回归模型对待处理遥感图像进行反演，得到待处理遥感图像中水体需要监测项目的浓度值；根据其处于不同的浓度区间赋予不同的颜色，从而直观显示整个水体水质参数的浓度情况，实现水体水质遥感监测。

2.根据权利要求1所述的基于深度回归模型的遥感水质监测方法，其特征在于：第3)步采样点的实际水质参数包含监测项浓度较大，中等及较小的数据分布。

3.根据权利要求1所述的基于深度回归模型的遥感水质监测方法，其特征在于：第4)步多层自动编码机中每层自动编码机按如下方法训练，

假设单层自动编码机输入层x包含d个输入神经元，隐层y含有h个神经元；重构层z含有与输入层相同的神经元数目；网络中各个神经元的激活函数设定为sigmoid函数：

$f\left(x\right)=sigmoid\left(x\right)=\frac{1}{1+e^{-x}}---\left(1\right)$

在训练过程中，模型首先将输入向量x(x∈R^d)映射到隐层中，以得到隐层特征y(y∈R^h)，然后将隐层特征y(y∈R^h)映射到重构层，得到重构层特征z(z∈R^d)；这两个映射环节分别用如下两个公式表示：

y＝sigmoid(W_yx+b_y) (2)

z＝sigmoid(W_zy+b_z) (3)

其中W_y，W_z分别是输入层-隐层以及隐层-重构层的连接权值矩阵；b_y，b_z是隐层与重构层每个神经元的偏置；为了减少模型中需要训练的参数，给出如下约束：

W_y＝W_z＝W (4)

对于这样的自动编码机，训练就意味着学习模型中包含的三组参数W，b_y，b_z；

模型的训练则是一个对该目标函数的最优化过程，即调节W,b_y,b_z来最小化重构的输出与输入之间的差异；故确定这样的目标函数：

$idealreconstruction=\arg \underset{W,b_y,b_z}{min}\[\cos t(x,z)\]---\left(5\right)$

这里的z在x给定的情况下取决于模型的参数W,b_y,b_z；cost(x,z)是重建与输入的差异的度量，定义目标函数为：

$\cos t(x,z)=-\frac{1}{m}\underset{i=1}{\overset{m}{\Σ}}\underset{k=1}{\overset{d}{\Σ}}\[x_{ik}log\left(z_{ik}\right)+(1-x_{ik})log(1-z_{ik})\]---\left(6\right)$

式中，d代表输入向量的大小，m代表每一份小批量数据的大小；x_ik或z_ik表示小批量数据中的第i个输入或重建的第k个元素；里层的求和是在输入向量的维度上进行的，而外层的求和则是对于小批量数据的每一个样本进行；

在确定好了优化目标函数以后，通过随机梯度下降法来训练网络侧参数，参数更新的规则定义如下，其中η代表学习率：

$W=W-\mathrm{\η}\frac{\∂\cos t(x,z)}{\∂W}---\left(7\right)$

$b_y=b_y-\mathrm{\η}\frac{\∂\cos t(x,z)}{\∂b_y}---\left(8\right)$

$b_z=b_z-\mathrm{\η}\frac{\∂\cos t(x,z)}{\∂b_z}---\left(9\right)$

推导式(7)～(9)的具体形式；

首先，用标量形式重写重构的输出：

${\mathrm{net}}_{ip}^y=\underset{q=1}{\overset{d}{\mathrm{\Σ}}}{x}_{iq}{W}_{qp}+b_{yp}---\left(10\right)$

${\mathrm{net}}_{ik}^z=\underset{p=1}{\overset{h}{\mathrm{\Σ}}}{W}_{kp}f\left({\mathrm{net}}_{ip}^y\right)+b_{zk}---\left(11\right)$

$z_{ik}=f\left({\mathrm{net}}_{ik}^z\right)=f\left(\underset{p=1}{\overset{h}{\Σ}}{W}_{kp}f\right(\underset{q=1}{\overset{d}{\Σ}}{x}_{iq}{W}_{qp}+b_{yp})+b_{zk})---\left(12\right)$

式中(或)表示对于小批量数据中的第i个样本，第p(或k)个隐层(或输出层)神经元的净输入；

式(1)的一阶和二阶导数如下：

f'(x)＝f(x)[1-f(x)] (13)

f”(x)＝f(x)[1-f(x)][1-2f(x)] (14)

根据式(10)～(12)，就可以计算重构的输出z对于W，b_y，b_z三个参数的偏导数；使用标量形式来表示这些偏导数：

$\frac{\∂z_{ik}}{\∂W_{rs}}=\left\{\begin{array}{cc}{f}^{\′}\left({\mathrm{net}}_{ik}^z\right)\[W_{ks}{f}^{\′}\left({\mathrm{net}}_{is}^y\right)x_{ir}+f\left({\mathrm{net}}_{is}^y\right)\]& k=r\\ f^{\′}\left({\mathrm{net}}_{ik}^z\right)\[W_{ks}{f}^{\′}\left({\mathrm{net}}_{is}^y\right)x_{ir}\]& k\≠r\end{array}\right.---\left(15\right)$

$\frac{\∂z_{ik}}{\∂b_{yr}}=f^{\′}\left({\mathrm{net}}_{ik}^z\right)W_{kr}{f}^{\′}\left({\mathrm{net}}_{ir}^y\right)---\left(16\right)$

$\frac{\∂z_{ik}}{\∂b_{zr}}=f^{\′}\left({\mathrm{net}}_{ik}^z\right)---\left(17\right)$

式中W_rs表示连接第r个输入节点与第s个隐层节点的连接权值；b_yr、b_zr代表第r个隐层、输出层节点的偏置；

综合式(10)～(17)，代价函数式(6)对于W，b_y，b_z三个参数的偏导数表达式为：

$\left\{\begin{array}{c}\frac{\∂\cos t}{\∂W_{rs}}=-\frac{1}{m}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\left\{\underset{k=1}{\overset{d}{\mathrm{\Σ}}}\[\frac{x_{ik}-z_{ik}}{z_{ik}\left(1-z_{ik}\right)}{f}^{\′}\left({\mathrm{net}}_{ik}^z\right)W_{ks}{f}^{\′}\left({\mathrm{net}}_{is}^y\right)x_{ir}\]+f^{\′}\left({\mathrm{net}}_{ik}^z\right)f\left({\mathrm{net}}_{is}^y\right)\right\}\\ \frac{\∂\cos t}{\∂b_{yr}}=-\frac{1}{m}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{k=1}{\overset{d}{\mathrm{\Σ}}}\frac{x_{ik}-z_{ik}}{z_{ik}\left(1-z_{ik}\right)}{f}^{\′}\left({\mathrm{net}}_{ik}^z\right)W_{kr}{f}^{\′}\left({\mathrm{net}}_{ir}^y\right)\\ \frac{\∂\cos t}{\∂b_{zr}}=-\frac{1}{m}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{k=1}{\overset{d}{\mathrm{\Σ}}}\frac{x_{ik}-z_{ik}}{z_{ik}\left(1-z_{ik}\right)}{f}^{\′}\left({\mathrm{net}}_{ik}^z\right)\end{array}\right.---\left(18\right)$

将这些表达式代入式(7)～(9)，即可得到权值的更新公式；

在训练完模型后，把隐层-重构层挪去，隐层的输出y即作为下一层自动编码机的输入，或者作为支持向量机回归算法的输入。

4.根据权利要求1所述的基于深度回归模型的遥感水质监测方法，其特征在于：第4)步支持向量回归算法具体过程为，设训练样本集S为{(x_i,y_i),i＝1,2，...,n}，其中x_i∈R^d为输入向量，y_i为输出向量；对于非线性支持向量回归，通过核函数k(x_i,x_j)将输入样本空间非线性映射到高维特征空间进行线性回归，非线性回归估计函数f(x_i)可表示为：

式(19)中表示向量w与的内积，w为权值向量，为非线性映射函数，b为阈值，b∈R；

为使支持向量回归保持较好的稀疏性，引入不敏感损失函数ε最小化经验风险，得到损失函数为

c(x,y,f)＝|y-f(x)|_ε＝max{0,|y-f(x)|-ε} (20)

若点x的观察值y和回归预测函数值f(x)之间的差别小于ε时，损失为0；

为求解支持向量回归函数中的w与b，在满足不敏感损失函数ε前提下将最小化；引入松弛变量ξ,ξ^*≥0，进一步将式(19)转化为如下优化问题：

$\underset{\mathrm{\ω},\mathrm{\ξ}}{min}\left\{\frac{1}{2}\right|\left|\mathrm{\ω}\right|{|}^2+C\underset{i=1}{\overset{n}{\Σ}}(\mathrm{\ξ}+{\mathrm{\ξ}}^{*})\}---\left(21\right)$

其中约束条件为：

式(22)中，C为惩罚因子，用于平衡模型复杂项和训练误差项的权重参数，C越大表示对训练误差大于ε的样本惩罚越大；不敏感损失函数ε，规定了回归函数的误差要求；ξ,ξ^*都为松弛变量。

将Lagrange函数引入式(21)，用核函数K(x_i,x_j)将高维空间内积运算通过原空间运算来实现，有

α_i和α_i ^*为拉格朗日算子，(α_i-α_i ^*)≠0的训练样本为支持向量，得到回归函数：

式(24)中n为支持向量个数，核函数采用高斯径向基(RBF)函数：

K(x_i,x_j)＝exp{-|x_i-x_j|²/(2σ²)} (25)。