CN115951025B - 一种水体总磷浓度的反演方法 - Google Patents

Abstract

本发明提出一种水体总磷浓度的反演方法，该方法可以降低环境气象要素对总磷浓度反演的干扰，提高总磷浓度反演精度。具体包括：S1、采集测试点的水体反射光谱、气象要素数据；S2、对所述水体反射光谱和气象要素数据进行预处理；S3、将所述预处理后的数据输入改进的随机森林，获得测试点的水体总磷浓度。本发明提出的水体总磷浓度反演方法考虑了环境气象要素对总磷浓度反演的影响，适用性较高，具有高精度、高灵活性和强稳健性，可作为水体总磷监测的辅助手段，提高水体总磷监测的广度和质量，节省人力、财力，提高效率。

Description

一种水体总磷浓度的反演方法

技术领域

本发明属于水质监测领域技术领域，尤其涉及一种水体总磷浓度的反演方法。

背景技术

水是生命之源，水域环境的健康与人类的命运息息相关。总磷(TP)浓度是水体污染程度的主要评价指标，总磷浓度的准确获取，能够为磷类污染物精准防治和水体修复提供重要的技术支撑。

目前，水体总磷浓度监测有人工移动监测、固定站点监测和遥感监测等方式。人工移动监测，依赖野外采样和实验室分析，工作量大、效率低，只能对特定时间段某一有限区域进行监测，难以反映全区域水质情况；固定站点监测，通过特定波长处吸光度的换算，可以实现定点实时在线监测，无法满足大范围同步动态监测的需求；遥感监测，包括可见光、红外光谱等监测手段，具有获取数据快、覆盖范围广和同步观测等优势。光谱遥感能够获取丰富精细的光谱信息，对水体总磷、总氮等元素感知度更灵敏、估算精度更高，已成为当前最有效的遥感监测手段。卫星光谱遥感成本较高、分辨率较低、很难推广应用。地面光谱遥感具有成本低、操作简单、受天气影响较小、无需大气校正等特点，是水体参数定量化研究的有效获取手段。

研究表明风速的变化会打破底泥与水之间的平衡关系，引起底泥再悬浮，进而导致底泥中的磷释放进入水中。不同类型的风浪会导致不同程度的沉积物再悬浮和磷释放，强风和中风条件下，悬浮固体和总磷均显著增加，且从沉积物中释放的磷主要以颗粒形式存在。短期强风条件导致更多的总磷释放，但风停止后总磷随悬浮固体迅速沉降到底部，并且没有促进藻类的生长。长期中风条件下，悬浮颗粒磷被磷酸酶水解为可溶性反应性磷，并促进藻类生长。藻类的生长反过来会分泌更多的磷酸酶，并促进颗粒物磷的再生。有研究认为气象因素可以解释太湖总磷浓度的58％年际变化，其中平均风速为主导气象因素，较低的风速和较长的低风持续时间导致湖底附近的溶解氧浓度较低，通过反硝化作用可以增加沉积物中磷的释放，增加氮损失，导致表层总磷浓度变高。

温度的变化会响水体中藻类、水生植物及微生物的生长，此类物质的生长消亡会影响总磷浓度的变化，进而对总磷的季节变化也有一定的影响。水体中各形态磷浓度在水平空间上沿水流方向递减，在垂直方向上各形态磷浓度向下递增，并受水温分层的影响。当湖水温度升高时，藻类植物繁殖加快，减少了上覆水体中磷的浓度，增加了对磷的需求，从而使平衡向着有利于磷释放的一方移动，促进沉积物中磷的释放，同时水生生物活动加剧，也促进了沉积物磷的释放。研究表明黑臭河沉积物中可溶性磷的浓度随着温度的升高而显着增加，与15℃时的浓度相比，在25℃和35℃时分别增加28％和87％，羟基氧化铁的还原和溶解性有机物竞争性吸收矿物质是沉积物中磷释放原因，而温度与这两者呈现正相关。

