函数型数据方法在水质分析中的应用

The Application of the Functional Data Method in Water Quality Analysis

摘要: 研究水质变化趋势是水质监测的重要内容。水质变化过程是一个连续的过程，只是我们监测到的数据是离散的。由于水质监测数据具有不等时间观测、非线性变化的特点以及其数据内部表现出的函数性特征，考虑采用函数型数据分析方法进行研究。在本文中，我们在对样本数据进行函数化处理的基础上，本文将函数型回归模型应用于松花江肇源段的水质分析中，预测效果良好，为该地区的水质监测提供参考。

Abstract: Studying the variation trend of water quality is an important part of water quality monitoring. Though the data we got via monitoring is discrete, the variation process of water quality is continuous. Considering the monitoring data of water quality has the characteristics of unequal time observation, nonlinear change and functional feature, we selected functional data analysis. Based on the functional processing of sample data, we used functional multiple regression method to predict water quality of Zhaoyuan section of Songhua River. And the cluster was carried out according to principal components scores. The results show that the functional data analysis method is effective. This method provides a reference for water quality monitoring in this area.

1. 引言

大多数水质变化过程是一个连续的过程，这个过程生成的统计数据可以用一个类似函数的特征表达式来描述，而通过现有的统计手段所获取的信息往往是一个不连续的、片段的、离散的有界、有序的数列。函数型数据分析方法则能较好地处理这一类的数据：它将观测数据的产生当作一个函数过程，认为样本数据之中存在着某种函数型特征，采用连续函数的方法将原本的离散数据有效地联系起来，更好地探究了数据本身的非线性变化趋势 [1]；针对函数数据的研究的目的与其他传统的统计学一样：在统计学思想和分析的指导下阐述问题；研究能够凸显数据重要特征的表现方式；为观测得到的数据建立统计模型等等 [2] [3] [4] [5]。水质预测是一个经典问题，使用的方法很多，回归分析是其中最常用的方法之一。但这些方法都是基于离散数据的方法，忽略了数据背后隐藏的连续性，采用函数型数据分析可以挖掘数据的隐藏信息。实际上，函数数据分析中的各个方法比如回归分析、聚类分析等在水质数据分析、空气质量数据分析等领域得到较为广泛地应用 [6] [7] [8] [9]。

本文将应用函数型数据分析方法进行黑龙江肇源水质数据的预测问题，实现更有效地监测水质的目的。

2. 函数型数据及函数型多元回归模型

2.1. 函数型数据

函数型数据(Function Data)是指一个集合，该集合中的元素均为定义在某个连续区间上的函数。

$F\left(x_1,x_2,\cdots ,x_p\right)=\left(\begin{array}{c}\begin{array}{c}{f}_1\left(x_1,x_2,\cdots ,x_p\right)\\ f_2\left(x_1,x_2,\cdots ,x_p\right)\\ ⋮\end{array}\\ f_n\left(x_1,x_2,\cdots ,x_p\right)\end{array}\right),$ (1)

$x_j\in \left(-\infty ,+\infty \right),j=1,2,\cdots ,p$

函数型数据是连续的数据，而一般的数据采集都是离散的。因为通过观察得到的原始数据通常以表格形式存储，这可以理解为对函数数据的截取，是自变量取一些特定的值所对应的函数值，他们是函数离散化的记录形式。显然，Ramsay所提出的函数型数据的定义域是整个区间，所以函数型数据可以包含的信息，比常见的以数据表形式出现的离散数据所包含的信息更多。

从离散的观测数据中提取连续的函数数据，我们可以用基函数法、小波变换、核函数等 [10]。本文中，我们介绍基函数法，其核心是用离散的观测值来估计其函数模型，可用下列公式表示

$\hat{x}\left(t\right)={\displaystyle {\sum }_{k=1}^K{c}_k{\phi }_k\left(t\right)}$ (2)

