CN103913422A - 一种基于hj-ccd影像的浅水湖泊水生植物快速监测方法 - Google Patents

Abstract

本发明提供了一种基于资源环境一号小卫星（HJ-CCD）影像的浅水湖泊水生植物快速监测方法，其包括以下步骤：1）获取研究区的HJ-CCD影像，并进行预处理；2）根据研究区HJ-CCD影像，计算挺水植物敏感特征指数；3）对研究区HJ-CCD影像进行主成分变换，获取浮叶植物敏感特征；4）对研究区HJ-CCD影像进行缨帽变换，获取沉水植物敏感特征指数；5）根据研究区影像同步的实测数据，确定、和识别不同水生植物的阈值，提出判别条件；6）根据判别条件建立决策树分类模型，生成浅水湖泊水生植物监测结果图。本发明能够快速、实时、大范围和高精度的识别和监测不同种类水生植物的空间分布。

Description

一种基于HJ-CCD影像的浅水湖泊水生植物快速监测方法

技术领域

本发明属于遥感应用领域，涉及一种浅水湖泊水生植物的遥感快速监测方法。

背景技术

水生植物是湖泊生态系统重要的组成部分，也是湖泊生态系统的重要调节者和指示者。水生植物的实时、快速监测对湖泊生态系统生态功能的评估具有指导意义，同时，对浊水态湖泊生态系统的修复和水生植物的打捞等具有重要的现实意义。

目前水生植物的监测方法仍然沿用人工样方、样线调查法，其费时费力、监测范围和样本量有限，监测难度大（很多水生植物生长的地方人很难到达），无法快速、实时、大面积的获取水生植物的类型和分布状况；因此，提出一种基于环境资源小卫星（HJ-CCD）影像的浅水湖泊水生植物快速监测方法，，开展浅水湖泊水生植物类型和分布的快速、实时、大范围和高精度的识别和监测势在必行。

发明内容

为了解决传统水生植物监测方法代表性、时效性差；费时费力、监测范围和样本量有限，监测难度大等问题和不足，本发明的提供了一种基于HJ-CCD影像的浅水湖泊水生植物快速监测方法，所述水生植物包括挺水植物、浮叶植物和沉水植物，其包括以下步骤：

步骤一、获取研究区HJ-CCD影像，并进行预处理；

步骤二、根据所述研究区HJ-CCD影像，计算挺水植物敏感特征指数EVSI；

步骤三、分别对所述研究区HJ-CCD影像进行主成分变换，获取浮叶植物敏感特征指数FVSI；

步骤四、对所述研究区HJ-CCD影像进行缨帽变换，获取沉水植物敏感特征指数SVSI；

步骤五、根据所述研究区HJ-CCD影像同步的实测数据，确定EVSI、FVSI和SVSI的阈值，并得到判别条件；

步骤六、根据判别条件建立决策树分类模型，生成浅水湖泊水生植物监测结果，并利用实测数据进行验证。

其中，所述步骤一中的影像预处理包括辐射校正、几何校正和研究区裁剪，其中所述辐射校正指辐射定标和大气校正。

其中，所述步骤二中的EVSI计算公式如下：

$\mathrm{EVSI}=\frac{R_{\mathrm{nir}}-(R_{\mathrm{blue}}+R_{\mathrm{green}}+R_{\mathrm{red}})}{R_{\mathrm{nir}}+(R_{\mathrm{blue}}+R_{\mathrm{green}}+R_{\mathrm{red}})}$

其中,R_blue为蓝波段，其波长范围为0.45-0.52微米；R_green为绿波段，其波长范围为0.52-0.60微米；R_red为红波段，其波长范围为0.63-0.69微米；R_nir为近红外波段，其波长范围为0.76-0.90微米。

其中，所述步骤三中的FVSI计算公式如下：

FVSI=PC₂，

其中，PC₂为研究区HJ-CCD影像主成分变换后的第二主成分。

其中，所述的PC₂获取步骤如下：利用获取的研究区HJ-CCD影像，进行主成分变换，所述主成分变换的步骤如下：

步骤（a）、以矩阵的形式来表示多光谱影像的原始数据，则设获取的研究区HJ-CCD影像数据矩阵为X，如下：