但是，当前总磷浓度遥感反演使用的数据多局限于水体反射光谱，结合环境气象要素的研究甚少，而且水体总磷浓度受到环境气象要素的影响，具有明显的区域性差异。因此，迫切需要一种考虑环境气象要素影响的水体总磷浓度遥感监测方法。

发明内容

基于上述存在的问题，本发明的目的是提出一种水体总磷浓度的反演方法，可作为传统水质监测手段的辅助方法，节省人力、财力，提高水体总磷浓度监测的广度、质量和效率。

本发明提供的一种水体总磷浓度的反演方法，包括：

S1、采集测试点的水体反射光谱、气象要素数据；

S2、对所述水体反射光谱和气象要素数据进行预处理；

S3、将所述预处理后的数据输入改进的随机森林，获得测试点的水体总磷浓度。

进一步地，所述气象数据为温度和风速。

进一步地，所述步骤S2中，所述预处理的方法包括：

S21、对同一地点采集的水体反射光谱构成的曲线进行平滑降噪，并获取A个光谱波段特征；

S22、对所述水体反射光谱进行归一化，公式如下：

式中，X_norm为归一化后的数值，X为水体反射光谱和气象要素，min表示最小值，max表示最大值，将水体反射光谱的归一化结果和气象要素的归一化结果拼接后构成所述同一地点的特征向量。

进一步地，所述步骤S2中，对同一地点采集多次数据，使用平均值作为最终数据。

进一步地，通过Savitzky-Golay平滑法进行所述平滑降噪。

进一步地，所述步骤S3中，改进的随机森林，利用相关系数对决策树结果进行加权以获取最终总磷预测结果并评价。

进一步地，所述步骤S3中，所述改进的随机森林的训练方法包括：

S31、训练样本数据集X，样本总数为p，样本包含光谱波段特征和气象要素特征，样本对应总磷浓度Y，A为光谱波段特征的数量，B为气象要素特征的数量，计算每个特征的皮尔逊相关系数，得到皮尔逊相关系数集合r＝{r₁,r₂,…,r_m},m∈{1,2,…,A+B}；

S32、使用训练样本数据集X与总磷浓度Y作为输入，采用bootstrap抽样技术从原始训练集随机产生t个训练子集，利用每个训练子集，生成对应的决策树{h₁(x),h₂(x),…,h_t(x)}；

S33、从光谱波段特征中随机选择q(q≤A)个特征与B个气象要素特征共同作为候选特征，在决策树的每个节点按照预定的规则选择最优属性进行分裂,每个决策树都得到最大限度的生长，过程中完全分裂不剪枝；

S34、计算每棵决策树h_i(x)的权值p_i，计算公式为：

其中，q+B为构建该决策树使用到的特征数量，R_ij为决策树h_i(x)的第j个特征所对应的相关系数；

然后对p_i进行归一化，公式为：

S35、对于特征向量X，总磷浓度预测的结果表示为：

其中，

为决策树i的预测结果。

进一步地，计算每个特征的皮尔逊相关系数的方法为：

训练样本数据集X，样本总数为p，特征数量为m，样本对应总磷浓度Y，根据公式

其中r_j为第j个特征与总磷的相关系数，/>

和/>

分别是X中的第j个特征和Y的统计数据的平均值。

进一步地，对反演模型性能进行评价的指标为决定系数R²，选择R²最高时的参数作为最终的结果，R²的公式为：

式中：p代表样本数；y_i表示第i个样本总磷实测值；

表示第i个样本总磷估计值；

表示p个样本的总磷实测值的平均值。

进一步地，所述训练样本数据集通过步骤S1、S2构建。

本发明所述的一种水体总磷浓度的反演方法，具有如下有益效果：

(1)本发明对光谱进行了平滑、归一化，降低了背景噪声对总磷反演的影响。

(2)本发明针对当前水质监测局限于水体反射光谱，而忽略环境要素对反演精度的影响的问题，在建模过程中考虑环境要素，解决了环境要素对总磷浓度反演效果的影响，所构建模型具有预测精度高，对环境差异具有较好适应能力等特点。

