CN102829739B - 一种面向对象的作物叶面积指数遥感反演方法 - Google Patents

Abstract

本发明公开了一种面向对象的作物叶面积指数遥感反演方法，包括：获取多光谱遥感数据；利用获取的多光谱遥感数据计算作物群体生物量光谱指数NDVI、作物养分光谱指数BRI和水分敏感光谱指数NDWI；依据作物群体生物量光谱指数NDVI、作物养分光谱指数BRI和水分敏感光谱指数NDWI，利用均值漂移算法进行面向对象分割及编码；按照编码顺序依次对各对象进行像元原始光谱均值并得出与LAI敏感的光谱指数SAVI，以及纹理结构的计算；建立地面LAI观测数据、与LAI敏感的光谱指数SAVI和纹理结构计算的回归模型；利用回归模型对没有地面LAI观测数据的对象进行反演计算，获得没有地面LAI观测数据对象的LAI。

Description

一种面向对象的作物叶面积指数遥感反演方法

技术领域

本发明涉及农业信息技术领域，特别涉及一种面向对象的作物叶面积指数遥感反演方法。

背景技术

利用遥感技术进行大范围作物叶面积指数（LAI,Leaf areaindex）反演已经成为较为成熟的技术手段和方法,主要利用多光谱遥感筛选出与LAI敏感波段，进而利用敏感波段计算出光谱指数，再利用实地测量的LAI数据和光谱指数建立模型，然后将模型应用到其他像元对应的光谱指数上，就可以得到完整影像对应区域的所有像元LAI。可以看出，LAI的反演计算是逐个像元进行的，在反演大范围LAI或多个时相LAI时，计算效率比较低。对于作物来说，由于在局地相同的地理、气象条件下，相邻像元间作物长势的差异并不是十分明显，而是以斑块状同质像元呈现，斑块内部作物长势总体上是一致的，所以再按照逐个像元反演LAI显然是没有必要的。如果将作物LAI反演由逐像元改为逐个对象计算，则计算效率将会极大提高。这里所说的对象既包括了斑块的空间范围，也含有斑块区域内像元的空间纹理结构特征和光谱特征。关于面向对象的图像分割技术已经较为成熟，并且已成功应用到遥感领域，在考虑像元光谱特征基础上，还要考虑图像纹理结构、上下文关系等，实现多尺度遥感图像分割。但目前关于面向对象的遥感图像分割侧重于地表分类或者地物识别精度的提高，针对作物长势状况的面向对象分割还欠缺。所以，把面向对象的遥感图像分割技术引入到LAI反演中来更有实际应用价值。

对于背景技术方案，存在的主要缺陷为：

（1）目前，所有利用遥感数据反演作物LAI的算法都是基于逐像元计算的，由于遥感器观测存在一定噪声，且遥感数据处理过程中也存在一定误差，不可避免的存在部分像元反演出的LAI是错误的，这给LAI的后续分析和应用带来较大问题；

（2）在进行作物长势诊断及变量肥水实施过程中，需要对以‘像元’为基础进行LAI空间插值处理，形成作物长势空间分布的处方图。因此，插值处理会带来一定的误差；而且这种单独依赖LAI的长势诊断及处方图可靠性不高，主要是由于LAI只是体现了作物群体生物量的大小，必须还要结合表征作物养分、水分的相关信息，才能全面的表达作物长势情况。

发明内容

（一）要解决的技术问题

本发明要解决的技术问题是利用遥感数据事先得到关于作物生物量、养分及水分的空间分布特征，并对其进行面向对象的分割和提取，得到具有不同生物量、养分及水分组合特征的空间单元，这种单元虽然大小各异，但相同类别单元对应的作物具有相同的生物量、养分及水分信息，在此基础上针对各对象再进行LAI反演。

（二）技术方案

一种面向对象的作物叶面积指数遥感反演方法，包括：

获取多光谱遥感数据；

利用获取的多光谱遥感数据计算作物群体生物量光谱指数NDVI、作物养分光谱指数BRI和水分敏感光谱指数NDWI；

