分段回归突变点确定浅水湖泊水质管理目标方法与流程

本发明涉及分段回归突变点确定浅水湖泊水质管理目标方法，属于湖泊领域。

背景技术：

我国大部分淡水湖泊都是浅水湖泊且富营养化问题较为严重，水质趋于恶化；而氮和磷是湖泊水体中最主要的营养盐，控制湖泊特别是浅水湖泊中氮和磷的含量对于降低湖泊富营养化程度，保护湖泊水体质量具有决定性的作用。制定恰当的氮磷营养盐水质管理目标是有效保护浅水湖泊水环境质量的前提条件，对于保护浅水湖泊水环境具有十分重要的意义。目前，浅水湖泊营养盐水质目标的制定多数依靠专家和领导的经验作出决策，其人为主观性较强；这导致制定的湖泊营养盐水质管理目标时常脱离湖泊水质的实际情况。事实上，若浅水湖泊水环境管理目标过于宽松则不利于湖泊水环境保护，导致浅水湖泊进一步富营养化，水质更加恶化；若管理目标过于严格则既无必要又可能妨碍当地经济社会发展。

技术实现要素：

发明目的：为了克服现有技术中存在的不足，本发明提供一种分段回归突变点确定浅水湖泊水质管理目标方法，完全利用了客观观测数据中自然展现的统计突变点，精确的确定水质等级区间，易于在各种受污染浅水湖泊水质管理目标的确定中推广使用。

技术方案：为解决上述技术问题，本发明的一种分段回归突变点确定浅水湖泊水质管理目标方法，包括以下步骤：

(1)搜集需要确定水质管理目标的浅水湖泊n组叶绿素a的浓度、总磷浓度以及总氮浓度，n为大于50的整数，Z_i、X_i和Y_i分别表示第i组叶绿素a的浓度、总磷浓度和总氮浓度的值，统一叶绿素a的浓度，总磷浓度和总氮浓度的单位；并将叶绿素a的浓度、总磷浓度以及总氮浓度分别取自然对数e为底的对数，获得对应的值logZ_i、logX_i及logY_i；

(2)以叶绿素a的浓度对数值为因变量，以总磷浓度对数值为自变量做非参数回归，并初步确定总磷突变点；

(3)以叶绿素a浓度对数值为因变量，以总磷浓度对数值为自变量做分段线性回归，用步骤(2)中确定的总磷突变点作为求解初始值，精确确定总磷突变点；

(4)用总氮浓度代替总磷浓度，重复步骤(2)和(3)精确确定总氮的突变点；

(5)根据总氮突变点和总磷突变点及湖泊所在区域的社会经济发展情况和国家水质标准，确定浅水湖泊候选营养盐水质管理目标；

(6)确定浅水湖泊水体用途所要求水质并与候选营养盐水质管理目标比较，最终确定具体湖泊营养盐水质管理目标。

作为优选，所述步骤(2)具体包括以下步骤：

a.根据湖泊总磷浓度和叶绿素a的浓度的观测值数据数量n，用h＝n^-1/5计算最佳带宽值h；

b.在总磷浓度对数值的取值范围内均匀取1000个数值作为自变量logM的值logM_p，即logM₁、logM₂、logM₃...logM₁₀₀₀；

c.从logM₁开始计算得到所有满足条件logX_i≤logM₁+3h及logX_i≥logM₁-3h的总磷对数值logX_i，对每一个满足条件的logX_i及logM₁计算相应的核函数K(logX_i-logM₁)，作为logX_i在点logM₁上的权重，不满足条件logX_i≤logM₁+3h及logX_i≥logM₁-3h的logX_i权重为0；其中：

$K(\log X_i-\log M_1)=\frac{1}{\sqrt{2\mathrm{\π}}h}{e}^{-\frac{1}{2h^2}{({\mathrm{logX}}_i-{\mathrm{logM}}_1)}^2}$

d.利用已有数值求解法，对logM₁依次求解两组实数a(logM₁)，b(logM₁)；具体即要求a(logM₁)，b(logM₁)使下式最小，

$\underset{i=1}{\overset{n}{\Σ}}{\{\log Z_i-a\left(\log M_1\right)-b\left(\log M_1\right)\log X_i\}}^{2}K(\log X_i-\log M_1)$

e.重复步骤c和d，依次将logM₂、logM₃...logM₁₀₀₀替换logM₁计算出两组实数a(logM_p)，b(logM_p)；