$X=\left[\begin{array}{cccc}{x}_{11}& x_{12}& \·\·\·& x_{1n}\\ x_{21}& x_{22}& \·\·\·& x_{2n}\\ \·& \·& & \·\\ \·& \·& & \·\\ \·& \·& & \·\\ x_{p1}& x_{p2}& \·\·\·& x_{\mathrm{pn}}\end{array}\right]={\left[x_{\mathrm{ik}}\right]}_{p\×n}$

其中，n为HJ-CCD影像的波段数，n=1,2,3,4；p HJ-CCD影像每个波段影像中的像元数，矩阵中每一行矢量表示一个波段的图像；

步骤（b）、基于矩阵X，通过Y=TX对图像进行变换得到主成分变换后的图像数据矩阵Y，T为变换矩阵，求解变换矩阵T；

所述的变换矩阵T为X空间协方差∑x_i的特征向量矩阵的转置矩阵，其求解步骤如下：

1）根据原始图像数据矩阵X，求出X的协方差矩阵S为：

$S=\left[\begin{array}{cccc}{s}_{11}& s_{12}& \·\·\·& s_{1p}\\ s_{21}& s_{22}& \·\·\·& x_{2p}\\ \·& \·& & \·\\ \·& \·& & \·\\ \·& \·& & \·\\ s_{p1}& s_{p2}& \·\·\·& s_{\mathrm{pp}}\end{array}\right]={\left[s_{\mathrm{ij}}\right]}_{p\×p}$

所述的S_ij表达式为：

$S_{\mathrm{ij}}=\frac{1}{n}\underset{k=1}{\overset{n}{\mathrm{\Σ}}}(x_{\mathrm{ik}}-\stackrel{\‾}{x_i})(x_{\mathrm{ij}}-\stackrel{\‾}{x_j})\stackrel{\‾}{x_i}=\frac{1}{n}\underset{k=1}{\overset{n}{\mathrm{\Σ}}}{x}_{\mathrm{ik}},$ 其中，k=1,2,3,….n

所述的为图像各波段的均值，表达式如下：

其中，i=1,2，…，n；j=1,2，…，p

2)协方差S的特征值λ_i和特征向量U_j，并组成变换矩阵T，，求解特征方程(λ_i-S)U=0；然后将特征值λ_i由小到大排列，求出对应特征值的单位特征向量U_j：

U_j=[u_1j u_2j ... u_pj]^T，

以各特征向量为列构成矩阵，即：

$U=\left[\begin{array}{cccc}{U}_1& U_2& ...& U_p\end{array}\right]={\left[u_{\mathrm{ij}}\right]}_{p\×p};$

3）U矩阵的转置矩阵，即U^T为所求的变换矩阵T；

步骤（c）、将变换矩阵T带入Y=TX，求解主成分变换后的矩阵Y：

$Y=\left[\begin{array}{cccc}{u}_{11}& u_{21}& \·\·\·& u_{p1}\\ u_{12}& u_{22}& \·\·\·& u_{p2}\\ \·& \·& & \·\\ \·& \·& & \·\\ \·& \·& & \·\\ u_{1p}& u_{2p}& \·\·\·& u_{\mathrm{pp}}\end{array}\right]X=U^{T}X,$

步骤（d）、求解PC₂，

令Y₂=PC₂，Y矩阵的行向量Y₂=[y_j1 y_j2 ... y_jp]为第j主成分；Y₁，Y₂….Y_p即为研究区HJ-CCD影像经过主成分变换后的新变量，以此被称为第一主成分，第二主成分第p主成分，将新变量恢复为二维图像，即得到p个主成分图像，因此，

Y₂=[y₂₁ y₂₂ ... y_2p]

将Y₂恢复为二维图像，即得到PC₂图像。

其中，所述的步骤四的SVSI计算公式如下：

SVSI=BI-GVI，

其中，BI为亮度指数，其为缨帽变换后的第一分量；GVI为绿度指数，其为缨帽变换后的第二分量。

其中，亮度指数BI和绿度指数GVI求算步骤如下：