其中基函数 ${\phi }_k\left(t\right)\left(k=1,2,3,\cdots ,K\right)$ 的选择和系数向量 $c={\left(c_1,\cdots ,c_k\right)}^{\prime }$ 的估计是两个难点。我们可以选择傅里叶变换和样条插值作为基函数，前者针对周期性数据，后者针对非周期性数据，系数向量则通过最小二乘法来解决。

2.2. 函数型多元回归模型

多元线性回归分析是一个被广泛应用的重要方法。在自然科学和社会科学的诸多领域都有很好的应用 [11] [12]。根据回归模型的不同，多元回归分析又可以分为线性回归和非线性回归两个大类。其中线性回归模型最为成熟，其应用也是最为广泛；此外，许多非线性回归模型也可以转化为线性回归模型来求解。

$Y\left(t\right)={\beta }_1{X}_1\left(t\right)+{\beta }_2{X}_2\left(t\right)+\cdots +{\beta }_p{X}_p\left(t\right)+\epsilon \left(t\right)$ (3)

在函数型数据分析领域内，与普通离散数据的多元线性回归分析技术相对应的就是函数型数据的常系数多元线性回归问题。与普通的多元线性回归分析的不同之处在于其因变量Y和自变量 $X_1,\cdots ,X_p$，以及随机误差项 $\epsilon$ 都是以函数曲线形式存在的函数数据。而函数数据线性回归分析的目的就是通过研究因变量Y的曲线形态与自变量 $X_1,X_2,\cdots ,X_p$ 的曲线形态之间的关系，建立因变量曲线与自变量曲线的线性回归模型 [13]。

函数型数据的常系数多元回归模型的建立步骤如下：

首先，我们给出点积的定义。对于普通的离散数据而言，两个p维变量x和y的点积定义为： $\langle x,y\rangle ={\displaystyle {\sum }_{i=1}^p{x}_i\cdot y_i}$。显然，因为函数型数据的特殊性，这样的定义形式不适合定义数据型数据的点积，特此引用以下积分形式定义函数数据的点积，如下：

在函数数据空间中，对于 $x\left(t\right)\in L^2\left[a,b\right]$， $y\left(t\right)\in L^2\left[a,b\right]$， $x\left(t\right)$ 和 $y\left(t\right)$ 的点积定义为：

$\langle x\left(t\right),y\left(t\right)\rangle ={\displaystyle {\int }_a^{b}x\left(t\right)y\left(t\right)\text{d}t}$ (4)

简记为：

$\langle x\left(t\right),y\left(t\right)\rangle ={\displaystyle \int x\left(t\right)y\left(t\right)\text{d}t}$ (5)

另外，称 $I\left(t\right)=1,\forall t\in \left[a,b\right]$ 为单位函数。

假设因变量 $y\left(t\right)$ 和p个自变量 $x_1\left(t\right),x_2\left(t\right),\cdots ,x_p\left(t\right)$ 均属于 $L^2\left[a,b\right]$。若观测了n个样本，则多元线性回归总体模型可以写成：

$y_i\left(t\right)={\beta }_{0}I\left(t\right)+{\beta }_1{x}_{i1}\left(t\right)+{\beta }_2{x}_{i2}\left(t\right)+\cdots +{\beta }_p{x}_{ip}\left(t\right)+{\epsilon }_i\left(t\right),i=1,2,\cdots ,n$ (6)

其中， ${\epsilon }_i\left(t\right)$ 为随机误差项，并且 $\epsilon \left(t\right)~N\left(0,{\sigma }^2\right),\forall t\in [a,b]$ ； ${\beta }_0,{\beta }_1,\cdots ,{\beta }_p$ 是模型的待估参数，记其估计值为 ${\hat{\beta }}_0,{\hat{\beta }}_1,\cdots ,{\hat{\beta }}_p$，则多元线性回归模型为：

${\hat{y}}_i\left(t\right)={\hat{\beta }}_{0}I\left(t\right)+{\hat{\beta }}_1{x}_{i1}\left(t\right)+{\hat{\beta }}_2{x}_{i2}\left(t\right)+\cdots +{\hat{\beta }}_p{x}_{ip}\left(t\right)$ (7)