f.画出回归方程logZ＝a(logM_p)+b(logM_p)logM的曲线，其中logZ代表回归方程的因变量，logM是回归方程的自变量；再次注意a(logM_p)，b(logM_p)是回归方程的系数，是随着自变量logM的取值logM_p的不同而不同的两组实数，观察方程logZ＝a(logM_p)+b(logM_p)logM的曲线中的转折点，并记录该变量的突变点logγ对应的值作为分段回归求解迭代的初始值，如有多个总磷突变点，取最小的总磷突变点。

作为优选，所述步骤(3)的求解方法是求解四个实数β₀，β₁，δ和使下式最小，求解方法可见参考文献(Ulm,1991)：

其中是示性函数，也即满足下式：

上面方程可以用已有标准非线性回归求解程序迭代计算，即为所需求解的总磷精确的突变点，迭代计算初始值的好坏直接影响求解的效率和结果的准确性；将logγ代入方程作为的迭代初始值能很好的提高计算效率和结果的准确性。的初始值为logγ。

作为优选，所述步骤(5)步骤具体包括以下步骤：

a.根据湖泊所在区域的经济发展指标，确定该区域经济发展等级，我国各地区经济发展等级包括：

A.发达地区：区域人均国民生产总值超过全国人均国民生产总值的2倍，即区域人均国民生产总值大大超过全国人均国民生产总值，人民生活富足；

B.较发达地区：区域人均国民生产总值是全国人均国民生产总值的1.3倍至2倍，即区域人均国民生产总值超过全国人均国民生产总值，人民生活比较富足；

C.中等地区：区域人均国民生产总值是全国人均国民生产总值的0.8至1.3倍，即区域人均国民生产总值与全国人均国民生产总值差不多，人民生活水平中等；

D.落后地区：区域人均国民生产总值低于全国人均国民生产总值0.8倍，即区域人均国民生产总值大大低于全国人均国民生产总值，人民生活较贫困。

b.将总氮、总磷和叶绿素a的突变点值与国家《地表水环境质量标准(GB3838-2002)》比较，确定突变点所在水质等级区间；

c.发达地区应制定更严格的浅水湖泊营养盐管理目标，因此在突变点值的基础上，升高一个等级作为营养盐水质管理目标的候选；较发达地区取突变点所在水质等级区间对应的高等级作为营养盐水质管理目标的候选；中等地区取突变点所在水质等级区间对应的低等级作为营养盐水质管理目标的候选；落后地区在突变点基础上，降低一个等级作为营养盐水质管理目标的候选。

作为优选，所述步骤(6)步骤具体包括以下步骤：

a.确定湖泊水体用途及所要求水质：按照国家水质标准水体用途可分为自然保护区、饮用水源地、渔业、工业和农业五大类；依次对应于《地表水环境质量标准(GB3838-2002)》中的I、II、III、IV和V类水；

b.将候选营养盐水质管理目标与该湖泊水体用途所要求水质进行比较，最终取两者中水质较为严格的结果作为具体湖泊营养盐水质管理的目标。

在本发明中，环境学的研究指出，浅水湖泊富营养化的主要指标叶绿素a浓度与营养盐总氮、总磷浓度之间存在着密切关系；且这些环境因子彼此关系中存在突变点。本发明依据以上观点，利用分段线性回归方法精确确定总氮、总磷和叶绿素a关系中的突变点；结合国家地表水环境质量标准、浅水湖泊所在区域的社会经济发展情况及湖泊水体的具体用途客观精确确定浅水湖泊营养盐水质管理目标的方法。本发明结合湖泊当地经济社会发展水平和湖泊水体具体用途，客观精确确定浅水湖泊水质管理目标的方法。该方法能帮助环境管理部门决定正确客观的管理目标，达到既能充分利用湖泊水资源，为经济社会发展提供必要的条件又能有效保护和改善湖泊水质，降低浅水湖泊富营养化程度。

有益效果：本发明的分段回归突变点确定浅水湖泊水质管理目标方法，完全利用了客观观测数据中自然展现的统计学突变点和国家权威的标准，几乎没有人为主观干预的步骤，克服了已有方法主观性和人为干预较强的缺点，可精确的确定湖泊区域的水质标准，且综合考虑了经济社会发展，湖泊水体用途因素，易于在各种受污染浅水湖泊水质管理目标的确定中推广使用。

附图说明

图1为本发明的流程图；

图2为总磷－叶绿素a局部线性回归图；

图3为总磷－叶绿素a分段线性回归图；

图4为总氮－叶绿素a局部线性回归图；