通过以下公式进行缨帽变换：

y=Cx+α，

α为表示常数偏移量，是为避免在变换过程中出现的负值；x为变换前的HJ-CCD多光谱影像的像元矢量；y为HJ-CCD变换后的多光谱影像的像元矢量；C为变换矩阵，如下：

$C=\left[\begin{array}{cccc}0.433& 0.632& 0.586& 0.264\\ -0.290& -0.562& -0.600& 0.491\\ -0.829& 0.522& -0.039& 0.194\\ 0.223& 0.012& 0.543& 0.810\end{array}\right]$

因此，可得：

BI=0.433R_b+0.632R_g+0.586R_r+0.264R_Nir；

GVI=-0.290R_b+(-0.562R_g)+(-0.600R_r)+0.491R_Nir。

其中，所述步骤六包括以下步骤：

步骤6.1、当研究区像元满足EVSI>a and D<500m时，则被判别为挺水植物；

步骤6.2、在不满足步骤6.1的像元中进行判别，当FVSI<b时，则被判别为浮叶植物；

步骤6.3、在不满足步骤6.1和6.2中条件的像元中进行判别，当SVSI<c时，则被判别为沉水植物；

步骤6.4、不满足以上三个步骤的条件的像元被判别为水体；

其中，D为离岸的距离；a为EVSI识别挺水植物的阈值，b为FVSI识别浮叶植物的阈值，c为SVSI识别沉水植物的阈值；a，b和c通过与影像同步的实测数据训练获取。

其中，水生植物监测精度利用总体精度和Kappa系数来评价，其计算公式如下：

总体分类精度 $\left(\mathrm{OA}\right)=\underset{k=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&beta;}}_{\mathrm{kk}}/\mathrm{\&beta;},$

其中，为与实际样本在图像上被正确分类的样本总数，β为样本总数；

$\mathrm{kappa}=\frac{N\underset{i}{\overset{r}{\mathrm{\&Sigma;}}}{A}_{\mathrm{ii}}-\mathrm{\&Sigma;}(A_{i+}+...A_{+i})}{N^2-\mathrm{\&Sigma;}(A_{i+}+...A_{+i})}$

其中，r是错误矩阵中总列数，即总的类别数；A_ii是错误矩阵中的第i行、第j列上像元数量，即正确分类的数目；A_i+和A_+i分别是第i和第i列的总像元数量；N是总的用于精度评估的像元数量。

附图说明

图1为本发明技术方案的流程图；

图2为本发明所述的方法实施例中2013年5月13日的洪泽湖水生植物监测结果；

图3为本发明所述的方法实施例中2013年7月11日的太湖水生植物监测结果；

图4为本发明所述的方法实施例中2013年9月26日的太湖水生植物监测结果。

具体的实施方式

下面结合附图和实施例对本发明进行进一步说明。一种基于HJ-CCD影像的浅水湖泊水生植物快速监测方法，以中国典型的大型浅水湖泊（平均深度为1.9米）——太湖和洪泽湖为例进行实例分析。其中，太湖经纬度范围为119.55-120.34E，30.55-31.32N;；洪泽湖经纬度范围为118.18-118.52E，北纬33.22～33.37N。经调研和资料显示，太湖东部和洪泽湖分布有大面积的水生植物，主要为以下三大类：挺水植物、浮叶植物和沉水植物。目前，由于近年来人类活动剧烈，湖泊水生植物和湿地环境造成严重的干扰和破坏，利用遥感技术对太湖和洪泽湖的水生植物分布区域及面积进行实时快速的监测，对湖泊生态修复和水生植物打捞具有重要的指导意义。下面以太湖和洪泽湖为研究区，进行实例分析：