根据前面函数数据的点积定义，多元线性回归模型的残差平方和SSE可以由下列公式表示：

$\text{SSE}={\displaystyle {\sum }_{i=1}^n{\Vert y_i\left(t\right)-{\hat{\beta }}_{0}I\left(t\right)-{\displaystyle {\sum }_{j=1}^p{\hat{\beta }}_j{x}_{ij}\left(t\right)}\Vert }^2}={\displaystyle {\sum }_{i=1}^n{\displaystyle \int {\left[y_i\left(t\right)-{\hat{\beta }}_{0}I\left(t\right)-{\displaystyle {\sum }_{j=1}^p{\hat{\beta }}_j{x}_{ij}\left(t\right)}\right]}^2\text{d}t}}$ (8)

根据最小二乘原则，对上式求偏导，得：

$\frac{\partial \text{SSE}}{\partial {\hat{\beta }}_0}=-2{\displaystyle {\sum }_{i=1}^n{\displaystyle \int I\left(t\right)\left[y_i\left(t\right)-{\hat{\beta }}_{0}I\left(t\right)-{\displaystyle {\sum }_{j=1}^p{\hat{\beta }}_j{x}_{ij}\left(t\right)}\right]\text{d}t}}=0$ (9)

$\frac{\partial \text{SSE}}{\partial {\hat{\beta }}_k}=-2{\displaystyle {\sum }_{i=1}^n{\displaystyle \int x_{ik}\left(t\right)}}\left[y_i\left(t\right)-{\hat{\beta }}_{0}I\left(t\right)-{\displaystyle {\sum }_{j=1}^p{\hat{\beta }}_j{x}_{ij}\left(t\right)}\right]\text{d}t=0,k=1,2,\cdots ,p$ (10)

整理后得到正则方程如下：

${\hat{\beta }}_0{\displaystyle {\sum }_{i=1}^n{\displaystyle \int I^2\left(t\right)\text{d}t}}+{\displaystyle {\sum }_{j=1}^p{\hat{\beta }}_j}{\displaystyle {\sum }_{i=1}^n{\displaystyle \int I\left(t\right)x_{ij}\left(t\right)\text{d}t}}={\displaystyle {\sum }_{i=1}^n{\displaystyle \int I\left(t\right)y_i\left(t\right)\text{d}t}}$ (11)

${\hat{\beta }}_0{\displaystyle {\sum }_{i=1}^n{\displaystyle \int x_{ik}\left(t\right)I\left(t\right)\text{d}t}}+{\displaystyle {\sum }_{j=1}^p{\hat{\beta }}_j}{\displaystyle {\sum }_{i=1}^n{\displaystyle \int x_{ik}\left(t\right)x_{ij}\left(t\right)\text{d}t}}={\displaystyle {\sum }_{i=1}^n{\displaystyle \int x_{ik}\left(t\right)y_i\left(t\right)\text{d}t}},k=1,2,\cdots ,p$ (12)

用矩阵表示上述方程，有，