图5为总氮－叶绿素a分段线性回归图。

具体实施方式

下面结合附图和我国东部某浅水湖泊实际观测数据对本发明作更进一步的说明。

(1)按照图1所示流程图，首先从文献及实测中获得该浅水湖泊总磷、总氮和叶绿素a的浓度值共151组；总氮和总磷浓度的单位均为毫克/升，叶绿素a浓度单位一般为微克/升，将其除以1000后统一与总氮、总磷浓度单位一致的数值。

(2)将总磷、总氮和叶绿素a的浓度取对数值，同时将n＝151，代入h＝n^-1/5得到理论最佳带宽h＝0.366，本次计算中取h＝0.37。

设定在总磷对数值的取值范围内要计算的自变量点的个数，为了最后显示的局部回归线较为光滑，本次在总磷对数值的取值范围-4.6至-0.5内均匀取1000个数值作为logMp，即为logM₁、logM₂、logM₃…logM₁₀₀₀，对于每一个logM_p对比总磷观测值的对数logX_i，所有满足条件logX_i≤logM_p+3h及logX_i≥logM_p-3h的logX_i，分别计算核函数K(logX_i-logM_p)，其中：

$K(\log X_i-\log M_p)=\frac{1}{\sqrt{2\mathrm{\π}}h}{e}^{-\frac{1}{2h^2}{({\mathrm{logX}}_i-{\mathrm{logM}}_p)}^2}$

求解a(logM_p)，b(logM_p)使下式取值最小：

$\underset{i=1}{\overset{n}{\Σ}}{\{\log Z_i-a\left(\log M_p\right)-b\left(\log M_p\right)\log X_i\}}^{2}K(\log X_i-\log M_p)$

这个求解用已有的Newton切线法完成，对于每一个logM_p，在满足上式和最小的条件下得到一组a(logM_p)，b(logM_p)，如图2所示画出方程logZ＝a(logM_p)+b(logM_p)logM的曲线，这条曲线的横纵坐标分别代表总磷浓度的对数值和对应的叶绿素a的对数值，图上有较明显的多个突变点，出于保护湖泊环境的目的，其目标浓度应该较小，所以应该取总磷值最小的突变点，在这个例子中突变点对应的总磷值约为0.045毫克/升，叶绿素a的值约为0.0037毫克/升。

(3)将上一步求得的总磷0.045毫克/升，作为突变点迭代的初始值代入中，求解β₀，β₁，δ和使这个和最小，具体求解方法参考Ulm,K.,1991.A statistical method for assessing a threshold in epidemiological studies.Statistics in medicine 10,341-349；其中是示性函数，也即下式：

对应的分段线性回归方程为：

如图3，在本例中，最终求出各参数的值：β₀＝-6.12；β₁＝-0.32；δ＝0.50；最重要的是精确估计的总磷突变点为0.061毫克/升。

(4)如图4和图5所示，将研究对象换为总氮－叶绿素a重复步骤(2)-(3)得到总氮突变点对应浓度为0.71毫克/升。

(5)对照《地表水环境质量标准(GB3838-2002)》可以发现总磷突变点0.061毫克/升，位于III类水和IV类水标准之间；而总氮突变点0.71毫克/升位于II类水和III类水标准之间。该浅水湖泊位于江苏南部的苏州和无锡两市辖区内，该区域人均国民生产总值远高于全国人均国民生产总值(区域人均国民生产总值超过全国人均国民生产总值的2倍，人民生活富足)，人民生活富裕，是我国经济最发达的区域之一。因此，升高一个等级作为营养盐水质管理目标的候选，即总氮的候选水质管理目标定为I类水水质标准，总磷的水质管理目标定为II类水。

(6)经调查该湖泊主要的用途为农业和渔业用水，其对应的水质标准为III类水；因此，经过与候选的总氮和总磷水质管理目标比较，最终该湖泊的总氮水质管理目标是I类水水质标准，总磷的水质管理目标定为II类水水质标准。

本发明完全利用了客观观测数据中自然展现的统计学突变点和国家权威的标准，几乎没有人为主观干预的步骤，克服了已有方法主观性和人为干预较强的缺点，能够精确的确定目标区域的水质标准，便于目标区域的水质管理。该方法对所用数据无任何特殊要求，且综合考虑了经济社会发展，湖泊水体用途因素，易于在各种受污染浅水湖泊水质管理目标的确定中推广使用。

以上所述仅是本发明的优选实施方式，应当指出：对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。