CN105701483A - 一种融合多光谱遥感数据和夜间灯光遥感数据的城市边界提取方法 - Google Patents

Abstract

本发明公开了一种融合多光谱遥感数据和夜间灯光遥感数据的城市边界提取方法，通过联合使用不透水面指数和夜间灯光遥感数据中提取的灯光指数可以更为精确地定位城区边界，降低误提取现象。同时，通过融合较夜间灯光遥感数据空间分辨率更高的多光谱遥感数据，将两者数据特征同时作为城市区域和非城市区域的区别特征，保证了提取的城区边界较单独使用夜间灯光遥感数据更为精细。大量实验结果表明，本发明城市边界正确提取率达到了90％以上，较现有技术提高了近10％，误提取率低至10％以下，解决了现有技术存在的问题，适用于土地规划、地理国情监测、城市扩张分析等城市遥感应用。

Description

一种融合多光谱遥感数据和夜间灯光遥感数据的城市边界提取方法

技术领域

本发明属于城市规划领域，特别涉及一种融合多光谱遥感数据和夜间灯光遥感数据的城市边界提取方法。

背景技术

城市边界自动提取是实现及时掌握城市扩展范围的关键，而城市扩展范围已成为城市规划领域内的一个重点关注问题。同时，快速自动的提取城市边界能够减少传统监测方法的人力物力。然而，当遥感影像的分辨率不够高时，自动提取的方法容易将一些地物误划入城市边界，因此，中低分辨率遥感影像城市边界提取更具有挑战性。目前较新且效果较好的城市边界提取方法以借助不透水面或夜间灯光遥感影像提取方法为主。这些方法的主要思路是：通过选取特征考虑，利用分类器(多数为支持向量机)训练城市区域检测模型并检测城市区域，从而提取出城市边界。其中从特征考虑，主要有利用多光谱遥感影像光谱特征和夜间灯光遥感影像(参考文献：XiaoP,WangX,FengX,etal.DetectingChina’sUrbanExpansionOverthePastThreeDecadesUsingNighttimeLightData[J].IEEEJournalofSelectedTopicsinAppliedEarthObservations&RemoteSensing,2014,7(10):4095-4106.)两种。由于仅利用光谱特征，很容易将非城市区域的不透水面误提取出来，从而影响城市边界提取，而仅利用夜间灯光遥感影像，由于其分辨率较低(1km)，同样会将非城市区域误提取出来，导致上述两种方法自动提取城市边界的效果不理想。

发明内容

本发明提出了一种融合多光谱遥感数据和夜间灯光遥感数据的城市边界提取方法，其目的在于，解决现有技术未能较快速准确提取城市边界这一问题，通过融合多种数据，利用多种特征，以达到快速自动提取城市边界，并提高城市边界的准确率，降低误提取率的效果。

一种融合多光谱遥感数据和夜间灯光遥感数据的城市边界提取方法，包括以下几个步骤：

步骤1：影像预处理；

对同一场景的多光谱遥感影像和夜间灯光遥感影像进行坐标转换，统一转为UTM投影坐标系；

对坐标转换后的夜间灯光遥感影像进行重采样至与坐标转换后的多光谱遥感影像的空间分辨率相同，获得Landsat图像数据I₀；

采用一致的空间分辨率30米，实现多源数据的配准；

以使得灯光遥感数据同多光谱遥感影像数据的分辨率一致，便于正负样本提取；

步骤2：自动提取正负训练样本；

首先利用Landsat图像数据I₀提供的光谱信息计算不透水面指数得到初始不透水面图像I₁，然后结合坐标转换后的夜间灯光遥感数据提供的初始灯光指数图像I₂，自动提取正负分类训练样本；

所述正负样本的自动提取过程如下：

A1：对初始不透水面图像I₁进行Otsu阈值分割，得到二值分割后的二值化不透水面图像I₃；