$\left[\begin{array}{cccc}{\displaystyle {\sum }_{i=1}^n{\displaystyle \int I^2\left(t\right)\text{d}t}}& {\displaystyle {\sum }_{i=1}^n{\displaystyle \int x_{i1}\left(t\right)I\left(t\right)\text{d}t}}& \cdots & {\displaystyle {\sum }_{i=1}^n{\displaystyle \int x_{ip}\left(t\right)I\left(t\right)\text{d}t}}\\ {\displaystyle {\sum }_{i=1}^n{\displaystyle \int x_{i1}\left(t\right)I\left(t\right)\text{d}t}}& {\displaystyle {\sum }_{i=1}^n{\displaystyle \int x_{i1}^2\left(t\right)\text{d}t}}& \cdots & {\displaystyle {\sum }_{i=1}^n{\displaystyle \int x_{ip}\left(t\right)x_{i1}\left(t\right)\text{d}t}}\\ ⋮& ⋮& & ⋮\\ {\displaystyle {\sum }_{i=1}^n{\displaystyle \int x_{ip}\left(t\right)I\left(t\right)\text{d}t}}& {\displaystyle {\sum }_{i=1}^n{\displaystyle \int x_{ip}\left(t\right)x_{i1}\left(t\right)\text{d}t}}& \cdots & {\displaystyle {\sum }_{i=1}^n{\displaystyle \int x_{ip}^2\left(t\right)\text{d}t}}\end{array}\right]\left[\begin{array}{c}{\hat{\beta }}_0\\ {\hat{\beta }}_1\\ ⋮\\ {\hat{\beta }}_p\end{array}\right]=\left[\begin{array}{c}{\displaystyle {\sum }_{i=1}^n{\displaystyle \int I\left(t\right)y_i\left(t\right)\text{d}t}}\\ {\displaystyle {\sum }_{i=1}^n{\displaystyle \int x_{i1}\left(t\right)y_i\left(t\right)\text{d}t}}\\ ⋮\\ {\displaystyle {\sum }_{i=1}^n{\displaystyle \int x_{ip}\left(t\right)y_i\left(t\right)\text{d}t}}\end{array}\right]$ (13)

这是一个典型的线性方程组问题，可以很方便地采用经典的高斯消元法来求解，从而得到回归系数 $\hat{\beta }={\left[\begin{array}{cccc}{\hat{\beta }}_0& {\hat{\beta }}_1& \cdots & {\hat{\beta }}_p\end{array}\right]}^{\prime }$。

回归模型即为：

$\hat{y}\left(t\right)={\hat{\beta }}^{\prime }\cdot \left(\begin{array}{c}I\left(t\right)\\ x_1\left(t\right)\\ x_2\left(t\right)\\ ⋮\\ x_p\left(t\right)\end{array}\right)$ (14)

特别的，若函数数据 $x_i\left(t\right)$ 和 $y\left(t\right)$ 是由标准正交基 $\phi$ 展开得到的，那我们的计算将会得到很大程度的简化，因为：

$\begin{array}{c}\langle x\left(t\right),y\left(t\right)\rangle ={\displaystyle \int x\left(t\right)y\left(t\right)\text{d}t}\\ ={{\displaystyle \int }}^{\text{}}\left(a_1{\phi }_q+a_2{\phi }_2+\dots +a_p{\phi }_p\right)\cdot \left(b_1{\phi }_1+b_2{\phi }_2+\cdots +b_p{\phi }_p\right)\text{d}t\\ ={{\displaystyle \int }}^{\text{}}\left(a_1{b}_1{\phi }_1^2+a_2{b}_2{\phi }_2^2+\cdots +a_p{b}_p{\phi }_p^2+a_1{b}_2{\phi }_1{\phi }_2+\cdots +a_p{b}_{p-1}{\phi }_p{\phi }_{p-1}\right)\text{d}t\\ ={\displaystyle {\sum }_{i=1}^p{a}_i{b}_i}\end{array}$ (15)

3. 案例研究

3.1. 研究区域及数据

松花江肇源江段位于黑龙江省西南部的松嫩平原第二松花江与嫩江汇合口以下，是松花江哈尔滨江段的门户。松花江肇源江段的地理位置使水环境对其影响很大。松花江上游的吉林化工企业及嫩江流域齐齐哈尔等工业城市的废水排放量很大，大庆油田的废水经古恰闸门也排入松花江，同时，吉林油田，大庆油田沿江在泛洪区和江中岛上的采油作业也对肇源地区的水环境产生很大影响，肇源江段水环境问题是松花江重要的环境问题之一。并在松花江污染防治中占有重要作用。水质数据选取4个变量，y——PH，x₁——溶解氧(mg/L)，x₂——化学需氧量(mg/L)，x₃——氨氮(mg/L)，原始数据见表1。

我们选取黑龙江肇源2006~2011年水质监测的数据，如下表所示：

Table 1. Water quality monitoring data of Zhaoyuan section of Songhua River from 2006 to 2011

表1. 松花江肇源段 2006~2011年水质监测数据

3.2. 水质预测模型

3.2.1. 模型建立

