CN104851113A - 多分辨率遥感影像的城市植被自动提取方法 - Google Patents

Abstract

本发明公开了一种多分辨率遥感影像的城市植被自动提取方法，包括基于高空间分辨率的多光谱影像植被初始斑块提取、对高空间分辨率的多光谱影像计算视觉感知参数、对高空间分辨率的全色波段影像计算纹理特征参数、多特征综合的植被区域自动增长，获得植被区域分布图等步骤。采用了图像分块分割处理，加快图像处理速度。根据NDVI特征，采用最大数学期望算法自适应的动态阈值的自动选择。将遥感影像的全色波段、多光谱波段充分利用，对从全色波段与多光谱波段中提取的NDVI特征、视觉特征、纹理特征综合，对植被进行植被区域判断，提高了准确性。

Description

多分辨率遥感影像的城市植被自动提取方法

技术领域

本发明涉及一种多分辨率遥感影像的城市植被自动提取方法，属于遥感影像领域。

背景技术

城市植被具有固碳、释养功能，在保持城市生态环境方面具有重要作用。基于高分辨率遥感影像提取的城市植被覆盖，是评估生态城市、园林城市的重要指标。现阶段利用高分辨率遥感影像提取城市植被的方法主要有以下几种：

(1)人工目视解译判读。该方法主要以人工判读植被区域，然后勾绘植被边界。在目视解译中主要依据植被在图像中表现的特征包括形状、大小、颜色和色调、阴影、位置、纹理关系等。目视解译采用人工作业屏幕数字化的方法，自动化程度较低，主观性强，需要有丰富的影像解译实践经验，并对所提取的专题信息有较好的理解和认识。但是由于能充分地利用遥感影像中各种信息，仍然是遥感影像处理的重要方法。

(2)面向对象的监督分类方法。该方法主要利用植被区域在影像中的色彩、空间纹理特性，在面向对象的图像分割基础上，计算每个对象的各种特征，包括光谱(其包括：光谱均值、亮度、均方差、比率、与领域相关的光谱特征、与父对象相关的光谱特征等等)形状(包括面积、周长、边长、紧密度、长宽比以及与子对象和父对象相关的形状特征)、纹理(根据灰度共生矩阵定义出的二阶矩(能量)、对比度、相关、熵、逆差矩)及拓扑特征等，最后选择相应的特征进行基于对象的分类提取。

上述这些提取方法都存在缺陷和问题，包括：

(1)在提取植被区域时，需要事先选择样本，生成先验知识，而非基于图像自身特性。这个过程中有人机交互，人为的影响因素多。

(2)自动化程度不高。在分类过程中分割的尺度，植被分类的阈值、特征都需要人工交互选择，这些中间过程严重影响了植被提取的自动化程度。另外，在分割过程中，整幅图像采用同一个阈值，没有充分考虑不同植被类型之间的差别(如草地、针叶林、阔叶林等)在影像上的特征差异。

(3)只利用了高空间分辨率影像中的多光谱数据，没有利用具有更高空间分辨率全色波段的特点，而实际中，像WorldView、QuickBird等高空间分辨率的遥感数据全色与多光谱影像是捆绑在一起的。只利用了高空间分辨率影像中的多光谱数据，没有利用具有更高空间分辨率全色波段的特点，而实际中，像WorldView、QuickBird等高空间分辨率的遥感数据全色与多光谱影像是捆绑在一起的。

发明内容

为了解决上述问题，本发明针对IKONOS、QuickBird、WorldView等具有全色和多光谱波段的高分辨率遥感影像，全自动提取城市植被区域，提出了一种多分辨率遥感影像的城市植被自动提取方法。

本发明提供一种多分辨率遥感影像的城市植被自动提取方法，主要包括以下步骤：

步骤一、基于高空间分辨率的多光谱影像植被初始斑块提取；

步骤二、对高空间分辨率的多光谱影像计算视觉感知参数；

步骤三、对高空间分辨率的全色波段影像计算纹理特征参数；

步骤四、多特征综合的植被区域自动增长，获得植被区域分布图。

优选的，上述步骤一具体包括以下步骤：

步骤1.1计算NDVI植被光谱指数，计算公式为

$\mathrm{NDVI}=\frac{\mathrm{NIR}-R}{\mathrm{NIR}+R}$

其中NIR表示近红外波段的像素值，R表示红色波段的像素值；

步骤1.2对NDVI影像进行分块，每块大小300*300像素，在影像边缘部分不足300*300像素时，以实际的大小为准；

步骤1.3对每一个300*300像素的NDVI分块数据自动获取分割阈值区分植被与非植被，获得初始植被区域；

步骤1.4重复步骤1.3，得到所有分块的植被初始区域，合并所有分块中的植被区域，得到初始植被区域，对整个图像植被区域进行连通性标记，获得每一个植被斑块的像素坐标。

优选的，上述步骤1.3中分割阈值的获取通过以下方法获取：

(1)每一个分块NDVI中对应若干种地表覆盖类型，以ω_i表示，假设每个地物类别的NDVI分布服从高斯密度函数分布，则每一个分块中NDVI总的概率分布函数p(x)可表示为：

$p\left(x\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}p\left({\mathrm{\ω}}_i\right)p\left(x\right|{\mathrm{\ω}}_i)$

其中，n表示类别的总数量，ω_i表示第i个地物类别，均值和方差分别用m_i、表示，其概率密度分布函数表示为