$I_3(i,j)=\left\{\begin{array}{cc}1& if\left(I_1\right(i,j)\&GreaterEqual;{\mathrm{th}}_1\\ 0& if\left(I_1\right(i,j)<{\mathrm{th}}_1\end{array}\right.$

其中，(i，j)表示图像中的像素坐标，th₁为Otsu阈值分割方法自适应选取的分割阈值；

A2：对初始灯光指数图像I₂进行Otsu阈值分割，得到二值分割后的二值化灯光指数图像I₄；

$I_4(i,j)=\left\{\begin{array}{cc}1& if\left(I_1\right(i,j)\&GreaterEqual;{\mathrm{th}}_2\\ 0& if\left(I_1\right(i,j)<{\mathrm{th}}_2\end{array}\right.$

其中，(i，j)表示图像中的像素坐标，th₂为Otsu阈值分割方法自适应选取的分割阈值；

A3：计算二值化灯光图像I₄中像素值为1的所有像素在初始灯光指数图像I₂中对应像素的像素值的平均值th₃，以及像素值为0的所有像素在初始灯光指数图像I₂中对应像素的像素值的平均值th₄；

A4：对二值化不透水面图像I₃进行格网划分，然后计算各个格网单位所包含像素对应的不透水面指数平均值a₁和灯光指数平均值a₂：

如果a₁≥th₁且a₂≥th₃，则将该网格中的所有像素提取至候选正样本集合P₁；

如果a₂≤th₄，则将该网格中的所有像素提取至候选负样本集合N₁；

考虑了两种数据的各自特征的优势互补，实现自动提取正负样本；

A5：从候选正样本集合P₁和候选负样本集合N₁中分别随机提取50％的正负样本构造最终的正样本集合P₂和负样本集合N₂；

步骤3：利用步骤2获得的正负样本进行城市边界提取；

步骤3.1：首先提取用于描述正负样本的分类特征F＝P₂(f)+N₂(f),f＝{f₁,f₂,f₃,f₄}；

其中，f为影像特征集合，f₁为从Landsat图像数据I₀中提取的像素的7维光谱特征，f₂为从Landsat图像数据I₀中提取的像素的2维纹理特征，f₃为从初始不透水面图像I₁中提取的像素的1维的不透水指数特征，f₄为从初始灯光指数图像I₂中提取的像素的1维的灯光指数特征，P₂为正样本集合中所有像素对应影像中的位置，N₂为负样本集合中所有像素对应影像中的位置；

然后利用正负样本的分类特征训练SVM分类器，再用训练好的分类器对待分类的场景的影像特征集合f进行分类，得到城市区域和非城市区域；

步骤3.2：利用获得的城市区域，进行分类后处理，确定出城市边界；

所述分类后处理包括对城市区域填充内部空洞以及剔除小图斑操作，其中，剔除图斑的阈值为3×3像素。

分类是以像素为单元的影像特征集合进行处理的，故需要对分类后的像素进行处理才能得到连续的闭合边界；

所述灯光遥感数据的重采样处理是指重采样方法采用三次卷积内插法。

所述利用多光谱遥感数据提供的光谱信息计算不透水面指数得到初始不透水面图像I₁的步骤如下：

步骤2.1：遥感影像辐射校正；

利用ENVI软件，对多光谱遥感影像进行辐射校正，将影像上的DN值转换为反射率；

步骤2.2：不透水面指数计算；

采用归一化差值不透水面指数NDISI来提取步骤2中初始不透水面I₁，该指数计算公式为：

$NDISI=\frac{TIR-(MNDWI+NIR+{\mathrm{MIR}}_1)/3}{TIR+(MNDWI+NIR+{\mathrm{MIR}}_1)/3}$

其中，NIR、MIR₁和TIR分别为近红外、中红外1和热红外波段的反射率，MNDWI为归一化水体指数，计算公式为：

$MNDWI=\frac{Green-{\mathrm{MIR}}_1}{Green+{\mathrm{MIR}}_1}$

其中，Green为绿波段的反射率。

通过上述步骤计算得到多光谱遥感影像中对应的各个像素的不透水面指数，从而得到初始不透水面图像I₁。

所述对二值化不透水面图像I₃进行格网划分时，每个格网单元大小为8*8。

所述f₂为从Landsat图像数据I₀中提取的像素的2维纹理特征是指在3*3窗口内的均值和方差。

有益效果

本发明提出了一种融合多光谱遥感数据和夜间灯光遥感数据的城市边界提取方法，相对于传统方法单独使用不透水面指数提取城区，通过联合使用不透水面指数和夜间灯光遥感数据中提取的灯光指数可以更为精确地定位城区边界，降低误提取现象。同时，通过融合较夜间灯光遥感数据(1km)空间分辨率更高的多光谱遥感数据(30m)，将两者数据特征同时作为城市区域和非城市区域的区别特征，对近红外、中红外1、热红外波段和MNDWI进行归一化处理，可以保证提取的城区边界较单独使用夜间灯光遥感数据更为精细。大量实验结果表明，本发明城市边界正确提取率达到了90％以上，较现有技术提高了近10％，误提取率低至10％以下，解决了现有技术存在的问题，适用于土地规划、地理国情监测、城市扩张分析等城市遥感应用。

附图说明

图1是本发明的流程图；

图2是实施例一应用本发明所述方法进行城区边界提取的遥感影像过程图，其中，(a)为原始多光谱遥感影像，(b)为原始夜间灯光遥感影像，(c)为预处理后的夜间灯光遥感影像，(d)为不透水面指数图，(e)为正样本，(f)为负样本，(g)为本发明所述的方法所得的城区边界结果示意图，(h)为现有技术的检测结果；

图3是实施例二应用本发明所述方法进行城区边界提取的遥感影像过程图，其中，(a)为原始多光谱遥感影像，(b)为原始夜间灯光遥感影像，(c)为预处理后的夜间灯光遥感影像，(d)为不透水面指数图，(e)为正样本，(f)为负样本，(g)为本发明所述的方法所得的城区边界结果示意图，(h)为现有技术的检测结果。

具体实施方式

下面将结合附图和实施例对本发明做进一步的说明。

如图1所示，一种融合多光谱遥感数据和夜间灯光遥感数据的城市边界提取方法，包括以下几个步骤：

步骤1：影像预处理；

对同一场景的多光谱遥感影像和夜间灯光遥感影像进行坐标转换，统一转为UTM投影坐标系；

对坐标转换后的夜间灯光遥感影像进行重采样至与坐标转换后的多光谱遥感影像的空间分辨率相同，获得Landsat图像数据I₀；

采用一致的空间分辨率30米，实现多源数据的配准；

以使得灯光遥感数据同多光谱遥感影像数据的分辨率一致，便于正负样本提取；

步骤2：自动提取正负训练样本；

首先利用Landsat图像数据I₀提供的光谱信息计算不透水面指数得到初始不透水面图像I₁，然后结合坐标转换后的夜间灯光遥感数据提供的初始灯光指数图像I₂，自动提取正负分类训练样本；

所述利用多光谱遥感数据提供的光谱信息计算不透水面指数得到初始不透水面图像I₁的步骤如下：

步骤2.1：遥感影像辐射校正；

利用ENVI软件，对多光谱遥感影像进行辐射校正，将影像上的DN值转换为反射率；

步骤2.2：不透水面指数计算；

采用归一化差值不透水面指数NDISI来提取步骤2中初始不透水面I₁，该指数计算公式为：

$NDISI=\frac{TIR-(MNDWI+NIR+{\mathrm{MIR}}_1)/3}{TIR+(MNDWI+NIR+{\mathrm{MIR}}_1)/3}$

其中，NIR、MIR₁和TIR分别为近红外、中红外1和热红外波段的反射率，MNDWI为归一化水体指数，计算公式为：

$MNDWI=\frac{Green-{\mathrm{MIR}}_1}{Green+{\mathrm{MIR}}_1}$

其中，Green为绿波段的反射率。

通过上述步骤计算得到多光谱遥感影像中对应的各个像素的不透水面指数，从而得到初始不透水面图像I₁。

所述正负样本的自动提取过程如下：

A1：对初始不透水面图像I₁进行Otsu阈值分割，得到二值分割后的二值化不透水面图像I₃；

$I_3(i,j)=\left\{\begin{array}{cc}1& if\left(I_1\right(i,j)\&GreaterEqual;{\mathrm{th}}_1\\ 0& if\left(I_1\right(i,j)<{\mathrm{th}}_1\end{array}\right.$

其中，(i，j)表示图像中的像素坐标，th₁为Otsu阈值分割方法自适应选取的分割阈值；

A2：对初始灯光指数图像I₂进行Otsu阈值分割，得到二值分割后的二值化灯光指数图像I₄；

$I_4(i,j)=\left\{\begin{array}{cc}1& if\left(I_1\right(i,j)\&GreaterEqual;{\mathrm{th}}_2\\ 0& if\left(I_1\right(i,j)<{\mathrm{th}}_2\end{array}\right.$

其中，(i，j)表示图像中的像素坐标，th₂为Otsu阈值分割方法自适应选取的分割阈值；

A3：计算二值化灯光图像I₄中像素值为1的所有像素在初始灯光指数图像I₂中对应像素的像素值的平均值th₃，以及像素值为0的所有像素在初始灯光指数图像I₂中对应像素的像素值的平均值th₄；

A4：对二值化不透水面图像I₃进行格网划分，然后计算各个格网单位所包含像素对应的不透水面指数平均值a₁和灯光指数平均值a₂：

如果a₁≥th₁且a₂≥th₃，则将该网格中的所有像素提取至候选正样本集合P₁；

如果a₂≤th₄，则将该网格中的所有像素提取至候选负样本集合N₁；

考虑了两种数据的各自特征的优势互补，实现自动提取正负样本；

A5：从候选正样本集合P₁和候选负样本集合N₁中分别随机提取50％的正负样本构造最终的正样本集合P₂和负样本集合N₂；

步骤3：利用步骤2获得的正负样本进行城市边界提取；

步骤3.1：首先提取用于描述正负样本的分类特征F＝P₂(f)+N₂(f),f＝{f₁,f₂,f₃,f₄}；

其中，f为影像特征集合，f₁为从Landsat图像数据I₀中提取的像素的7维光谱特征，f₂为从Landsat图像数据I₀中提取的像素的2维纹理特征，f₃为从初始不透水面图像I₁中提取的像素的1维的不透水指数特征，f₄为从初始灯光指数图像I₂中提取的像素的1维的灯光指数特征，P₂为正样本集合中所有像素对应影像中的位置，N₂为负样本集合中所有像素对应影像中的位置；

然后利用正负样本的分类特征训练SVM分类器，再用训练好的分类器对待分类的场景的影像特征集合f进行分类，得到城市区域和非城市区域；

步骤3.2：利用获得的城市区域，进行分类后处理，确定出城市边界；

所述分类后处理包括对城市区域填充内部空洞以及剔除小图斑操作，其中，剔除图斑的阈值为3×3像素。

分类是以像素为单元的影像特征集合进行处理的，故需要对分类后的像素进行处理才能得到连续的闭合边界；

所述灯光遥感数据的重采样处理是指重采样方法采用三次卷积内插法。

所述对二值化不透水面图像I₃进行格网划分时，每个格网单元大小为8*8。

所述f₂为从Landsat图像数据I₀中提取的像素的2维纹理特征是指在3*3窗口内的均值和方差。

实施例1：过程示意图见图2。

(1)影像预处理

对图2(a)和(b)，利用ENVI软件将多光谱遥感影像和夜间灯光遥感影像进行坐标转换，统一转换为UTM投影坐标系，并将夜间灯光遥感影像重采样至30米，如图2(c)所示。

(2)不透水面提取

多光谱遥感影像经步骤1处理后，可进行不透水面提取，提取过程如下：

①利用ENVI软件，对多光谱遥感影像进行辐射校正，将影像上的DN值转换为反射率；

②计算归一化水体指数MNDWI，计算公式为：

$MNDWI=\frac{Green-{\mathrm{MIR}}_1}{Green+{\mathrm{MIR}}_1}$

其中，Green和MIR₁为绿波段和中红外1的反射率。

③对近红外、中红外1、热红外波段和MNDWI进行归一化处理。

④采用归一化差值不透水面指数NDISI来提取不透水面，该指数计算公式为：

$NDISI=\frac{TIR-(MNDWI+NIR+{\mathrm{MIR}}_1)/3}{TIR+(MNDWI+NIR+{\mathrm{MIR}}_1)/3}$

其中，NIR和TIR分别为近红外和热红外波段的反射率。

提取结果如图2(d)所示。

(3)正负样本自动提取

利用步骤1处理的图2(c)和步骤2处理的图2(d)，进行正负样本的提取，具体提取过程如下：

①对初始不透水面图像I₁(图2(d))进行Otsu阈值分割，得到二值分割后的不透水面图像I₃；

$I_3(i,j)=\left\{\begin{array}{cc}1& if\left(I_1\right(i,j)\&GreaterEqual;{\mathrm{th}}_1\\ 0& if\left(I_1\right(i,j)<{\mathrm{th}}_1\end{array}\right.$

其中th₁为Otsu阈值分割方法自适应选取的分割阈值；

②对初始灯光指数图像I₂(图2(c))进行Otsu阈值分割，得到二值分割后的灯光指数图像I₄；

$I_4(i,j)=\left\{\begin{array}{cc}1& if\left(I_1\right(i,j)\&GreaterEqual;{\mathrm{th}}_2\\ 0& if\left(I_1\right(i,j)<{\mathrm{th}}_2\end{array}\right.$

其中th₂为Otsu阈值分割方法自适应选取的分割阈值；

③计算二值分割后的灯光图像I₄中像素值为1的所有像素在灯光指数图像I₂中对应像素的像素值的平均值th₃，以及像素值为0的所有像素在灯光指数图像I₂中对应像素的像素值的平均值th₄；

④对二值化不透水面图像I₃进行格网划分(每个格网单元大小为8*8)，然后计算各个格网单位所包含像素对应的不透水面指数平均值a₁和灯光指数平均值a₂，如果a₁≥th₁且a₂≥th₃，则将该网格中的所有像素提取至正样本集合P。如果a₂≤th₄，则将该网格中的所有像素提取至负样本集合N。

⑤从候选正样本集合P₁和候选负样本集合N₁中分别随机提取50％的正负样本构造最终的正样本集合P₂和负样本集合N₂；

提取正样本结果如图2(e)所示，提取负样本结果如图2(f)所示。

(4)城区边界提取

①首先提取用于描述正负样本的分类特征F＝P₂(f)+N₂(f),f＝{f₁,f₂,f₃,f₄}；

其中，f为影像特征集合，f₁为从Landsat图像数据I₀中提取的像素的7维光谱特征，f₂为从Landsat图像数据I₀中提取的像素的2维纹理特征，f₃为从初始不透水面图像I₁中提取的像素的1维的不透水指数特征，f₄为从初始灯光指数图像I₂中提取的像素的1维的灯光指数特征，P₂为正样本集合中所有像素对应影像中的位置，N₂为负样本集合中所有像素对应影像中的位置；

然后利用正负样本的分类特征训练SVM分类器，再用训练好的分类器对待分类的场景的影像特征集合f进行分类，得到城市区域和非城市区域；

②利用获得的城市区域，进行分类后处理，填充内部空洞，剔除小图斑，确定出城市边界，如图2(g)所示。

实施例2：过程示意图见图3。

其处理步骤和实施例1相同，最终结果图如图3(g)所示。

实施例1和实施例2的检测结果均由实线标出在图中。

由此可见，本发明整体检测效果理想。此外，将现有技术(参考文献：XiaoP,WangX,FengX,etal.DetectingChina’sUrbanExpansionOverthePastThreeDecadesUsingNighttimeLightData[J].IEEEJournalofSelectedTopicsinAppliedEarthObservations&RemoteSensing,2014,7(10):4095-4106.)的检测结果(图2(h)、图3(h))和本发明结果(图2(g)、图3(g))作对比，通过统计正确率、错误率及kappa系数，比较分析本发明同现有技术的优势，统计结果见下表：

实施例1统计结果

	正确率(％)	错误率(％)	Kappa系数
				本发明方法	90.3	9.7	0.61
现有技术方法	82.6	17.4	0.45

实施例2统计结果

	正确率(％)	错误率(％)	Kappa系数
				本发明方法	93.4	6.6	0.70
现有技术方法	78.3	21.7	0.38

综合以上分析，本发明方法在正确率、错误率及Kappa系数上都要优于现有技术方法，说明本发明方法的可行性。

Claims (5)

1.一种融合多光谱遥感数据和夜间灯光遥感数据的城市边界提取方法，其特征在于，包括以下几个步骤：

步骤1：影像预处理；

对同一场景的多光谱遥感影像和夜间灯光遥感影像进行坐标转换，统一转为UTM投影坐标系；

将坐标转换后的夜间灯光遥感影像重采样至与坐标转换后的多光谱遥感影像相同的空间分辨率，获得Landsat图像数据I₀；

步骤2：自动提取正负训练样本；

首先利用Landsat图像数据I₀提供的光谱信息计算不透水面指数得到初始不透水面图像I₁，然后结合坐标转换后的夜间灯光遥感数据提供的初始灯光指数图像I₂，自动提取正负分类训练样本；

所述正负样本的自动提取过程如下：

A1：对初始不透水面图像I₁进行Otsu阈值分割，得到二值分割后的二值化不透水面图像I₃；

$I_3(i,j)=\left\{\begin{array}{cc}1& if\left(I_1\right(i,j)\&GreaterEqual;{\mathrm{th}}_1\\ 0& if\left(I_1\right(i,j)<{\mathrm{th}}_1\end{array}\right.$

其中，(i，j)表示图像中的像素坐标，th₁为Otsu阈值分割方法自适应选取的分割阈值；

A2：对初始灯光指数图像I₂进行Otsu阈值分割，得到二值分割后的二值化灯光指数图像I₄；

$I_4(i,j)=\left\{\begin{array}{cc}1& if\left(I_1\right(i,j)\&GreaterEqual;{\mathrm{th}}_2\\ 0& if\left(I_1\right(i,j)<{\mathrm{th}}_2\end{array}\right.$

其中，(i，j)表示图像中的像素坐标，th₂为Otsu阈值分割方法自适应选取的分割阈值；

A3：计算二值化灯光图像I₄中像素值为1的所有像素在初始灯光指数图像I₂中对应像素的像素值的平均值th₃，以及像素值为0的所有像素在初始灯光指数图像I₂中对应像素的像素值的平均值th₄；

A4：对二值化不透水面图像I₃进行格网划分，然后计算各个格网单位所包含像素对应的不透水面指数平均值a₁和灯光指数平均值a₂：

如果a₁≥th₁且a₂≥th₃，则将该网格中的所有像素提取至候选正样本集合P₁；

如果a₂≤th₄，则将该网格中的所有像素提取至候选负样本集合N₁；

A5：从候选正样本集合P₁和候选负样本集合N₁中分别随机提取50％的正负样本构造最终的正样本集合P₂和负样本集合N₂；

步骤3：利用步骤2获得的正负样本进行城市边界提取；

步骤3.1：首先提取用于描述正负样本的分类特征F＝P₂(f)+N₂(f),f＝{f₁,f₂,f₃,f₄}；

其中，f为影像特征集合，f₁为从Landsat图像数据I₀中提取的像素的7维光谱特征，f₂为从Landsat图像数据I₀中提取的像素的2维纹理特征，f₃为从初始不透水面图像I₁中提取的像素的1维的不透水指数特征，f₄为从初始灯光指数图像I₂中提取的像素的1维的灯光指数特征，P₂为正样本集合中所有像素对应影像中的位置，N₂为负样本集合中所有像素对应影像中的位置；

然后利用正负样本的分类特征训练SVM分类器，再用训练好的分类器对待分类的场景的影像特征集合f进行分类，得到城市区域和非城市区域；

步骤3.2：利用获得的城市区域，进行分类后处理，确定出城市边界；

所述分类后处理包括对城市区域填充内部空洞以及剔除小图斑操作，其中，剔除图斑的阈值为3×3像素。

2.根据权利要求1所述的方法，其特征在于，所述灯光遥感数据的重采样处理是指重采样方法采用三次卷积内插法。

3.根据权利要求1所述的方法，其特征在于，所述利用多光谱遥感数据提供的光谱信息计算不透水面指数得到初始不透水面图像I₁的步骤如下：

步骤2.1：遥感影像辐射校正；

利用ENVI软件，对Landsat图像数据I₀进行辐射校正，将影像上的DN值转换为反射率；

步骤2.2：不透水面指数计算；

采用归一化差值不透水面指数NDISI来提取步骤2中初始不透水面I₁，该指数计算公式为：

$NDISI=\frac{TIR-(MNDWI+NIR+{\mathrm{MIR}}_1)/3}{TIR+(MNDWI+NIR+{\mathrm{MIR}}_1)/3}$

其中，NIR、MIR₁和TIR分别为近红外、中红外1和热红外波段的反射率，MNDWI为归一化水体指数，计算公式为：

$MNDWI=\frac{Green-{\mathrm{MIR}}_1}{Green+{\mathrm{MIR}}_1}$

其中，Green为绿波段的反射率。

通过上述步骤计算得到多光谱遥感影像中对应的各个像素的不透水面指数，从而得到初始不透水面图像I₁。

4.根据权利要求1-3任一项所述的方法，其特征在于，所述对二值化不透水面图像I₃进行格网划分时，每个格网单元大小为8*8。

5.根据权利要求1-3任一项所述的方法，其特征在于，所述f₂为从Landsat图像数据I₀中提取的像素的2维纹理特征是指在3*3窗口内的均值和方差。