一种水体水质变化检测方法及系统、存储介质和装置与流程

1.本发明涉及水质监测技术领域，具体涉及一种水体水质变化检测方法及系统、存储介质和装置。


背景技术：

2.近年来，随着我国城镇化的快速发展，工业生产、农业灌溉和人民日常用水量日趋增加。人们生产、生活产生的工业废水和生活污水被排放到城市河道中，造成城市水体水质日益恶化，严重破坏城市河流生态系统，对城市景观和城市发展造成较大影响，对区域的生态环境、人居环境和人们的身体健康造成了严重的危害。监测和治理水体污染是目前亟待解决的城市环境问题。传统水质监测方法主要采用地面水站监测，监测范围有限，难以实现长时间连续的大范围监测。遥感技术可以从全空间领域，多时间序列快速、连续地体现整个水域的水质状况。
3.水体的水质参数对于遥感影像的红外和近红外波段敏感，在此波段内污染水体与非污染水体光谱特征明显不同，可辅助识别水体的污染状况。水体中的叶绿素a浓度、悬浮物浓度和透明度是水质监测的重要指标参数。利用中高分辨率卫星遥感影像，基于水质参数反演模型对研究区的水质参数进行反演可以精确和客观的得到水体污染状况分布情况。在此基础上，利用遥感时间序列数据结合分析方法可以对水体水质的变化情况进行进一步的描述。随着人们对变化检测的时间精度要求不断提高以及输入数据的时间分辨率的不断增加，双时相、多时相的变化检测方法逐渐被淘汰，基于稠密时间序列（指一年内多景影像）的变化分析可以对水质的变化情况进行更细致地分析。
4.如何基于稠密遥感时间序列数据构建水质参数变化分析方法是亟待解决的问题。当前算法在对突变信号的检测上取得了一定成功，但是在渐变的描述上，当前算法得到的分段线性的结果过于简单，还不能与水质变化建立紧密的联系，其关键点在于水质指标的年内季节性变化干扰——季相波动干扰了对时间序列整体趋势的定量表达，影响了对指标整体趋势的判断。


技术实现要素：