(3)本发明针对传统随机森林中不同泛化能力决策树具有相同权重的问题，将提出的改进随机森林方法应用于光谱数据处理领域，实现了总磷浓度的预测，提高了总磷浓度反演精度和效率。

(4)本发明的方法适用性较高，具有高精度、高灵活性和强稳健性，为水体总磷浓度监测提供了新思路，提高总磷监测的广度、质量和效率。

附图说明

为了更清楚地说明本公开实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本公开的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些图获得其他的附图。

图1为本发明一个实施例的反演方法的流程示意图；

图2为本发明一个实施例的改进的随机森林的训练流程示意图。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整的描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。

本发明的思路是：首先采集样本点的水体反射光谱、环境气象要素和水样。其中水样需要及时送至实验室进行化学分析，获取水体总磷浓度。然后对样本的反射光谱和环境气象要素进行预处理，使用预处理结果和总磷浓度作为输入，基于改进的随机森林算法建立总磷反演模型。模型建立完成后，按照相同流程对未知水体反射光谱和环境气象要素进行预处理，输入改进随机森林反演模型即可获取水体总磷浓度。

本发明提出的水体总磷浓度的反演方法，如图1所示，包括：

S1、采集测试点的水体反射光谱、气象要素数据；

S2、将所述水体反射光谱和气象数据进行预处理；

S3、将所述预处理后的数据输入改进的随机森林，获得测试点的水体总磷浓度。

在步骤S1中，光谱中包括多个波段，如270个波段，气象数据包括风速、温度等，一般来说波段数量远远大于气象因素，因为对总磷浓度有影响的气象要素是有限的。

训练改进的随机森林所需的数据集可以在水体处进行采集。比如，采集多个样点的水体反射光谱、环境气象要素和水样，其中同一样本采样点获取多条光谱曲线，水样当天送至实验室内进行分析获取总磷浓度。

在步骤S2中，要在对获取的数据(包括测试数据和待分析数据)进行预处理和归一化。首先，将同一样本点获取的多条光谱曲线，使用Savitzky-Golay平滑方法(当然，还有很多方法可以平滑去噪)对每条光谱曲线进行平滑，取这些曲线的平均值作为该样本的最终水体反射光谱曲线。

归一化处理使用线性函数归一化方法，该方法实现对原始数据的等比例缩放,线性函数将原始数据线性化的方法转换到[0,1]的范围，归一化公式如下：

式中X_norm为归一化后的数据，X分别为水体反射光谱P、气象要素(例如，风速U、温度T，下面以此为例进行描述，但这并不是对本发明的限定)，X_max分别为水体反射光谱P中每条光谱的最大值、风速U的最大值、温度T的最大值，X_min分别为水体反射光谱P中每条光谱的最小值、风速U的最小值、温度T的最小值。水体反射光谱P、风速U、温度T的归一化结果分别表示为P_norm、U_norm、T_norm。

在步骤S3中，将水体反射光谱P、风速U、温度T的归一化特征P_norm、U_norm、T_norm拼接组成特征向量X和总磷浓度Y作为训练样本集输入，随机森林由100～200个决策树组成，决策树就是由不同特征经过各种组合构建成的)，数量可根据实际情况进行调整。

传统随机森林由多个决策树组成，对于随机森林{h₁(x),h₂(x),…,h_t(x)}，有输入向量X，每个决策树相对独立的对X进行结果预测，随机森林模型获得所有决策树的预测结果之后，通过设定的统计规则给出模型整体的预测结果，一般地，对于回归问题，随机森林通过计算所有决策树给出的预测结果的平均数作为最终的预测结果。本发明针对具有不同预测能力的决策树却具有相同的权重的不合理机制，通过相关系数指标来分配决策树的权重，从而提高模型整体的预测准确率。具体步骤包括：