依据作物群体生物量光谱指数NDVI、作物养分光谱指数BRI和水分敏感光谱指数NDWI，利用均值漂移算法进行面向对象分割及编码；

按照编码顺序依次对各对象进行像元原始光谱均值并得出与LAI敏感的光谱指数SAVI，以及纹理结构的计算；

建立地面LAI观测数据、与LAI敏感的光谱指数SAVI和纹理结构计算的回归模型；

利用回归模型对没有地面LAI观测数据的对象进行反演计算，获得没有地面LAI观测数据对象的LAI。

其中，

所述作物群体生物量光谱指数NDVI的计算公式为 $\mathrm{NDVI}=\frac{R_{\mathrm{nir}}-R_{\mathrm{red}}}{R_{\mathrm{nir}}+R_{\mathrm{red}}},$

所述作物养分光谱指数BRI的计算公式为

所述水分敏感光谱指数NDWI的计算公式为

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

其中，所述利用均值漂移算法进行面相对象分割的过程为：

依据均值漂移分割算法核函数公式：

$K_{h_s,h_r}\left(x\right)=\frac{C}{h_s^2{h}_r^p}k\left({\left|\right|\frac{x^s}{h_s}\left|\right|}^2\right)k\left({\left|\right|\frac{x^r}{h_r}\left|\right|}^2\right),$

其中，C—归一化常数；k—核函数；h_s—MS算法的空间带宽系数；h_r—MS算法的光谱指数带宽系数；通过控制核带宽参数h＝(h_s,h_r)来决定分割精度，通常设置h_s＝7，h_r＝6.5；

x^s表示空间坐标，x^r表示r维特征向量；分量x^s∈{0＜i≤w,0＜j≤h}，h、w分别是图像行高与列宽；x^r∈{L,u,v}，L、u、v分别对应NDVI、BRI及NDWI经过色彩空间LUV变换后的三个色彩分量。

其中，所述光谱敏感指数SAVI的计算公式为，

$\mathrm{SAVI}=\frac{(R_{\mathrm{nir}}-R_{\mathrm{red}})}{(R_{\mathrm{nir}}+R_{\mathrm{red}}+L)}(1+L);$

其中，R_red为红波段，其波长范围为0.63-0.69微米；R_nir为近红外波段，其波长范围为0.76-0.90微米；L为对应NDVI经过色彩空间LUV变换后的色彩分量。

其中，所述纹理结构计算包括熵值计算和同质度计算；

所述熵值计算的公式为：熵值 $\mathrm{Ent}=-\underset{i=0}{\overset{N}{\mathrm{\Σ}}}\underset{j=0}{\overset{N}{\mathrm{\Σ}}}{C}_{\mathrm{ij}}\log \frac{P_{\mathrm{ij}}(d,\mathrm{\θ})}{\underset{i=0}{\overset{N}{\mathrm{\Σ}}}\underset{j=0}{\overset{N}{\mathrm{\Σ}}}{P}_{\mathrm{ij}}(d,\mathrm{\θ})};$

所述同质度的计算公式为：

同质度 $\mathrm{Homo}=\underset{i=0}{\overset{N}{\mathrm{\Σ}}}\underset{j=0}{\overset{N}{\mathrm{\Σ}}}\frac{1}{1+{(i-j)}^2}\frac{P_{\mathrm{ij}}(d,\mathrm{\θ})}{\underset{i=0}{\overset{N}{\mathrm{\Σ}}}\underset{j=0}{\overset{N}{\mathrm{\Σ}}}{P}_{\mathrm{ij}}(d,\mathrm{\θ})};$

其中，P_i,j(d,θ)为利用对象内像元NDVI、BRI及NDWI计算的灰度共生矩阵，N为对象内部像元数目；i、j为自然数。

其中，所述回归模型为

LAI＝a_m·[SAVI Ent_NDVI Ent_BRI Ent_NDWI Homo_NDVI Homo_BRI Homo_NDWI]^T其中a_m＝[a₀ a₁ a₂ a₃ a₄ a₅ a₆]为回归系数，SAVI为光谱敏感指数，R_blue为蓝波段,R_green为绿波段,R_red为红波，；R_nir为近红外波段，Ent为熵值，Homo为同质度。