$p\left(x\right|{\mathrm{\ω}}_i)=\frac{1}{\sqrt{2\mathrm{\π}{\mathrm{\σ}}_i^2}}\exp \{-\frac{{(x-m_i)}^2}{2{\mathrm{\σ}}_i^2}\}$

其中，p(ω_i)为各个地物类别的初始概率密度，其满足的条件为：

$\underset{i=1}{\overset{n}{\mathrm{\Σ}}}p\left({\mathrm{\ω}}_i\right)=1$

在上述条件下，求解植被与非植被的阈值可转化为估算地物类别的数量n，各个类别的ω_i的均值和方差分m_i、具体步骤如下：

(2)统计分块影像中NDVI的频率直方图f(x)，搜索分块影像中NDVI的最大值NDVI_max与最小值NDVI_min，将NDVI_min至NDVI_max平均分为256份，统计NDVI值出现在每一份区间的像素数量，得到NDVI的直方图f(x)；

(3)计算NDVI的直方图f(x)的一阶倒数，根据一阶倒数为0的位置寻找NDVI频率直方图中的极大值点、极小值点；

(4)统计极大值点的数量，确定地表覆盖类别数量n，以极大值处的NDVI值作为每一个地表覆盖类别的初始均值m_i，以相邻两个极小值点计算初始σ_i的公式如下：

σ_i＝(NDVI₂-NDVI₁)/4

其中，NDVI₂，NDVI₁分别表示相邻两个极小值处对应的NDVI值；

$p\left({\mathrm{\ω}}_i\right)=\frac{\mathrm{ncount}}{n\_\mathrm{total}}$

其中，其中ncount为NDVI值在相邻两个极小值之间的所有像素的个数，n_total为分块的总像素个数；

(5)基于最大数学期望(EM)算法，采用循环迭代估算p(ω_i)、m_i、EM算法通过循环迭代，每次迭代由求期望值和期望最大化两个步骤组成；前者根据待估计参数的当前值，从观测数据中直接估计概率密度的期望值，后者通过最大化这一期望来更新参数的估计量，这两步在整个迭代过程中依次交替进行，直至迭代过程收敛；具体计算公式如下：

$p^{t+1}\left({\mathrm{\ω}}_k\right)=\frac{\underset{X(i,j)\∈M_D}{\mathrm{\Σ}}\frac{p^t\left({\mathrm{\ω}}_k\right)p^t(X(i,j)/{\mathrm{\ω}}_k)}{p^t\left(X(i,j)\right)}}{I*I}$

$m_k^{t+1}\left({\mathrm{\ω}}_k\right)=\frac{\underset{X(i,j)\∈M_D}{\mathrm{\Σ}}\frac{p^t\left({\mathrm{\ω}}_k\right)p^t(X(i,j)/{\mathrm{\ω}}_k)}{p^t\left(X(i,j)\right)}X(i,j)}{\underset{X(i,j)\∈M_D}{\mathrm{\Σ}}\frac{p^t\left({\mathrm{\ω}}_k\right)p^t(X(i,j)/{\mathrm{\ω}}_k)}{p^t\left(X(i,j)\right)}}$

${\left({\mathrm{\σ}}_k^2\right)}^{t+1}\left({\mathrm{\ω}}_k\right)=\frac{\underset{X(i,j)\∈M_D}{\mathrm{\Σ}}\frac{p^t\left({\mathrm{\ω}}_k\right)p^t(X(i,j)/{\mathrm{\ω}}_k)}{p^t\left(X(i,j)\right)}{[X(i,j)-m_k^t]}^2}{\underset{X(i,j)\∈M_D}{\mathrm{\Σ}}\frac{p^t\left({\mathrm{\ω}}_k\right)p^t(X(i,j)/{\mathrm{\ω}}_k)}{p^t\left(X(i,j)\right)}}$

上述三式估计的分别是先验概率、均值和标准差，式中k代表第k个地物类别，t和t+1分别代表了当前和下一次迭代所用的估计值，i，j分别代表了NDVI影像的行数和列数，X(i，j)表示NDVI影像中第i行j列的NDVI值，条件概率p(X(i，j)|ω_i)的计算和全概率p(X(i，j))的值由步骤(1)中的相应公式得出；当相邻两次迭代计算的p(ω_i)、m_i和σ_i的值小于给定的阈值ε(ε＝10^-8)时迭代终止；

(6)确定分块的初始植被区域，假设上步中得到的每个类别的NDVI均值值按从大到小排序，m₁，m₂，m_n，对应的方差为σ_i，则初始植被分布图可表示为VI(i，j)，其中1表示植被区域，0表示非植被区域，则