S31、计算得到每个特征变量与总磷浓度之间的相关关系，作为特征权重的分配依据，因为相关系数取值在-1到1之间，当相关系数小于0时，说明两者呈现出负相关性，两者之间的影响仍然存在，因此在利用特征的相关系数计算特征权重的时候应该对相关系数求绝对值，利用相关系数的绝对值进行对应的特征权重的计算。具体来说，训练样本数据集X，样本总数为p，样本包含水体反射光谱波段特征和气象要素特征，样本对应总磷浓度Y，A为光谱波段特征的数量，B为气象要素特征的数量，计算每个特征的皮尔逊相关系数，得到皮尔逊相关系数集合r＝{r₁,r₂,…,r_m},m∈{1,2,…,A+B}；其中，训练样本数据集X可以通过步骤S1和S2采集、处理。

S32、使用训练样本数据集X与总磷浓度Y作为输入，采用bootstrap抽样技术从原始训练集随机产生t个训练子集，利用每个训练子集，生成对应的决策树{h₁(x),h₂(x),…,h_t(x)}；

S33、从光谱波段特征中随机选择q(q≤A)个特征与B个气象要素特征共同作为候选特征，在决策树的每个节点按照预定的规则选择最优属性进行分裂,每个决策树都得到最大限度的生长，过程中完全分裂不剪枝；

步骤S32、S33与传统的随机森林的方法一样，此处不再赘述。

S34、对于每棵决策树h_i(x)，所选特征对应的相关系数R_i＝{R_i1,R_i2,…,R_iq+B}，R_i∈r，计算权值：

上式中p_i为构建决策树h_i(x)的所有特征对应的相关系数的绝对值之和，q+B为构建该决策树使用到的特征数量，R_ij为决策树h_i(x)第j个特征所对应的相关系数；pn_i为归一化后的权值。对于每个决策树，使用的训练子集不相同，虽然均使用了q+B个特征，但是每个决策树使用的特征是不相同的，即相当于是不同的样本组合、特征组合。

S35、对于特征向量X，总磷浓度预测的结果表示为：

上式中，各决策树的预测结果

分别乘上对应的权重后相加，即获得最终的总磷浓度。

在训练时，对随机森林算法中较为重要的两个超参数：决策树棵数t和候选特征数q进行调节优化。在所有候选的参数选择中，通过循环遍历，尝试每一种可能性，将输出的总磷浓度预测值与真实值进行比较，选择出表现最好的参数作为最终的结果(调节优化属于现有技术)。

对模型性能进行评价，评价指标为决定系数R²。

式中：p代表样本数；y_i表示第i个样本总磷实测值；

表示第i个样本总磷估计值；

表示p个样本的总磷实测值的平均值。训练过程中，需要对随机森林算法的两个重要参数决策树棵数t和候选光谱特征数q进行调节优化，循环遍历两种参数的候选值，选择R²最高时的参数作为最终的结果。

模型训练建立后，在应用时，采集目标水体数据，对目标水体反射光谱以及环境气象要素作为输入，即可进行预处理和归一化，经改进的随机森林预测后，即可得到该水体的总磷浓度。

实施例1

依据本发明对某一未知总磷浓度的水体进行监测。光谱数据获取时间为北京时间上午9点-11点或下午14点-16点，波长范围350～1000nm，通道数为270，光谱分辨率优于1nm。

(1)数据采集。采集水体反射光谱曲线、环境气象要素(风速、温度)以及水样。

步骤1：采集水体反射光谱。采集水体反射光谱时需要选择晴朗无云的天气，采集时间在上午9点-11点或下午14点-16点，采集时要确保水体、太阳和设备探头之间无遮挡，选择合适的几何观测角度，每个样点重复采集10条光谱。

步骤2：采集环境气象要素。采样点处风速和温度通过手持风速仪获取，采集时需确保没有其他物体挡住风口。

步骤3：总磷浓度获取。通过塑料水瓶获取采样点处水样，水样需要及时送往实验室进行化学分析，继而获取水体的总磷浓度。

(2)数据预处理。

步骤1：对每个样点重复采集的10条光谱使用Savitzky-Golay平滑法进行平滑处理，降低水体反射光谱中的噪声，使光谱更加平滑，水体反射光谱记为P，每条水体反射光谱包含270个光谱特征。

步骤2：取步骤1获取的10条平滑后的平均曲线作为该样点的反射光谱曲线，以获取可靠的水体反射光谱。