（三）有益效果

本发明提出的技术方案直接面向农田作业变量处方图生成目标，利用多光谱遥感数据提取与作物生物量、养分及水分相关的光谱指数，并基于面向对象的图像分割方法实现三种组合光谱指数的分割，通过对分割对象进行统一编码即可获得作物长势状况空间分布图。在此基础上充分利用遥感数据的光谱信息和空间纹理信息，并建立作物LAI与光谱指数(SAVI)、纹理参数（熵、同质度）间的回归模型，将模型应用于所有分割对象实现面向对象的LAI反演。此方法适合大多数具有相同光谱波段的卫星、航空及地面遥感数据，由于同时考虑了光谱信息和纹理信息，避免了逐像元反演中由于像元本身误差造成的精度较低等问题；另外，作物LAI反演是基于对象进行，大幅提高了LAI反演的计算效率，而且反演结果可直接作物农田管理的作业处方图，避免了逐像元反演再空间插值带来的误差。

附图说明

图1是本发明技术方案的流程图。

具体实施方式

下面结合附图和实施例，对本发明的具体实施方式作进一步详细描述。以下实施例用于说明本发明，但不用来限制本发明的范围。

图1是本发明技术方案的主流程图。如图1所示，本发明提出面向对象的作物LAI反演方法包括下列步骤：

（1）获取多光谱遥感数据，数据中包含以下四个波段反射率数据：蓝波段（Rblue）,波长范围为0.45-0.52微米；绿波段(Rgreen),波长范围为0.52-0.60微米；红波段（Rred），波长范围为0.63-0.69微米；近红外波段（Rnir），波长范围为0.76-0.90微米。同时获取地面观测LAI数据，并利用差分GPS进行观测点精准定位，保证定位误差小于5米。

（2）利用遥感数据四个波段反射率光谱计算与作物群体生物量、作物养分及水分敏感的光谱指数NDVI、BRI及NDWI，具体公式为：

$\mathrm{NDVI}=\frac{R_{\mathrm{nir}}-R_{\mathrm{red}}}{R_{\mathrm{nir}}+R_{\mathrm{red}}}---\left(1\right)$

$\mathrm{BRI}=\frac{R_{\mathrm{blue}}}{R_{\mathrm{red}}}---\left(2\right)$

$\mathrm{NDWI}=\frac{R_{\mathrm{nir}}-R_{\mathrm{green}}}{R_{\mathrm{nir}}+R_{\mathrm{green}}}---\left(3\right)$

（3）利用计算的NDVI、BRI及NDWI,选择均值漂移（Mean-Shift，MS）算法进行面向对象分割。由于MS算法完全依靠特征空间中的样本点进行分析，不需要任何先验知识，且收敛速度快，近年来被广泛应用于图像分割。

考虑NDVI、BRI及NDWI三个光谱指数（r=3），同时考虑三者的空间信息（s=2），就组成一个r+s=5维的向量x＝(x^s,x^r)。其中x^s表示空间坐标，x^r表示r维特征向量，则MS分割算法核函数公式：

$K_{h_s,h_r}\left(x\right)=\frac{C}{h_s^2{h}_r^p}k\left({\left|\right|\frac{x^s}{h_s}\left|\right|}^2\right)k\left({\left|\right|\frac{x^r}{h_r}\left|\right|}^2\right)---\left(4\right)$

式中：C—归一化常数；k—核函数；h_s—MS算法的空间带宽系数；h_r—MS算法的光谱指数带宽系数；通过控制核带宽参数h＝(h_s,h_r)来决定分割精度，通常设置h_s＝7，h_r＝6.5。分量x^s∈{0＜i≤w,0＜j≤h}，h、w分别是图像行高与列宽；x^r∈{L,u,v}，L、u、v分别对应NDVI、BRI及NDWI经过色彩空间LUV变换后的三个色彩分量。

利用公式（4）经过逐像元迭代计算，即可得到分割后的对象。在此基础上，对不同类别对象进行编码，并保证相同类别的对象具有相同编码，这就意味着具有相同编码的对象所对应的NDVI、BRI及NDWI相同，对应的作物生物量、养分及水分状况相同。