对这4个变量进行8项傅里叶级数曲线拟合。然后，按照函数型数据多元回归建模方法实现下列线性方程组的计算：

$\left[\begin{array}{cccc}{\displaystyle {\sum }_{i=1}^n{\displaystyle \int I^2\left(t\right)\text{d}t}}& {\displaystyle {\sum }_{i=1}^n{\displaystyle \int x_{i1}\left(t\right)I\left(t\right)\text{d}t}}& \cdots & {\displaystyle {\sum }_{i=1}^n{\displaystyle \int x_{ip}\left(t\right)I\left(t\right)\text{d}t}}\\ {\displaystyle {\sum }_{i=1}^n{\displaystyle \int x_{i1}\left(t\right)I\left(t\right)\text{d}t}}& {\displaystyle {\sum }_{i=1}^n{\displaystyle \int x_{i1}^2\left(t\right)\text{d}t}}& \cdots & {\displaystyle {\sum }_{i=1}^n{\displaystyle \int x_{ip}\left(t\right)x_{i1}\left(t\right)\text{d}t}}\\ ⋮& ⋮& \ddots & ⋮\\ {\displaystyle {\sum }_{i=1}^n{\displaystyle \int x_{ip}\left(t\right)I\left(t\right)\text{d}t}}& {\displaystyle {\sum }_{i=1}^n{\displaystyle \int x_{ip}\left(t\right)x_{i1}\left(t\right)\text{d}t}}& \cdots & {\displaystyle {\sum }_{i=1}^n{\displaystyle \int x_{ip}^2\left(t\right)\text{d}t}}\end{array}\right]\left[\begin{array}{c}{\hat{\beta }}_0\\ {\hat{\beta }}_1\\ ⋮\\ {\hat{\beta }}_p\end{array}\right]=\left[\begin{array}{c}{\displaystyle {\sum }_{i=1}^n{\displaystyle \int I\left(t\right)y_i\left(t\right)\text{d}t}}\\ {\displaystyle {\sum }_{i=1}^n{\displaystyle \int x_{i1}\left(t\right)y_i\left(t\right)\text{d}t}}\\ ⋮\\ {\displaystyle {\sum }_{i=1}^n{\displaystyle \int x_{ip}\left(t\right)y_i\left(t\right)\text{d}t}}\end{array}\right]$ (16)

得到的结果是：

$\left[\begin{array}{cccc}271& 2054& 1645& 170\\ 2054& 16014& 12217& 1360\\ 1645& 12217& 10446& 983\\ 170& 1360& 983& 137\end{array}\right]\left[\begin{array}{c}{\hat{\beta }}_0\\ {\hat{\beta }}_1\\ {\hat{\beta }}_2\\ {\hat{\beta }}_3\\ {\hat{\beta }}_4\end{array}\right]=\left[\begin{array}{c}2020\\ 15324\\ 12238\\ 1276\end{array}\right]$

计算得到：

$\hat{\beta }={\left(\begin{array}{cccc}7.8367& -0.0454& -0.0407& 0.3327\end{array}\right)}^{\prime }$

所以最终得到的回归模型是：

$y=7.8367-0.0454x_1-0.0407x_2+0.3227x_3$

3.2.2. 模型拟合优度

拟合优度 $R^2$ 是检验回归模型优劣的重要指标，是检验模型拟合实际数据的接近程度。该模型的拟合优度为( $R^2$ 越接近1，效果越好)：

$R^2=\frac{\text{SSR}}{\text{SST}}=0.9089$

从理论上讲，模型考虑了数据的连续性，能挖掘数据的隐藏信息，实际计算的拟合优度 $R^2>0.9$，说明模型的效果优良，可以放入更大型的数据中进行水质数据的预测工作，为水质监测工作提供参考。

4. 结论

水质数据本质上是连续性数据，本文利用函数型数据中的回归方法进行黑龙江省松花江肇源段水质数据的预测工作。该方法具有理论优势，模型精确度也较高。该方法的引入有利于对水质监测工作进行科学合理地改进。