步骤一、根据实测样点获取时间，从中国资源卫星应用中心的数据产品查询订购网站(http://218.247.138.121/DSSPlatform/index.html)订购下载与实测样点同步的覆盖太湖的2013年7月11日和覆盖洪泽湖的2013年5月13日的环境资源小卫星影像数据(HJ-CCD)。

步骤二、利用ENVI软件分别对两景影像进行大气校正（辐射定标和大气校正）和几何校正，然后利用太湖和洪泽湖矢量边界从影像上裁剪出相应的研究区影像；

步骤三、根据EVSI的计算公式获取两景影像的EVSI图像；

步骤四、对2013年5月13日、2013年7月11日和2013年9月26日的研究区影像分别进行主成分变换，并根据公式FVSI=PC₂，分别获取其FVSI图像；

步骤五、对2013年5月13日、2013年7月11日和2013年9月26日研究区影像分别进行缨帽变换，得到研究区的亮度指数BI和绿度指数GVI，并根据公式SVSI=BI-GVI，分别获取其SVSI图像；

步骤六、利用EVSI、FVSI和SVSI，构建识别不同水生植物的判别式，如下：

条件1、当研究区像元满足EVSI>a and D<500m时，则被判别为挺水植物；

条件2、在不满足条件1的像元中，当FVSI<b时，则被判别为浮叶植物；

条件3、在不满足1和2的像元中，当SVSI<c时，则被判别为沉水植物；

条件4、不满足以上步骤中个条件的像元被判别为水体。

其中，D为离岸的距离；a为EVSI识别挺水植物的阈值，b为FVSI识别浮叶植物的阈值，c为SVSI识别浮叶植物的阈值；a，b和c通过与影像同步的实测数据训练获取；表1为2013年5月13日洪泽湖和2013年7月11日与2013年9月26日太湖水生植物监测的各特征阈值。

表1洪泽湖和太湖各特征阈值

步骤七、在ArcGIS里将分类结果进行出图，如图3，并根据7月11日和8月16日的实测样点和分类结果，计算分类总体精度和Kappa系数，评价分类结果，如表2。

表2是本发明实施例所述的方法获取的太湖水生植物监测结果的精度评价表

表2是本发明实施例所述的方法获取的太湖水生植物监测结果的精度评价表，通过监测结果与实测结果的对比可知，本发明方法对浅水湖泊（太湖和洪泽湖）水生植物的监测精度均高于80%，Kappa系数大于0.8，符合大部分应用需求。此外，图2是本发明实施例所述的方法获取的洪泽湖水生植物监测结果示意图；图3和图4是本发明实施例所述的方法获取的太湖水生植物监测结果示意图。可以看出，本发明相比现有的传统方法的优势为：利用有限的实测样点和遥感影像，能够快速、准确的获取研究区水生植物的分布类型和分布区域。

尽管已经示出和描述了本发明的实施例，本领域的普通技术人员可以理解：在不脱离本发明的原理和宗旨的情况下可以对这些实施例进行多种变化、修改、替换和变型，本发明的范围由权利要求及其等同物限定。

Claims (9)

1.一种基于HJ-CCD影像的浅水湖泊水生植物快速监测方法，所述水生植物包括挺水植物、浮叶植物和沉水植物，其包括以下步骤：

步骤一、获取研究区HJ-CCD影像，并进行预处理；

步骤二、根据所述研究区HJ-CCD影像，计算挺水植物敏感特征指数EVSI；

步骤三、分别对所述研究区HJ-CCD影像进行主成分变换，获取浮叶植物敏感特征指数FVSI；

步骤四、对所述研究区HJ-CCD影像进行缨帽变换，获取沉水植物敏感特征指数SVSI；

步骤五、根据所述研究区HJ-CCD影像同步的实测数据，确定EVSI、FVSI和SVSI的阈值，并得到判别条件；

步骤六、根据判别条件建立决策树分类模型，生成浅水湖泊水生植物监测结果，并利用实测数据进行验证。

2.根据权利要求1所述的方法，所述步骤一中的影像预处理包括辐射校正、几何校正和研究区裁剪，其中所述辐射校正指辐射定标和大气校正。

3.根据权利要求1所述的方法，所述步骤二中的EVSI计算公式如下：

$\mathrm{EVSI}=\frac{R_{\mathrm{nir}}-(R_{\mathrm{blue}}+R_{\mathrm{green}}+R_{\mathrm{red}})}{R_{\mathrm{nir}}+(R_{\mathrm{blue}}+R_{\mathrm{green}}+R_{\mathrm{red}})}$

其中,R_blue为蓝波段，其波长范围为0.45-0.52微米；R_green为绿波段，其波长范围为0.52-0.60微米；R_red为红波段，其波长范围为0.63-0.69微米；R_nir为近红外波段，其波长范围为0.76-0.90微米。

4.根据权利要求1所述方法，所述步骤三中的FVSI计算公式如下：

FVSI=PC₂，

其中，PC₂为研究区HJ-CCD影像主成分变换后的第二主成分。

5.根据权利要求4所述的方法，所述的PC₂获取步骤如下：利用获取的研究区HJ-CCD影像，进行主成分变换，所述主成分变换的步骤如下：

步骤（a）、以矩阵的形式来表示多光谱影像的原始数据，则设获取的研究区HJ-CCD影像数据矩阵为X，如下：

$X=\left[\begin{array}{cccc}{x}_{11}& x_{12}& \&CenterDot;\&CenterDot;\&CenterDot;& x_{1n}\\ x_{21}& x_{22}& \&CenterDot;\&CenterDot;\&CenterDot;& x_{2n}\\ \&CenterDot;& \&CenterDot;& & \&CenterDot;\\ \&CenterDot;& \&CenterDot;& & \&CenterDot;\\ \&CenterDot;& \&CenterDot;& & \&CenterDot;\\ x_{p1}& x_{p2}& \&CenterDot;\&CenterDot;\&CenterDot;& x_{\mathrm{pn}}\end{array}\right]={\left[x_{\mathrm{ik}}\right]}_{p\&times;n}$

其中，n为HJ-CCD影像的波段数，n=1,2,3,4；p HJ-CCD影像每个波段影像中的像元数，矩阵中每一行矢量表示一个波段的图像；

步骤（b）、基于矩阵X，通过Y=TX对图像进行变换得到主成分变换后的图像数据矩阵Y，T为变换矩阵，求解变换矩阵T；

所述的变换矩阵T为X空间协方差∑x_i的特征向量矩阵的转置矩阵，其求解步骤如下：

1）根据原始图像数据矩阵X，求出X的协方差矩阵S为：

$S=\left[\begin{array}{cccc}{s}_{11}& s_{12}& \&CenterDot;\&CenterDot;\&CenterDot;& s_{1p}\\ s_{21}& s_{22}& \&CenterDot;\&CenterDot;\&CenterDot;& x_{2p}\\ \&CenterDot;& \&CenterDot;& & \&CenterDot;\\ \&CenterDot;& \&CenterDot;& & \&CenterDot;\\ \&CenterDot;& \&CenterDot;& & \&CenterDot;\\ s_{p1}& s_{p2}& \&CenterDot;\&CenterDot;\&CenterDot;& s_{\mathrm{pp}}\end{array}\right]={\left[s_{\mathrm{ij}}\right]}_{p\&times;p}$

所述的S_ij表达式为：

$S_{\mathrm{ij}}=\frac{1}{n}\underset{k=1}{\overset{n}{\mathrm{\&Sigma;}}}(x_{\mathrm{ik}}-\stackrel{\&OverBar;}{x_i})(x_{\mathrm{ij}}-\stackrel{\&OverBar;}{x_j})\stackrel{\&OverBar;}{x_i}=\frac{1}{n}\underset{k=1}{\overset{n}{\mathrm{\&Sigma;}}}{x}_{\mathrm{ik}},$ 其中，k=1,2,3,….n

所述的为图像各波段的均值，表达式如下：

其中，i=1,2，…，n；j=1,2，…，p

2)协方差S的特征值λ_i和特征向量U_j，并组成变换矩阵T，，求解特征方程(λ_i-S)U=0；然后将特征值λ_i由小到大排列，求出对应特征值的单位特征向量U_j：