（4）经过对分割对象及编码后，可以按照编码顺序依次对各对象进行光谱及空间特征计算。首先根据各对象所包含像元对应的原始多光谱反射率进行均值计算，得到代表对象的平均反射光谱，进而计算与LAI敏感的光谱指数SAVI:

$\mathrm{SAVI}=\frac{(R_{\mathrm{nir}}-R_{\mathrm{red}})}{(R_{\mathrm{nir}}+R_{\mathrm{red}}+L)}(1+L)---\left(5\right)$

另外，对分割对象内部各像元NDVI、BRI及NDWI值进行纹理结构计算，包括熵和同质度。

熵（Entropy，Ent）是反映图像中纹理的复杂程度或非均匀度。熵大时，纹理复杂；反之，纹理简单。公式为：

$\mathrm{Ent}=-\underset{i=0}{\overset{N}{\mathrm{\Σ}}}\underset{j=0}{\overset{N}{\mathrm{\Σ}}}{C}_{\mathrm{ij}}\log \frac{P_{\mathrm{ij}}(d,\mathrm{\θ})}{\underset{i=0}{\overset{N}{\mathrm{\Σ}}}\underset{j=0}{\overset{N}{\mathrm{\Σ}}}{P}_{\mathrm{ij}}(d,\mathrm{\θ})}---\left(6\right)$

同质度（Homogeneity，Homo），又称逆差矩，用来度量图像纹理局部变化的多少，反映图像的均质性。Homo大时，纹理不均匀；反之，Homo小时，纹理均匀。公式为：

$\mathrm{Homo}=\underset{i=0}{\overset{N}{\mathrm{\Σ}}}\underset{j=0}{\overset{N}{\mathrm{\Σ}}}\frac{1}{1+{(i-j)}^2}\frac{P_{\mathrm{ij}}(d,\mathrm{\θ})}{\underset{i=0}{\overset{N}{\mathrm{\Σ}}}\underset{j=0}{\overset{N}{\mathrm{\Σ}}}{P}_{\mathrm{ij}}(d,\mathrm{\θ})}---\left(7\right)$

其中P_i,j(d,θ)为利用对象内像元NDVI、BRI及NDWI计算的灰度共生矩阵，N为对象内部像元数目。

（5）利用地面观测LAI数据，并根据观测点位置查找到对应的分割对象，提取对象对应的SAVI及Ent_NDVI、Ent_BRI、Ent_NDWI、Homo_NDVI、Homo_BRI、Homo_NDWI值，并建立相应的回归模型：

LAI＝a_m·[SAVI Ent_NDVI Ent_BRI Ent_NDWI Homo_NDVI Homo_BRI Homo_NDWI]^T（8）其中a_m＝[a₀ a₁ a₂ a₃ a₄ a₅ a₆]为回归系数。

（6）利用公式（8）对其他没有地面观测的对象进行LAI计算，即可得到所有对象的LAI值，完成面向对象的LAI反演。

以上实施方式仅用于说明本发明，而并非对本发明的限制，有关技术领域的普通技术人员，在不脱离本发明的精神和范围的情况下，还可以做出各种变化和变型，因此所有等同的技术方案也属于本发明的范畴，本发明的专利保护范围应由权利要求限定。

Claims (5)

1.一种面向对象的作物叶面积指数遥感反演方法，其特征在于，包括：

获取多光谱遥感数据；

利用获取的多光谱遥感数据计算作物群体生物量光谱指数NDVI、作物养分光谱指数BRI和水分敏感光谱指数NDWI；

依据作物群体生物量光谱指数NDVI、作物养分光谱指数BRI和水分敏感光谱指数NDWI，利用均值漂移算法进行面向对象分割及编码；

按照编码顺序依次对各对象进行像元原始光谱均值并得出与叶面积指数LAI敏感的光谱指数SAVI，以及纹理结构的计算；

建立地面LAI观测数据、与LAI敏感的光谱指数SAVI和纹理结构计算的回归模型；

利用回归模型对没有地面LAI观测数据的对象进行反演计算，获得没有地面LAI观测数据对象的LAI；

