[发明专利]多分辨率遥感影像的城市植被自动提取方法

摘要

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

主权项

一种多分辨率遥感影像的城市植被自动提取方法，其特征在于所述方法主要包括以下步骤：步骤一、基于高空间分辨率的多光谱影像植被初始斑块提取；具体包括以下步骤：步骤1.1计算NDVI植被光谱指数，计算公式为NDVI=NIR-RNIR+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(x)=Σi=1np(ωi)p(x|ωi)]]>其中，n表示类别的总数量，ωi表示第i个地物类别，均值和方差分别用mi、表示，其概率密度分布函数表示为p(x|ωi)=12πσi2exp{-(x-mi)22σi2}]]>其中，p(ωi)为各个地物类别的初始概率密度，其满足的条件为：Σi=1np(ωi)=1]]>在上述条件下，求解植被与非植被的阈值可转化为估算地物类别的数量n，各个类别的ωi的均值和方差分mi、具体步骤如下：(2)统计分块影像中NDVI的频率直方图f(x)，搜索分块影像中NDVI的最大值NDVImax与最小值NDVImin，将NDVImin至NDVImax平均分为256份，统计NDVI值出现在每一份区间的像素数量，得到NDVI的直方图f(x)；(3)计算NDVI的直方图f(x)的一阶倒数，根据一阶倒数为0的位置寻找NDVI频率直方图中的极大值点、极小值点；(4)统计极大值点的数量，确定地表覆盖类别数量n，以极大值处的NDVI值作为每一个地表覆盖类别的初始均值mi，以相邻两个极小值点计算初始σi的公式如下：σi＝(NDVI2‑NDVI1)/4其中，NDVI2，NDVIi分别表示相邻两个极小值处对应的NDVI值；p(ωi)=ncountn_total]]>其中，其中ncount为NDVI值在相邻两个极小值之间的所有像素的个数，n_total为分块的总像素个数；(5)基于最大数学期望(EM)算法，采用循环迭代估算p(ωi)、mi、通过循环迭代，每次迭代由求期望值和期望最大化两个步骤组成；求期望值的步骤根据待估计参数的当前值，从观测数据中直接估计概率密度的期望值，期望最大化的步骤通过最大化这一期望来更新参数的估计量，这两步在整个迭代过程中依次交替进行，直至迭代过程收敛；具体计算公式如下：pt+1(ωk)=ΣX(i,j)∈MDpt(ωk)pt(X(i,j)/ωk)pt(X(i,j))I*I]]>mkt+1(ωk)=ΣX(i,j)∈MDpt(ωk)pt(X(i,j)/ωk)pt(X(i,j))X(i,j)ΣX(i,j)∈MDpt(ωk)pt(X(i,j)/ωk)pt(X(i,j))]]>(σk2)t+1(ωk)=ΣX(i,j)∈MDpt(ωk)pt(X(i,j)/ωk)pt(X(i,j))[X(i,j)-mkt]2ΣX(i,j)∈MDpt(ωk)pt(X(i,j)/ωk)pt(X(i,j))]]>上述三式估计的分别是先验概率、均值和标准差，式中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)、mi和σi的值小于给定的阈值ε时迭代终止，所述阈值ε＝10‑8；(6)确定分块的初始植被区域，假设上步中得到的每个类别的NDVI均值值按从大到小排序，m1，m2，mn，对应的方差为σi，则初始植被分布图可表示为VI(i，j)，其中1表示植被区域，0表示非植被区域，则VI(i,j)=1NDVI(i,j)>T0NDVI(i,j)≤T]]>其中阈值T通过如下公式求得T=m1-1.5σ1m1-m2>1.5(σ1+σ2)m2+1.5σ2m1-m2≤1.5(σ1+σ2);]]>步骤二、对高空间分辨率的多光谱影像计算视觉感知参数；步骤三、对高空间分辨率的全色波段影像计算纹理特征参数；步骤四、多特征综合的植被区域自动增长，获得植被区域分布图。