步骤3：重复步骤2，对所有样点的最终光谱特征进行归一化处理，水体反射光谱数据集为：

其中p为样点数，n∈{1,2…270}为光谱特征数。对P进行归一化处理，环境气象要素(风速/>

温度/>

)也分别进行归一化，以消除指标之间的量纲影响。

归一化处理使用线性函数归一化方法，该方法实现对原始数据的等比例缩放,线性函数将原始数据线性化的方法转换到[0,1]的范围，归一化公式如下：

式中，X_norm为归一化后的数据，X分别为水体反射光谱P、风速U、温度T，X_min分别为水体反射光谱P中每条光谱的最小值、风速U的最小值、温度T的最小值，X_max分别为水体反射光谱P中每条光谱的最大值、风速U的最大值、温度T的最大值。水体反射光谱P的归一化结果

风速的归一化结果/>

温度的归一化结果

其中p为样点数，n为光谱特征个数。

将p个样点获取的所有归一化特征P_norm、U_norm、T_norm组成自变量

其中p为样点总数，包含n个光谱特征以及温度和风速两个环境特征，共n+2个特征。为了方便表示，设m＝n+2，则X为p行m列的二维矩阵，X_ij即表示第i个样本的第j个特征。

步骤4：将步骤3获取的所有样点的归一化结果X，与(1)中获取的总磷浓度(Y＝{y₁,y₂,…,y_p}，其中p为样点总数)，一起作为改进随机森林的输入数据。相较于神经网络的大数据特性，随机森林作为一种非线性模型，更加适合小样本建模，简单易实现，且不易陷入过拟合。

(3)将(2)中的归一化数据X和总磷浓度Y输入到改进随机森林中进行训练，利用相关系数对决策树结果进行加权以获取最终总磷预测结果并评价。

改进的随机森林的训练方法包括：

步骤1：训练样本数据集X，样本总数为p，特征数量为m，样本对应总磷浓度Y，根据公式

其中r_j为第j个特征与总磷的相关系数，/>

和/>

分别是X中的第j个特征和Y的统计数据的平均值，计算得到的皮尔逊相关系数集合r＝{r₁,r₂,…,r_m},m∈{1,2,…,n+2}。

步骤2：使用训练样本集X与总磷浓度Y作为输入，采用bootstrap抽样技术从原始训练集随机产生t个训练子集，利用每个训练子集，生成对应的决策树{h₁(x),h₂(x),…,h_t(x)}。

步骤3：从光谱特征中随机选择q(q≤n)个光谱特征，与两个环境特征一起作为候选特征，在决策树的每个节点按照预设的规则(基尼指数、信息增益率等)选择最优属性进行分裂,每个决策树都得到最大限度的生长，过程中完全分裂不剪枝。

步骤4：传统随机森林回归，一般通过计算所有决策树给出的预测结果的平均数作为最终的预测结果。改进随机森林考虑到不同决策树的不同泛化能力，结合特征的相关系数赋予每个决策树不同的权重，将所有决策树的加权和作为最终的预测结果。

对于生成的随机森林模型为{h₁(x),h₂(x),…,h_t(x)}，t为决策树个数，对于第i棵决策树h_i(x)，所选特征对应的相关系数R_i＝{R_i1,R_i2,…,R_iq+2}，R_i∈r，定义权值：

上式中p_i为构建决策树h_i(x)的所有特征对应的相关系数的绝对值之和，q+2为构建该决策树使用到的特征数量，R_ij为决策树h_i(x)第j个特征所对应的相关系数，pn_i为归一化后的权重。

随机森林模型{h₁(x),h₂(x),…,h_t(x)}拥有加权值{pn₁,pn₂,…,pn_t}，对于任意输入向量X，模型预测的结果表示为：

上式中，各决策树的预测结果

分别乘上对应的权重后相加，即获得最终的总磷浓度。

对模型性能进行评价，评价指标为决定系数R²。

式中：p代表样本数；y_i表示第i个样本总磷实测值；

表示第i个样本总磷估计值；