5.本发明的目的在于提供一种实现稠密时间序列下的水体水质变化检测方法，解决现有水质变化检测无法排除年内季节性变化干扰的问题，导致水质指标整体趋势的判断不准确。
6.本发明通过下述技术方案实现：一种水体水质变化检测方法，包括以下步骤：s1、对获取的水质指标历史遥感数据进行预处理，获得监测水体范围内的反射率影像；s2、基于步骤s1获得的反射率影像，采用反演模型对水质指标进行反演；s3、以年积日（day of year，doy）作为自变量（不考虑年份），以经过步骤s2反演后
的水质指标作为因变量构建拟合函数来估计水质指标的长期平均年度变化，采用最小二乘法对拟合函数进行处理，获得水质指标的拟合时间序列；s4、用逐像元对所有经过步骤s2反演后的水质指标减去对应的经过步骤s3拟合的水质指标获得残差值，基于残差值构建残差时间序列；s5、基于步骤s4获得的残差时间序列，采用mk检验构建uf
k
时间序列，采用uf
k
时间序列曲线判断监测水体范围内的水质指标变化。
7.由于水质指标本身存在季节性波动，对于一年内有多期的时间序列影像，现有技术直接比较年内两个时刻的指标是不科学的。
8.本发明的目的在于：如何在比较年内多个时间点的指标值时能屏蔽掉季节性变化的干扰。
9.本发明的残差是原始时间序列曲线减掉了“季节性变化”得到的，比如图2中的残差时间序列曲线已经不存在季节性变化了，在图2和图3中，虽然以年作为刻度，但是还是以日期作为自变量的，比如1988代表的是1988年1月1日；其次，doy和真实的日期不是对立的概念，1988年1月15日就是1988年的第15天，即doy=15，求残差时间序列和后面的mk检验区别是求残差时间序列抛弃了“年”的概念，因为要考虑多年平均水平，所以就把多年的时间序列点放到了一年里。
10.即本发明通过构建残差时间序列，能够有效屏蔽季节性变化的干扰，得到准确的水质指标变化；进而能够实现面向稠密时间序列的水质指标变化监测。
11.本发明通过残差时间序列和mk时间序列检验方法与水质指标相结合以实现水质变化的探测，具有计算简单，而且可以明确突变开始的位置，并指出突变区域。
12.进一步，步骤s1中，所述预处理依次包括辐射校正、正射校正和影像镶嵌。
13.首先进行辐射校正，消除传感器成像条件、大气条件、太阳位置和角度条件及某些不可避免的噪声引起的传感器所得到的目标测量值与目标的光谱发射率或光谱辐亮度等物理量之间的差异。然后进行正射校正，即对影像空间和几何畸变进行校正，消除影像倾斜和投影误差，从而输出多中心投影平面正射影像的处理过程。最后进行影像镶嵌，将覆盖研究范围的相同数据源的遥感影像进行镶嵌，并按照研究区矢量范围对镶嵌后影像进行裁剪，得到研究区范围内的反射率影像。
14.进一步，步骤s2中，采用遥感影像的波段差异处的绿光波段和红光波段作为特征波段进行反演。
15.反演模型都是现有技术，根据不同的水质指标采用不同的反演模型进行反演，例如：选用遥感影像的波段差异处的绿光波段和红光波段作为反演水体悬浮物浓度的主要特征波段。
16.c
tss
=119.62
·
（b3/b2）
6.0823
其中，c
tss
为悬浮物浓度，单位为mg
·
l
‑1；b2和b3分别为遥感影像的绿光波段和红光波段。
17.透明度与悬浮物浓度呈负相关，悬浮物浓度越大的河段，其透明度反而越小。反演模型为：z
sd
=284.15
·
c
tss
‑
0.67
式中，c
tss
为悬浮物浓度，z
sd 为水体透明度，单位为cm。
18.水体的缺氧性和富营养为评价水体污染程度的重要指标，主要是基于水体的叶绿素a浓度和透明度进行建立综合营养指数模型：tli
sd
=51.18－19.4
·
lnz
sd
其中，tli
sd 基于透明度的营养状态指数。
19.进一步，步骤s3中，拟合时间序列的模型如下：其中，t是时间（以doy表示）, m1是背景ndvi值，m2是ndvi的最大值，m3和m4控制斜率，m5和m6控制相位。由此得到了整个时间序列的平均年内季相波动情况（图1）。
20.m1‑
m6是通过拟合得到的，这6个参数也代表拟合曲线的形状。在拟合中的初始数据有足够多的落在坐标系中的点，以及一个方程式。比如，给出一些散点做线性拟合y=ax b，得到了的结果是y=2x 1，这个2和1就是拟合得到的参数，它们也具有物理意义，比如2代表斜率，1代表截距，类似于本发明中的m1‑
m6，只不过本发明中采用的是稍微复杂一些的非线性拟合。
21.进一步，步骤s4中，残差的计算模型如下：i
r
（t）=i（t）－i
f
式中，i
f
为时间序列，i（t）为观测值，观测值指的是影像反演的水质指标。
22.进一步，步骤s5中，将残差时间序列曲线和uf
k
时间序列曲线，以及mk检验中当给定显著性水平α=0.05，u
0.05
=
±
1.96时的两条临界直线均绘在同一张图上，若uf
k
的值大于0，则表明序列呈上升趋势，小于0则表明呈下降趋势，当uf
k
的值的值超过临界直线时，表明上升或下降趋势显著。
23.其中，u
0.05
=
±
1.96分别为u
0.05
= 1.96和u
0.05
=
‑
1.96对应两个点值的临界直线。u
0.05
= 1.96和u
0.05
=
‑
1.96是通过查表获得，uf
k
时间序列曲线为根据mk检验算法对残差时间序列进行处理获得的曲线，可以理解为一个计算方法，比如mk的结果是原时间序列曲线各个点都 1，它的结果还是一条时间序列曲线，跟原曲线都是“对齐”的，其中的 uf
k
就是一个命名，以区分残差时间序列曲线，换句话说，mk序列趋势检验是对某一时间序列曲线趋势的统计学判断，在mk趋势检验的结果是一个uf
k
序列，其中某一个点的uf
k
值是原时间序列曲线从第1各点到该点趋势判断结果（比如ufk>1.96则称之为：从第一个点到当前点有显著增加的趋势），即本发明的uf
k
时间序列曲线是对残差时间序列曲线的检验结果。
24.mk检验是一种非参数统计检验，为现有技术，该方法用统计方法评估变量是否随时间存在单调增加或减少趋势。
25.在数学上，mk检验统计量s是通过从负差异的数量中减去数据中的正差异的数量来计算的，通过以下公式计算：其中x
j
和x
k
观察值分别在1，2，
……
，n处获得，如果后来收集的数据大于先前获得的数据，则s为正，反之亦然，结果表明，在零假设下，mk统计量可以通过正态分布近似，条件
是观测数超过10个。平均值和方差由下式给出：是观测数超过10个。平均值和方差由下式给出：是观测数超过10个。平均值和方差由下式给出：其中，n是时间序列中的观察总数。是调整序列相关后的有效样本量。是关于滞后i的数据等级之间的自相关函数。由此定义统计量：式中uf1=0，e(s
k
)、var(s
k
)是累计数s
k
的均值和方差。mk的结果以统计量给出，记为uf
i
，uf
i
为标准正态分布，其按时间序列x顺序x1,x2,
…
x
n
，计算出的统计量序列，给定显著性水平α，查正态分布表，若|uf
i
|>u
α
，则表明序列存在明显的趋势变化。这一方法的优点在于不仅计算简单，而且可以明确突变开始的位置，并指出突变区域。因此，是一种常用的监测方法。本研究给定显著性水平α=0.05，那么根据查表得临界值u
0.05
=
±
1.96。将uf
k
序列曲线和
±
1.96两条直线均绘在同一张图上，若uf
k
的值大于0，则表明序列呈上升趋势，小于0则表明呈下降趋势，当它们超过临界直线时，表明上升或下降趋势显著。
26.一种水体水质变化检测系统，包括：数据采集模块：用于获取水质指标历史遥感数据；预处理模块：用于对水质指标历史遥感数据进行预处理获得反射率影像；反演模块：基于反射率影像，采用存储的反演模型对水质指标进行反演；第一数据处理模块：用于对反演后的水质指标进行拟合获得拟合时间序列；第二数据处理模块：用于获取第一数据处理模块中的拟合时间序列和反演模块中的水质指标，获得残差值，基于残差值构建残差时间序列；判断模块：基于残差时间序列，采用mk检验构建uf
k
时间序列；输出模块：用于输出残差时间序列和uf
k
时间序列对应在同一坐标下的曲线图。
27.一种计算机可读存储介质，该存储介质存储有上述系统。
28.一种计算机可读存储介质，该存储介质存储有指令，该指令用于当被处理器执行时使所述处理器执行水体水质变化检测方法。
29.一种装置，包括处理器，所述处理器配置为执行水体水质变化检测方法。
30.一种装置，该装置包含计算机可读存储介质。
31.本发明与现有技术相比，具有如下的优点和有益效果：1、本发明能够有效屏蔽季节性变化的干扰，提高了水质指标整体趋势的判断的准确性。
32.2、本发明通过残差时间序列和mk时间序列检验方法与水质指标相结合以实现水质变化的探测，具有计算简单，而且可以明确突变开始的位置，并指出突变区域。
33.3、本发明能够实现稠密时间序列下水质指标变化的精细描述。
附图说明
34.此处所说明的附图用来提供对本发明实施例的进一步理解，构成本技术的一部分，并不构成对本发明实施例的限定。在附图中：图1为基于时间序列数据的双逻辑函数拟合示意图；图2为用原始指标值减拟合指标值残差时间序列示意图；图3为mk时间序列检验趋势判断结果图；图中，两条临界直线之间的曲线为残差时间序列曲线。
具体实施方式
35.为使本发明的目的、技术方案和优点更加清楚明白，下面结合实施例和附图，对本发明作进一步的详细说明，本发明的示意性实施方式及其说明仅用于解释本发明，并不作为对本发明的限定。
36.实施例1：一种水体水质变化检测方法，包括以下步骤：s1、对获取的水质指标历史遥感数据进行预处理，获得监测水体范围内的反射率影像；首先进行辐射校正，消除传感器成像条件、大气条件、太阳位置和角度条件及某些不可避免的噪声引起的传感器所得到的目标测量值与目标的光谱发射率或光谱辐亮度等物理量之间的差异。然后进行正射校正，即对影像空间和几何畸变进行校正，消除影像倾斜和投影误差，从而输出多中心投影平面正射影像的处理过程。最后进行影像镶嵌，将覆盖研究范围的相同数据源的遥感影像进行镶嵌，并按照研究区矢量范围对镶嵌后影像进行裁剪，得到研究区范围内的反射率影像；s2、基于步骤s1获得的反射率影像，采用反演模型对水质指标进行反演，具体采用遥感影像的波段差异处的绿光波段和红光波段作为特征波段进行反演；根据不同的水质指标采用不同的反演模型进行反演，例如：选用遥感影像的波段差异处的绿光波段和红光波段作为反演水体悬浮物浓度的主要特征波段。
37.c
tss
=119.62
·
（b3/b2）
6.0823
其中，c
tss
为悬浮物浓度，单位为mg
·
l
‑1；b2和b3分别为遥感影像的绿光波段和红光波段；透明度与悬浮物浓度呈负相关，悬浮物浓度越大的河段，其透明度反而越小。反演模型为：z
sd
=284.15
·
c
tss
‑
0.67 式中，c
tss
为悬浮物浓度，z
sd
为水体透明度，单位为cm。
38.水体的缺氧性和富营养为评价水体污染程度的重要指标，主要是基于水体的叶绿素a浓度和透明度进行建立综合营养指数模型：tli
sd
=51.18－19.4
·
lnz
sd
其中，tli
sd
为基于透明度的营养状态指数。
39.s3、以年积日作为自变量，以经过步骤s2反演后的水质指标作为因变量构建拟合
函数，采用最小二乘法对拟合函数进行处理，获得水质指标的拟合时间序列；拟合时间序列的模型如下：其中，t是时间（以doy表示）, m1是背景ndvi值，m2是ndvi的最大值，m3和m4控制斜率，m5和m6控制相位。由此得到了整个时间序列的平均年内季相波动情况，如图1所示；s4、用逐像元对所有经过步骤s2反演后的水质指标减去对应的经过步骤s3拟合的水质指标获得残差值，基于残差值构建残差时间序列，残差时间序列如图2所示；残差的计算模型如下：i
r
（t）=i（t）－i
f
式中，i
f
为时间序列，i（t）为观测值，观测值指的是影像反演的水质指标；s5、基于步骤s4获得的残差时间序列，采用mk检验构建uf
k
时间序列，采用uf
k
时间序列曲线判断监测水体范围内的水质指标变化；mk检验是一种非参数统计检验，该方法用统计方法评估变量是否随时间存在单调增加或减少趋势。
40.在数学上，mk检验统计量s是通过从负差异的数量中减去数据中的正差异的数量来计算的，通过以下公式计算：其中x
j
和x
k
观察值分别在1，2，
……
，n处获得，如果后来收集的数据大于先前获得的数据，则s为正，反之亦然，结果表明，在零假设下，mk统计量可以通过正态分布近似，条件是观测数超过10个。平均值和方差由下式给出：是观测数超过10个。平均值和方差由下式给出：是观测数超过10个。平均值和方差由下式给出：其中，n是时间序列中的观察总数。是调整序列相关后的有效样本量。是关于滞后i的数据等级之间的自相关函数。由此定义统计量：式中uf1=0，e(s
k
)、var(s
k
)是累计数s
k
的均值和方差。mk的结果以统计量给出，记为uf
i
，uf
i
为标准正态分布，其按时间序列x顺序x1,x2,
…
x
n
，计算出的统计量序列，给定显著性水平α，查正态分布表，若|uf
i
|>u
α
，则表明序列存在明显的趋势变化。这一方法的优点在于不仅计算简单，而且可以明确突变开始的位置，并指出突变区域。因此，是一种常用的监测方法。本研究给定显著性水平α=0.05，那么根据查表得临界值u
0.05
=
±
1.96。将uf
k
序列曲线和
±
1.96两条直线均绘在同一张图上，若uf
k
的值大于0，则表明序列呈上升趋势，小于0
则表明呈下降趋势，当它们超过临界直线时，表明上升或下降趋势显著。
41.如图3所示，两条临界直线之间的曲线为残差时间序列曲线，另外一条线为uf
k
序列曲线，uf是mk统计量，它也是一个时间序列，换言之每一个原始时间序列点对应一个残差时间序列点，也对应uf时间序列的一个点，假设我们想看看从1988年某月某日到2013年某月某日这个指标有没有表现出显著的趋势变化呢，我们可以看图3中2013年该日的mk统计量大于1.96，那么我们可以给出结论“从1988年某月某日到2013年某月某日水质指标表现出了显著的上升趋势”，而类似地，可以得出结论“从1988年某月某日到1995年某月某日水质指标没有表现出了显著的上升趋势”。当然，图3给出的是1988
‑
2017年这两种情况下，2013或1995年之后的数据就不用统计了。
42.实施例2：一种水体水质变化检测系统，包括：数据采集模块：用于获取水质指标历史遥感数据；预处理模块：用于对水质指标历史遥感数据进行预处理获得反射率影像；反演模块：基于反射率影像，采用存储的反演模型对水质指标进行反演；第一数据处理模块：用于对反演后的水质指标进行拟合获得拟合时间序列；第二数据处理模块：用于获取第一数据处理模块中的拟合时间序列和反演模块中的水质指标，获得残差值，基于残差值构建残差时间序列；判断模块：基于残差时间序列，采用mk检验构建uf
k
时间序列；输出模块：用于输出残差时间序列和uf
k
时间序列对应在同一坐标下的曲线图。
43.实施例3：一种计算机可读存储介质，该存储介质存储有实施例2所述系统。
44.实施例4：一种计算机可读存储介质，该存储介质存储有指令，该指令用于当被处理器执行时使所述处理器执行实施例1所述的水体水质变化检测方法。
45.实施例5：一种装置，包括处理器，所述处理器配置为执行实施例1所述的水体水质变化检测方法。
46.实施例6：一种装置，该装置包含权实施例3或实施例4所述计算机可读存储介质。
47.以上所述的具体实施方式，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施方式而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。