U_j=[u_1j u_2j ... u_pj]^T，

以各特征向量为列构成矩阵，即：

$U=\left[\begin{array}{cccc}{U}_1& U_2& ...& U_p\end{array}\right]={\left[u_{\mathrm{ij}}\right]}_{p\&times;p};$

3）U矩阵的转置矩阵，即U^T为所求的变换矩阵T；

步骤（c）、将变换矩阵T带入Y=TX，求解主成分变换后的矩阵Y：

$Y=\left[\begin{array}{cccc}{u}_{11}& u_{21}& \&CenterDot;\&CenterDot;\&CenterDot;& u_{p1}\\ u_{12}& u_{22}& \&CenterDot;\&CenterDot;\&CenterDot;& u_{p2}\\ \&CenterDot;& \&CenterDot;& & \&CenterDot;\\ \&CenterDot;& \&CenterDot;& & \&CenterDot;\\ \&CenterDot;& \&CenterDot;& & \&CenterDot;\\ u_{1p}& u_{2p}& \&CenterDot;\&CenterDot;\&CenterDot;& u_{\mathrm{pp}}\end{array}\right]X=U^{T}X,$

步骤（d）、求解PC₂

令Y₂=PC₂，Y矩阵的行向量Y₂=[y_j1 y_j2 ... y_jp]为第j主成分；Y₁，Y₂….Y_p即为研究区HJ-CCD影像经过主成分变换后的新变量，以此被称为第一主成分，第二主成分第p主成分，将新变量恢复为二维图像，即得到p个主成分图像，因此，