表示p个样本的总磷实测值的平均值。训练过程中，需要对随机森林算法的两个重要参数决策树棵数t和候选光谱特征数q进行调节优化，循环遍历两种参数的候选值，选择R²最高时的参数作为最终的结果。

(4)目标水体总磷浓度反演。按照(1)中的标准采集未知水体的相关数据，采用(2)中相同的流程对数据进行预处理，最终得到该水体的归一化数据，将归一化结果输入到(3)构建的改进随机森林模型中得到该水体的总磷浓度。

将本发明结果与不考虑环境特征的传统随机森林总磷反演模型以及考虑环境特征的传统随机森林总磷反演模型进行对比，模型精度如表1所示：

表1模型精度对比

本领域普通技术人员可以理解：以上各实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述各实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明权利要求所限定的范围。

Claims (6)

1.一种水体总磷浓度的反演方法，其特征在于，所述方法包括以下步骤：

S1、采集测试点的水体反射光谱、气象要素数据；所述气象要素数据包括温度和风速；

S2、对所述水体反射光谱和气象要素数据进行预处理；

S3、将所述预处理后的数据输入改进的随机森林，获得测试点的水体总磷浓度；

所述步骤S2中，所述预处理的方法包括：

S21、对同一地点采集的水体反射光谱构成的曲线进行平滑降噪，并获取A个光谱波段特征；

S22、对所述水体反射光谱和气象要素数据进行归一化，公式如下：

式中，X_norm为归一化后的数值，X为水体反射光谱和气象要素数据，min表示最小值，max表示最大值，将水体反射光谱的归一化结果和气象要素数据的归一化结果拼接后构成所述同一地点的特征向量；

所述步骤S3中，改进的随机森林，利用相关系数对决策树结果进行加权以获取最终总磷预测结果并评价，所述改进的随机森林的训练方法包括：

S31、训练样本数据集X，样本总数为p，样本包含光谱波段特征和气象要素特征，样本对应总磷浓度Y，A为光谱波段特征的数量，B为气象要素特征的数量，计算每个特征的皮尔逊相关系数，得到皮尔逊相关系数集合r＝{r₁,r₂,…,r_m},m∈{1,2,…,A+B}；

S32、使用训练样本数据集X与总磷浓度Y作为输入，采用bootstrap抽样技术从原始训练集随机产生t个训练子集，利用每个训练子集，生成对应的决策树{h₁(x),h₂(x),…,h_t(x)}；

S33、从光谱波段特征中随机选择q个特征与B个气象要素特征共同作为候选特征，其中，q≤A，在决策树的每个节点按照预定的规则选择最优属性进行分裂,每个决策树都得到最大限度的生长，过程中完全分裂不剪枝；

S34、计算每棵决策树h_i(x)的权值p_i，计算公式为：

其中，q+B为构建该决策树使用到的特征数量，R_ij为决策树h_i(x)的第j个特征所对应的相关系数；

然后对p_i进行归一化，公式为：

S35、对于特征向量X，总磷浓度预测的结果表示为：

其中，

为决策树h_i(x)的预测结果。

2.根据权利要求1所述的反演方法，其特征在于，所述步骤S2中，对同一地点采集多次数据，使用平均值作为最终数据。

3.根据权利要求1所述的反演方法，其特征在于，通过Savitzky-Golay平滑法进行所述平滑降噪。

4.根据权利要求1所述的反演方法，其特征在于，计算每个特征的皮尔逊相关系数的方法为：

训练样本数据集X，样本总数为p，特征数量为m，样本对应总磷浓度Y，根据公式

其中r_j为第j个特征与总磷的相关系数，/>

和/>

分别是X中的第j个特征和Y的统计数据的平均值。

5.根据权利要求1所述的反演方法，其特征在于，对反演模型性能进行评价的指标为决定系数R²，选择R²最高时的参数作为最终的结果，R²的公式为：

式中：p代表样本数；y_i表示第i个样本总磷实测值；

表示第i个样本总磷估计值；/>

表示p个样本的总磷实测值的平均值。

6.根据权利要求1所述的反演方法，其特征在于，所述训练样本数据集通过步骤S1、S2构建。

Images