$\mathrm{VI}(i,j)=\left\{\begin{array}{cc}1& \mathrm{NDVI}(i,j)>T\\ 0& \mathrm{NDVI}(i,j)\≤T\end{array}\right.$

其中阈值T通过如下公式求得

$T=\left\{\begin{array}{cc}{m}_1-1.5{\mathrm{\σ}}_1& m_1-m_2>1.5({\mathrm{\σ}}_1+{\mathrm{\σ}}_2)\\ m_2+1.5{\mathrm{\σ}}_2& m_1-m_2\≤1.5({\mathrm{\σ}}_1+{\mathrm{\σ}}_2)\end{array}\right..$

优选的，上述步骤1.4中的连通性标记采用4邻域或8邻域。

优选的，上述步骤二具体包括以下步骤：

步骤2.1对多光谱影像有RGB彩色空间变换到HSI色彩空间，通过以下公式计算植被亮度、色度、饱和度的变化结果

$\left\{\begin{array}{cc}H=\left\{\begin{array}{cc}\mathrm{\θ},& \begin{array}{cc}\mathrm{if}& G\≥B\end{array}\\ 2\mathrm{\π}-\mathrm{\θ},& \mathrm{otherwise}\end{array}\right.& \mathrm{\θ}=\arccos \frac{(R-G)+(R-B)}{2\sqrt{{(R-G)}^2+(R-B)(G-B)}}\\ S=1-\frac{3}{R+G+B}\min (R,G,B)& \\ I=(R+G+B)/3& \end{array}\right.;$

步骤2.2计算视觉特征感知参数，根据初始植被斑块作为样本，计算每一个斑块的植被的亮度、色度、饱和度的均值μ_x(H，S，I)，方差σ_x(H，S，I)，计算归一化色调饱和参数

$H^{\&prime;}(i,j)=\left\{\begin{array}{cc}1-\frac{|H(i,j)-{\mathrm{\&mu;}}_k|}{{\mathrm{\&mu;}}_k}& \\ 0,& \mathrm{otherwise}\end{array}\right.,\begin{array}{cc}\mathrm{if}& H_L<H(i,j)<H_U\end{array}$

其中，H(i，j)为色调，(i，j)代表像素坐标，μ_k第k个植被斑块的均值，其对应的方差为σ_k，H_L，H_U为植被色调的下限和上限

$\left\{\begin{array}{c}{H}_L={\mathrm{\&mu;}}_k-2{\mathrm{\&sigma;}}_k\\ H_U={\mathrm{\&mu;}}_k+2{\mathrm{\&sigma;}}_k\end{array}\right..$

优选的，上述步骤三具体为：

针对全色波段影像生成灰度共生矩阵，计算灰度共生矩阵的6个纹理参数：反差、相异性、角能量、熵，均一性、自相关系数，具体计算公式如下：

$\mathrm{CON}=\underset{i,j}{\mathrm{\&Sigma;}}{(i-j)}^{2}p(i,j)$

$\mathrm{DIS}=\underset{i,j}{\mathrm{\&Sigma;}}p(i-j)|i-j|$

$\mathrm{ASM}=\underset{i,j}{\mathrm{\&Sigma;}}{p(i-j)}^2$

$\mathrm{ENT}=\underset{i,j}{\mathrm{\&Sigma;}}p(i-j)\log \left(p(i,j)\right)$

$\mathrm{HOM}=\underset{i,j}{\mathrm{\&Sigma;}}\frac{p(i,j)}{1+{(i-j)}^2}$

$\mathrm{COR}=\underset{i,j}{\mathrm{\&Sigma;}}p(i,j)=\frac{(i-{\mathrm{\&mu;}}_i)(j-{\mathrm{\&mu;}}_j)}{{\mathrm{\&sigma;}}_x{\mathrm{\&sigma;}}_y}$

纹理特征计算时窗口的选择与全色影像与多光谱影像的分辨率有关系，假设全色波段的空间分辨率是多光谱影像的n倍，则纹理计算的窗口为2n-l，每隔n个像素计算一次。

优选的，上述步骤四具体为：以每个初始植被斑块为“种子”，采用区域增长的方法将相邻于“种子”斑块的像素附加到种子区域上。

优选的，上述区域增长的准则采用纹理特征、视觉感知特征综合采用多特征表决来确定。

优选的，上述区域增长的准则如下：

(1)对每一个初始植被斑块，分别统计亮度(I)、归一化色调(H)、饱和度(S)、反差(Contrast)、相异性(Dissimilarity)、角能量(Energy)、熵(Entropy)，均一性(Homogeneity)、自相关系数以及上述特征参数的均值、方差，假设其值为μ，σ；对初始植被斑块的8邻域像素，根据每一个特征，判断其是否为植被，判别函数定义如下

$T_2=\mathrm{\&mu;}-\sqrt{2}\mathrm{\&sigma;}$

$V(i,j)=\left\{\begin{array}{cc}1,\mathrm{if}& M(i,j)>T_2\\ 0,\mathrm{otherwise}& \end{array}\right.$

其中，M(i，j)为计算出的特征参数，在T₂表示判断阈值，i，j表示像素的行列号，1表示植被区域，0表示非植被区域，进而得到上述9个特征参数分别对应的潜在植被区域图；

(2)多特征表决确定最终植被区域，将这9幅潜在植被区域图进行投票表决，即在植被区域图上每个像素当有6个及其以上时特征都标记为1时，认为是最终的植被区域。

相对于现有技术，本发明提供的多分辨率遥感影像的城市植被自动提取方法有以下优点：

(1)本发明采用了基于最大数学期望的植被阈值自动选择方法，提高了数据处理中自动化程度。另外采用了分块处理技术，避免整幅影像采用同一个阈值判断，避免了不同植被类型(如草地、针叶林、阔叶林等)在遥感数据中表现的特征差异，增强了植被阈值的自适应性。

(2)本发明采用NDVI提出初始植被区域，在此基础上综合了高分辨率多光谱数据的视觉感知特征、全色波段的纹理特征，采用多特征投票表决的方法确定最终的植被区域。这个方法综合考虑高空间分辨率的视觉、纹理特征，提高了植被提取的准确性。

附图说明

图1为本发明流程示意图。

具体实施方式

为了便于本领域普通技术人员理解和实施本发明，下面结合附图及具体实施方式对本发明作进一步的详细描述。

本方案针对IKONOS、QuickBird、WorldView等具有全色和多光谱波段的高分辨率遥感影像，全自动提取城市植被区域。具体的步骤包括：

步骤一、基于高空间分辨率的多光谱影像植被初始斑块提取。

计算NDVI植被光谱指数。具体计算公式见(1)

$\mathrm{NDVI}=\frac{\mathrm{NIR}-R}{\mathrm{NIR}+R}---\left(1\right)$

其中NIR表示近红外波段的像素值，R表示红色波段的像素值。

对NDVI影像进行分块，每块大小300*300像素，在影像边缘部分不足300*300像素时，以实际的大小为准。

对每一个300*300像素的NDVI分块数据自动获取分割阈值区分植被与非植被，获得初始植被区域。获取分割阈值的方法如下：

由于植被在红色波段有强的吸收，而在近红外有强的反射，由红、近红外波段的计算的NDVI指数增强了植被的这种特性，能更好的反应植被与其他地物类型的差别。NDVI的值域范围[-1，1]，对于植被而言，NDVI一般都较大。

每一个分块NDVI中对应若干种地表覆盖类型(以ω_i表示)，假设每个地物类别的NDVI分布服从高斯密度函数分布，则每一个分块中NDVI总的概率分布函数p(x)可表示为

$p\left(x\right)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}p\left({\mathrm{\&omega;}}_i\right)p\left(x\right|{\mathrm{\&omega;}}_i)---\left(2\right)$

其中，n表示类别的总数量，ω_i表示第i个地物类别，均值和方差分别用m_i、表示，其概率密度分布函数表示为公式(3)，p(ω_i)为各个地物类别的初始概率密度，其满足的条件见公式(4)。

$p\left(x\right|{\mathrm{\&omega;}}_i)=\frac{1}{\sqrt{2\mathrm{\&pi;}{\mathrm{\&sigma;}}_i^2}}\exp \{-\frac{{(x-m_i)}^2}{2{\mathrm{\&sigma;}}_i^2}\}---\left(3\right)$

$\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}p\left({\mathrm{\&omega;}}_i\right)=1---\left(4\right)$

在上述条件下，求解植被与非植被的阈值可转化为估算地物类别的数量n，各个类别的ω_i的均值和方差分m_i、具体过程如下：

统计分块影像中NDVI的频率直方图f(x)。搜索分块影像中NDVI的最大值NDVI_max与最小值NDVI_min，将NDVI_min至NDVI_max平均分为256份，统计NDVI值出现在每一份区间的像素数量，得到NDVI的直方图f(x)。

计算NDVI的直方图f(x)的一阶倒数，根据一阶倒数为0的位置寻找NDVI频率直方图中的极大值点、极小值点。

统计极大值点的数量，确定地表覆盖类别数量n。以极大值处的NDVI值作为每一个地表覆盖类别的初始均值m_i，以相邻两个极小值点按公式(5)计算初始σ_i.

σ_i＝(NDVI₂-NDVI₁)/4   (5)

NDVI₂，NDVI₁分别表示相邻两个极小值处对应的NDVI值。

$p\left({\mathrm{\&omega;}}_i\right)=\frac{\mathrm{ncount}}{n\_\mathrm{total}}$

其中ncount为NDVI值在相邻两个极小值之间的所有像素的个数，n_total为分块的总像素个数。

基于最大数学期望(EM)算法，采用循环迭代估算p(ω_i)、m_i、EM算法通过循环迭代，每次迭代由求期望值和期望最大化两个步骤组成；前者根据待估计参数的当前值，从观测数据中直接估计概率密度的期望值，后者通过最大化这一期望来更新参数的估计量，这两步在整个迭代过程中依次交替进行，直至迭代过程收敛；具体计算公式如下：

$p^{t+1}\left({\mathrm{\&omega;}}_k\right)=\frac{\underset{X(i,j)\&Element;M_D}{\mathrm{\&Sigma;}}\frac{p^t\left({\mathrm{\&omega;}}_k\right)p^t(X(i,j)/{\mathrm{\&omega;}}_k)}{p^t\left(X(i,j)\right)}}{I*I}---\left(5\right)$

$m_k^{t+1}\left({\mathrm{\&omega;}}_k\right)=\frac{\underset{X(i,j)\&Element;M_D}{\mathrm{\&Sigma;}}\frac{p^t\left({\mathrm{\&omega;}}_k\right)p^t(X(i,j)/{\mathrm{\&omega;}}_k)}{p^t\left(X(i,j)\right)}X(i,j)}{\underset{X(i,j)\&Element;M_D}{\mathrm{\&Sigma;}}\frac{p^t\left({\mathrm{\&omega;}}_k\right)p^t(X(i,j)/{\mathrm{\&omega;}}_k)}{p^t\left(X(i,j)\right)}}---\left(6\right)$

${\left({\mathrm{\&sigma;}}_k^2\right)}^{t+1}\left({\mathrm{\&omega;}}_k\right)=\frac{\underset{X(i,j)\&Element;M_D}{\mathrm{\&Sigma;}}\frac{p^t\left({\mathrm{\&omega;}}_k\right)p^t(X(i,j)/{\mathrm{\&omega;}}_k)}{p^t\left(X(i,j)\right)}{[X(i,j)-m_k^t]}^2}{\underset{X(i,j)\&Element;M_D}{\mathrm{\&Sigma;}}\frac{p^t\left({\mathrm{\&omega;}}_k\right)p^t(X(i,j)/{\mathrm{\&omega;}}_k)}{p^t\left(X(i,j)\right)}}---\left(7\right)$

上面三式估计的分别是先验概率、均值和标准差，式中k代表第k个地物类别，t和t+1分别代表了当前和下一次迭代所用的估计值，i，j分别代表了NDVI影像的行数和列数，X(i，j)表示NDVI影像中第i行j列的NDVI值，条件概率p(X(i，j)|ω_i)的计算见(3)，全概率p(X(i，j))的值由式(2)给出；当相邻两次迭代计算的p(ω_i)、m_i和σ_i的值小于给定的阈值ε(ε＝10^-8)时迭代终止。

确定分块的初始植被区域。假设上步中得到的每个类别的NDVI均值值按从大到小排序，m₁，m₂，m_n，对应的方差为σ_i，则初始植被分布图可表示为VI(i，j)，其中1表示植被区域，0表示非植被区域。

则

$\mathrm{VI}(i,j)=\left\{\begin{array}{cc}1& \mathrm{NDVI}(i,j)>T\\ 0& \mathrm{NDVI}(i,j)\&le;T\end{array}\right.---\left(8\right)$

其中阈值T见公式

$T=\left\{\begin{array}{cc}{m}_1-1.5{\mathrm{\&sigma;}}_1& m_1-m_2>1.5({\mathrm{\&sigma;}}_1+{\mathrm{\&sigma;}}_2)\\ m_2+1.5{\mathrm{\&sigma;}}_2& m_1-m_2\&le;1.5({\mathrm{\&sigma;}}_1+{\mathrm{\&sigma;}}_2)\end{array}\right.---\left(9\right)$

重复1.3步，得到所有分块的植被初始区域。合并所有分块中的植被区域，得到初始植被区域。然后对整个图像植被区域进行连通性标记，获得每一个植被斑块的像素坐标。连通性标记可采用4邻域或8邻域。

步骤二、对高空间分辨率的多光谱影像计算视觉感知参数。

2.1对多光谱影像有RGB彩色空间变换到HSI色彩空间。

高分辨率多光谱影像中对应的红、绿、蓝波段的影像，按照公式(10)计算植被亮度、色度、饱和度的变化结果。

$\left\{\begin{array}{cc}H=\left\{\begin{array}{cc}\mathrm{\&theta;},& \begin{array}{cc}\mathrm{if}& G\&GreaterEqual;B\end{array}\\ 2\mathrm{\&pi;}-\mathrm{\&theta;},& \mathrm{otherwise}\end{array}\right.& \mathrm{\&theta;}=\arccos \frac{(R-G)+(R-B)}{2\sqrt{{(R-G)}^2+(R-B)(G-B)}}\\ S=1-\frac{3}{R+G+B}\min (R,G,B)& \\ I=(R+G+B)/3& \end{array}\right.---\left(10\right)$

2.2计算视觉特征感知参数。

根据初始植被斑块作为样本，计算每一个斑块的植被的亮度、色度、饱和度的均值μ_x(H，S，I)，方差σ_x(H，S，I)等统计参数值。计算归一化色调饱和参数

$H^{\&prime;}(i,j)=\left\{\begin{array}{cc}1-\frac{|H(i,j)-{\mathrm{\&mu;}}_k|}{{\mathrm{\&mu;}}_k}& \\ 0,& \mathrm{otherwise}\end{array}\right.,\begin{array}{cc}\mathrm{if}& H_L<H(i,j)<H_U---\left(11\right)\end{array}$

式中，(i，j)代表像素坐标，μ_k第k个植被斑块的均值，其对应的方差为σk，H_L，H_U为植被色调的下限和上限，见公式(11)

$\left\{\begin{array}{c}{H}_L={\mathrm{\&mu;}}_k-2{\mathrm{\&sigma;}}_k\\ H_U={\mathrm{\&mu;}}_k+2{\mathrm{\&sigma;}}_k\end{array}\right.---\left(11\right)$

步骤三、对高空间分辨率的全色波段影像计算纹理特征参数

针对全色波段影像生成灰度共生矩阵，计算灰度共生矩阵的6个纹理参数：反差(Con)、相异性(Dis)、角能量(ASM)、熵(Ent)，均一性(Hom)、自相关系数(Cor)，具体公式见(11-17)。纹理特征计算时窗口的选择与全色影像与多光谱影像的分辨率有关系，假色全色波段的空间分辨率是多光谱影像的n倍，则纹理计算的窗口为2n-1，每隔n个像素计算一次。例如quickbird影像全色波段影像的分辨率为0.6 1m，而多光谱影像的分辨率为2.4 4m，计算全色波段纹理时窗口的大小为7×7，每隔4个像素计算一次纹理特征。这样便可以将全色影像的纹理特征和多光谱影像的视觉感知特征及NDVI相对应。

$\mathrm{CON}=\underset{i,j}{\mathrm{\&Sigma;}}{(i-j)}^{2}p(i,j)---\left(11\right)$

$\mathrm{DIS}=\underset{i,j}{\mathrm{\&Sigma;}}p(i-j)|i-j|---\left(12\right)$

$\mathrm{ASM}=\underset{i,j}{\mathrm{\&Sigma;}}{p(i-j)}^2---\left(13\right)$

$\mathrm{ENT}=\underset{i,j}{\mathrm{\&Sigma;}}p(i-j)\log \left(p(i,j)\right)---\left(14\right)$

$\mathrm{HOM}=\underset{i,j}{\mathrm{\&Sigma;}}\frac{p(i,j)}{1+{(i-j)}^2}---\left(15\right)$

$\mathrm{COR}=\underset{i,j}{\mathrm{\&Sigma;}}p(i,j)=\frac{(i-{\mathrm{\&mu;}}_i)(j-{\mathrm{\&mu;}}_j)}{{\mathrm{\&sigma;}}_x{\mathrm{\&sigma;}}_y}---\left(16\right)$

步骤四、多特征综合的植被区域自动增长，获得植被区域分布图

以每个初始植被斑块为“种子”，采用区域增长的方法将相邻于“种子”斑块的像素附加到种子区域上。区域增长的准则采用纹理特征、视觉感知特征综合采用多特征表决来确定。

区域增长的准则如下

对每一个初始植被斑块，分别统计亮度(I)、归一化色调(H)、饱和度(S)、反差(Contrast)、相异性(Dissimilarity)、角能量(Energy)、熵(Entropy)，均一性(Homogeneity)、自相关系数，这9个特征参数的均值、方差，假设其值为μ，σ。对初始植被斑块的8邻域像素，根据每一个特征，判断其是否为植被，判别函数定义如下

$T_2=\mathrm{\&mu;}-\sqrt{2}\mathrm{\&sigma;}---\left(17\right)$

$V(i,j)=\left\{\begin{array}{cc}1,\mathrm{if}& M(i,j)>T_2\\ 0,\mathrm{otherwise}& \end{array}\right.---\left(18\right)$

其中，M(i，j)是指计算出的特征参数，就是(1)中所述的3个视觉特征参数，6个纹理特征参数，共有9个，即亮度、归一化色调、饱和度、反差、相异性、角能量、熵，均一性、自相关系数。此处用M来代替其中的任意一个变量。在T2表示判断阈值，i，j表示像素的行列号，1表示植被区域，0表示非植被区域。

这样，可以亮度、归一化色调、饱和度、反差(Contrast)、相异性(Dissimilarity)、角能量(Energy)、熵(Entropy)，均一性(Homogeneity)、自相关系数，这9个特征参数可以分别得到一幅潜在植被区域图。

多特征表决确定最终植被区域。将这9幅潜在植被区域图进行投票表决，即在植被区域图上每个像素当有6个及其以上特征都标记为1时，认为是最终的植被区域。

本发明提供的多分辨率遥感影像的城市植被自动提取方法，采用了图像分块分割处理，加快图像处理速度。根据NDVI特征，采用最大数学期望算法自适应的动态阈值的自动选择。将遥感影像的全色波段、多光谱波段充分利用，对从全色波段与多光谱波段中提取的NDVI特征、视觉特征、纹理特征综合，对植被进行植被区域判断，提高了准确性。本发明提供的方法还具有一定的适用性，如本方法中的NDVI植被指数可以换成EVI指数、SR指数等。

以上所述，仅是用以说明本发明的具体实施案例而已，并非用以限定本发明的可实施范围，举凡本领域熟练技术人员在未脱离本发明所指示的精神与原理下所完成的一切等效改变或修饰，仍应由本发明权利要求的范围所覆盖。

Claims (9)

1.一种多分辨率遥感影像的城市植被自动提取方法，其特征在于所述方法主要包括以下步骤：

步骤一、基于高空间分辨率的多光谱影像植被初始斑块提取；

步骤二、对高空间分辨率的多光谱影像计算视觉感知参数；

步骤三、对高空间分辨率的全色波段影像计算纹理特征参数；

步骤四、多特征综合的植被区域自动增长，获得植被区域分布图。

2.根据权利要求1所述的多分辨率遥感影像的城市植被自动提取方法，其特征在于：所述步骤一具体包括以下步骤：

步骤1.1计算NDVI植被光谱指数，计算公式为

$\mathrm{NDVI}=\frac{\mathrm{NIR}-R}{\mathrm{NIR}+R}$

其中NIR表示近红外波段的像素值，R表示红色波段的像素值；

步骤1.2对NDVI影像进行分块，每块大小300*300像素，在影像边缘部分不足300*300像素时，以实际的大小为准；

步骤1.3对每一个300*300像素的NDVI分块数据自动获取分割阈值区分植被与非植被，获得初始植被区域；

步骤1.4重复步骤1.3，得到所有分块的植被初始区域，合并所有分块中的植被区域，得到初始植被区域，对整个图像植被区域进行连通性标记，获得每一个植被斑块的像素坐标。

3.根据权利要求2所述的多分辨率遥感影像的城市植被自动提取方法，其特征在于：所述步骤1.3中分割阈值的获取通过以下方法获取：

(1)每一个分块NDVI中对应若干种地表覆盖类型，以ω_i表示，假设每个地物类别的NDVI分布服从高斯密度函数分布，则每一个分块中NDVI总的概率分布函数p(x)可表示为：

$p\left(x\right)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}p\left({\mathrm{\&omega;}}_i\right)p\left(x\right|{\mathrm{\&omega;}}_i)$

其中，n表示类别的总数量，ω_i表示第i个地物类别，均值和方差分别用m_i、表示，其概率密度分布函数表示为

$p\left(x\right|{\mathrm{\&omega;}}_i)=\frac{1}{\sqrt{2\mathrm{\&pi;}{\mathrm{\&sigma;}}_i^2}}\exp \{-\frac{{(x-m_i)}^2}{2{\mathrm{\&sigma;}}_i^2}\}$

其中，p(ω_i)为各个地物类别的初始概率密度，其满足的条件为：

$\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}p\left({\mathrm{\&omega;}}_i\right)=1$

在上述条件下，求解植被与非植被的阈值可转化为估算地物类别的数量n，各个类别的ω_i的均值和方差分m_i、具体步骤如下：

(2)统计分块影像中NDVI的频率直方图f(x)，搜索分块影像中NDVI的最大值NDVI_max与最小值NDVI_min，将NDVI_min至NDVI_max平均分为256份，统计NDVI值出现在每一份区间的像素数量，得到NDVI的直方图f(x)；

(3)计算NDVI的直方图f(x)的一阶倒数，根据一阶倒数为0的位置寻找NDVI频率直方图中的极大值点、极小值点；

(4)统计极大值点的数量，确定地表覆盖类别数量n，以极大值处的NDVI值作为每一个地表覆盖类别的初始均值m_i，以相邻两个极小值点计算初始σ_i的公式如下：

σ_i＝(NDVI₂-NDVI₁)/4

其中，NDVI₂，NDVI₁分别表示相邻两个极小值处对应的NDVI值；

$p\left({\mathrm{\&omega;}}_i\right)=\frac{\mathrm{ncount}}{n\_\mathrm{total}}$

其中，其中ncount为NDVI值在相邻两个极小值之间的所有像素的个数，n_total为分块的总像素个数；

(5)基于最大数学期望(EM)算法，采用循环迭代估算p(ω_i)、m_i、通过循环迭代，每次迭代由求期望值和期望最大化两个步骤组成；求期望值的步骤根据待估计参数的当前值，从观测数据中直接估计概率密度的期望值，期望最大化的步骤通过最大化这一期望来更新参数的估计量，这两步在整个迭代过程中依次交替进行，直至迭代过程收敛；具体计算公式如下：

$p^{t+1}\left({\mathrm{\&omega;}}_k\right)=\frac{\underset{X(i,j)\&Element;M_D}{\mathrm{\&Sigma;}}\frac{p^t\left({\mathrm{\&omega;}}_k\right)p^t(X(i,j)/{\mathrm{\&omega;}}_k)}{p^t\left(X(i,j)\right)}}{I*I}$

$m_k^{t+1}\left({\mathrm{\&omega;}}_k\right)=\frac{\underset{X(i,j)\&Element;M_D}{\mathrm{\&Sigma;}}\frac{p^t\left({\mathrm{\&omega;}}_k\right)p^t(X(i,j)/{\mathrm{\&omega;}}_k)}{p^t\left(X(i,j)\right)}X(i,j)}{\underset{X(i,j)\&Element;M_D}{\mathrm{\&Sigma;}}\frac{p^t\left({\mathrm{\&omega;}}_k\right)p^t(X(i,j)/{\mathrm{\&omega;}}_k)}{p^t\left(X(i,j)\right)}}$

${\left({\mathrm{\&sigma;}}_k^2\right)}^{t+1}\left({\mathrm{\&omega;}}_k\right)=\frac{\underset{X(i,j)\&Element;M_D}{\mathrm{\&Sigma;}}\frac{p^t\left({\mathrm{\&omega;}}_k\right)p^t(X(i,j)/{\mathrm{\&omega;}}_k)}{p^t\left(X(i,j)\right)}{[X(i,j)-m_k^t]}^2}{\underset{X(i,j)\&Element;M_D}{\mathrm{\&Sigma;}}\frac{p^t\left({\mathrm{\&omega;}}_k\right)p^t(X(i,j)/{\mathrm{\&omega;}}_k)}{p^t\left(X(i,j)\right)}}$

上述三式估计的分别是先验概率、均值和标准差，式中k代表第k个地物类别，t和t+1分别代表了当前和下一次迭代所用的估计值，i，j分别代表了NDVI影像的行数和列数，X(i，j)表示NDVI影像中第i行j列的NDVI值，条件概率p(X(i，j)|ω_i)的计算和全概率p(X(i，j))的值由步骤(1)中的相应公式得出；当相邻两次迭代计算的p(ω_i)、m_i和σ_i的值小于给定的阈值ε(ε＝10^-8)时迭代终止；

(6)确定分块的初始植被区域，假设上步中得到的每个类别的NDVI均值值按从大到小排序，m₁，m₂，m_n，对应的方差为σ_i，则初始植被分布图可表示为VI(i，j)，其中1表示植被区域，0表示非植被区域，则

$\mathrm{VI}(i,j)=\left\{\begin{array}{cc}1& \mathrm{NDVI}(i,j)>T\\ 0& \mathrm{NDVI}(i,j)\&le;T\end{array}\right.$

其中阈值T通过如下公式求得

$T=\left\{\begin{array}{cc}{m}_1-1.5{\mathrm{\&sigma;}}_1& m_1-m_2>1.5({\mathrm{\&sigma;}}_1+{\mathrm{\&sigma;}}_2)\\ m_2+1.5{\mathrm{\&sigma;}}_2& m_1-m_2\&le;1.5({\mathrm{\&sigma;}}_1+{\mathrm{\&sigma;}}_2)\end{array}\right..$

4.根据权利要求2所述的多分辨率遥感影像的城市植被自动提取方法，其特征在于：所述步骤1.4中的连通性标记采用4邻域或8邻域。

5.根据权利要求1所述的多分辨率遥感影像的城市植被自动提取方法，其特征在于：所述步骤二具体包括以下步骤：

步骤2.1对多光谱影像有RGB彩色空间变换到HSI色彩空间，通过以下公式计算植被亮度、色度、饱和度的变化结果

$\left\{\begin{array}{c}H=\left\{\begin{array}{cc}\mathrm{\&theta;},& \mathrm{ifG}\&GreaterEqual;B\\ 2\mathrm{\&pi;}-\mathrm{\&theta;},& \mathrm{otherwise}\end{array}\right.\mathrm{\&theta;}=\arccos \frac{(R-G)+(R-B)}{2\sqrt{{(R-G)}^2+(R-B)(G-B)}}\\ S=1-\frac{3}{R+G+B}\min (R,G,B)\\ I=(R+G+B)/3\end{array}\right.;$

步骤2.2计算视觉特征感知参数，根据初始植被斑块作为样本，计算每一个斑块的植被的亮度、色度、饱和度的均值μ_x(H，S，I)，方差σ_x(H，S，I)，计算归一化色调饱和参数

$H^{\&prime;}(i,j)=\left\{\begin{array}{cc}1-\frac{|H(i,j)-{\mathrm{\&mu;}}_k|}{{\mathrm{\&mu;}}_k},& \mathrm{if}{H}_L<H(i,j)<H_U\\ 0,& \mathrm{otherwise}\end{array}\right.$

其中，H(i，j)为色调，(i，j)代表像素坐标，μ_k第k个植被斑块的均值，其对应的方差为σ_k，H_L，H_U为植被色调的下限和上限

$\left\{\begin{array}{c}{H}_L={\mathrm{\&mu;}}_k-2{\mathrm{\&sigma;}}_k\\ H_U={\mathrm{\&mu;}}_k+2{\mathrm{\&sigma;}}_k\end{array}\right..$

6.根据权利要求1所述的多分辨率遥感影像的城市植被自动提取方法，其特征在于所述步骤三具体为：

针对全色波段影像生成灰度共生矩阵，计算灰度共生矩阵的6个纹理参数：反差、相异性、角能量、熵，均一性、自相关系数，具体计算公式如下：

$\mathrm{CON}=\underset{i,j}{\mathrm{\&Sigma;}}{(i-j)}^{2}p(i,j)$

$\mathrm{DIS}=\underset{i,j}{\mathrm{\&Sigma;}}p(i,j)|i-j|$

$\mathrm{ASM}=\underset{i,j}{\mathrm{\&Sigma;}}p{(i,j)}^2$

$\mathrm{ENT}=\underset{i,j}{\mathrm{\&Sigma;}}p(i,j)\log \left(p(i,j)\right)$

$\mathrm{HOM}=\underset{i,j}{\mathrm{\&Sigma;}}\frac{p(i,j)}{1+{(i-j)}^2}$

$\mathrm{COR}=\underset{i,j}{\mathrm{\&Sigma;}}p(i,j)\frac{(i-{\mathrm{\&mu;}}_i)(j-{\mathrm{\&mu;}}_j)}{{\mathrm{\&sigma;}}_x{\mathrm{\&sigma;}}_y}$

纹理特征计算时窗口的选择与全色影像与多光谱影像的分辨率有关系，假设全色波段的空间分辨率是多光谱影像的n倍，则纹理计算的窗口为2n-1，每隔n个像素计算一次。

7.根据权利要求1所述的多分辨率遥感影像的城市植被自动提取方法，其特征在于所述步骤四具体为：以每个初始植被斑块为“种子”，采用区域增长的方法将相邻于“种子”斑块的像素附加到种子区域上。

8.根据权利要求7所述的多分辨率遥感影像的城市植被自动提取方法，其特征在于所述区域增长的准则采用纹理特征、视觉感知特征综合采用多特征表决来确定。

9.根据权利要求8所述的多分辨率遥感影像的城市植被自动提取方法，其特征在于所述区域增长的准则如下：

(1)对每一个初始植被斑块，分别统计亮度、归一化色调、饱和度、反差、相异性、角能量、熵，均一性、自相关系数以及上述特征参数的均值、方差，假设其值为μ，σ；对初始植被斑块的8邻域像素，根据每一个特征，判断其是否为植被，判别函数定义如下

$T_2=\mathrm{\&mu;}-\sqrt{2}\mathrm{\&sigma;}$

$V(i,j)=\left\{\begin{array}{cc}1,& \mathrm{ifM}(i,j)>T_2\\ 0,& \mathrm{otherwise}\end{array}\right.$

其中，M(i，j)为计算出的特征参数，在T₂表示判断阈值，i，j表示像素的行列号，1表示植被区域，0表示非植被区域，进而得到上述9个特征参数分别对应的潜在植被区域图；

(2)多特征表决确定最终植被区域，将这9幅潜在植被区域图进行投票表决，即在植被区域图上每个像素当有6个及其以上时特征都标记为1时，认为是最终的植被区域。