Y₂=[y₂₁ y₂₂ ... y_2p]

将Y₂恢复为二维图像，即得到PC₂图像。

6.根据权利要求1所述的方法，其特征在于，所述的步骤四的SVSI计算公式如下：

SVSI=BI-GVI，

其中，BI为亮度指数，其为缨帽变换后的第一分量；GVI为绿度指数，其为缨帽变换后的第二分量。

7.根据权利要求6所述的方法，其特征在于，亮度指数BI和绿度指数GVI求算步骤如下：

通过以下公式进行缨帽变换：

y=Cx+α，

α为表示常数偏移量，是为避免在变换过程中出现的负值；x为变换前的HJ-CCD多光谱影像的像元矢量；y为HJ-CCD变换后的多光谱影像的像元矢量；C为变换矩阵，如下：

$C=\left[\begin{array}{cccc}0.433& 0.632& 0.586& 0.264\\ -0.290& -0.562& -0.600& 0.491\\ -0.829& 0.522& -0.039& 0.194\\ 0.223& 0.012& 0.543& 0.810\end{array}\right]$

因此，可得：

BI=0.433R_b+0.632R_g+0.586R_r+0.264R_Nir；

GVI=-0.290R_b+(-0.562R_g)+(-0.600R_r)+0.491R_Nir。

8.根据权利要求1所述的方法，所述步骤六包括以下步骤：

步骤6.1、当研究区像元满足EVSI>a and D<500m时，则被判别为挺水植物；

步骤6.2、在不满足步骤6.1的像元中进行判别，当FVSI<b时，则被判别为浮叶植物；

步骤6.3、在不满足步骤6.1和6.2中条件的像元中进行判别，当SVSI<c时，则被判别为沉水植物；

步骤6.4、不满足以上三个步骤的条件的像元被判别为水体；

其中，D为离岸的距离；a为EVSI识别挺水植物的阈值，b为FVSI识别浮叶植物的阈值，c为SVSI识别沉水植物的阈值；a，b和c通过与影像同步的实测数据训练获取。

9.根据权利要求1所述的方法，水生植物监测精度利用总体精度和Kappa系数来评价，其计算公式如下：

总体分类精度 $\left(\mathrm{OA}\right)=\underset{k=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&beta;}}_{\mathrm{kk}}/\mathrm{\&beta;},$

其中，为与实际样本在图像上被正确分类的样本总数，β为样本总数；

$\mathrm{kappa}=\frac{N\underset{i}{\overset{r}{\mathrm{\&Sigma;}}}{A}_{\mathrm{ii}}-\mathrm{\&Sigma;}(A_{i+}+...A_{+i})}{N^2-\mathrm{\&Sigma;}(A_{i+}+...A_{+i})}$

其中，r是错误矩阵中总列数，即总的类别数；A_ii是错误矩阵中的第i行、第j列上像元数量，即正确分类的数目；A_i+和A_+i分别是第i和第i列的总像元数量；N是总的用于精度评估的像元数量。