其中，所述利用均值漂移算法进行面向对象分割的过程为：

依据均值漂移分割算法核函数公式：

$K_{h_s,h_r}\left(x\right)=\frac{C}{h_s^2{h}_r^p}k\left({\left|\right|\frac{x^s}{h_s}\left|\right|}^2\right)k\left({\left|\right|\frac{x^r}{h_r}\left|\right|}^2\right),$

其中，C—归一化常数；k—核函数；h_s—MS算法的空间带宽系数；h_r—MS算法的光谱指数带宽系数；通过控制核带宽参数h＝(h_s,h_r)来决定分割精度，设置h_s＝7，h_r＝6.5；

x^s表示空间坐标，x^r表示r维特征向量；分量x^s∈{0<i≤w,0<j≤h}，h、w分别是图像行高与列宽；x^r∈{L,u,v}，L、u、v分别对应NDVI、BRI及NDWI经过色彩空间LUV变换后的三个色彩分量。

2.如权利要求1所述的方法，其特征在于，

所述作物群体生物量光谱指数NDVI的计算公式为 $\mathrm{NDVI}=\frac{R_{\mathrm{nir}}-R_{\mathrm{red}}}{R_{\mathrm{nir}}+R_{\mathrm{red}}},$

所述作物养分光谱指数BRI的计算公式为

所述水分敏感光谱指数NDWI的计算公式为

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

3.如权利要求1所述的方法，其特征在于，所述光谱指数SAVI的计算公式为，

$\mathrm{SAVI}=\frac{(R_{\mathrm{nir}}-R_{\mathrm{red}})}{(R_{\mathrm{nir}}+R_{\mathrm{red}}+L)}(1+L);$

其中，R_red为红波段，其波长范围为0.63-0.69微米；R_nir为近红外波段，其波长范围为0.76-0.90微米；L为对应NDVI经过色彩空间LUV变换后的色彩分量。

4.如权利要求1所述的方法，其特征在于，所述纹理结构计算包括熵值计算和同质度计算；

所述熵值计算的公式为：熵值 $\mathrm{Ent}=-\underset{i=0}{\overset{N}{\mathrm{\&Sigma;}}}\underset{j=0}{\overset{N}{\mathrm{\&Sigma;}}}{C}_{\mathrm{ij}}\log \frac{P_{\mathrm{ij}}(d,\mathrm{\&theta;})}{\underset{i=0}{\overset{N}{\mathrm{\&Sigma;}}}\underset{j=0}{\overset{N}{\mathrm{\&Sigma;}}}{P}_{\mathrm{ij}}(d,\mathrm{\&theta;})};$

所述同质度的计算公式为：

同质度 $\mathrm{Homo}=\underset{i=0}{\overset{N}{\mathrm{\&Sigma;}}}\underset{j=0}{\overset{N}{\mathrm{\&Sigma;}}}\frac{1}{1+{(i-j)}^2}\frac{P_{\mathrm{ij}}(d,\mathrm{\&theta;})}{\underset{i=0}{\overset{N}{\mathrm{\&Sigma;}}}\underset{j=0}{\overset{N}{\mathrm{\&Sigma;}}}{P}_{\mathrm{ij}}(d,\mathrm{\&theta;})};$

其中，P_i,j(d,θ)为利用对象内像元NDVI、BRI及NDWI计算的灰度共生矩阵，N为对象内部像元数目；i、j为自然数；d为计算灰度共生矩阵时相邻两个像素i与j的距离，θ为计算灰度共生矩阵时相邻两个像素i与j直线连接方向与直角坐标轴间的夹角。

5.如权利要求1所述的方法，其特征在于，所述回归模型为LAI＝a_m·[SAVI Ent_NDVI Ent_BRI Ent_NDWI Homo_NDVI Homo_BRI Homo_NDWI]^T

其中a_m＝[a₀ a₁ a₂ a₃ a₄ a₅ a₆]为回归系数，SAVI为与LAI敏感的光谱指数，R_blue为蓝波段,R_green为绿波段,R_red为红波，R_nir为近红外波段，Ent为熵值，Homo为同质度。