CN103824077A - 一种基于多源遥感数据的城市不透水层率息提取方法 - Google Patents

Abstract

一种基于多源遥感数据的城市不透水层信息提取方法，包括根据夜间灯光亮度影像，计算归一化灯光亮度并基于阈值分割方法对城市地表进行区域类型分割；利用灯光亮度、陆地表面温度与归一化土壤调节植被指数，构建增强型归一化不透水层指数；通过建立每一种类型分割区域的增强型归一化不透水层指数MNDISI与不透水层率之间的线性定量关系，进行不透水层的提取工作。本发明较好地结合了多源遥感数据的优势，有效抑制了光谱混淆对于不透水层率估算结果的不利影响，为城市土地利用和环境变化监测提供了一种新的途径。

Description

一种基于多源遥感数据的城市不透水层率息提取方法

技术领域

本发明属于遥感影像数据的信息提取领域，涉及一种基于多源遥感数据的城市不透水层信息提取方法。

背景技术

随着经济发展和城市化进程的加快，人类的各种生产建设活动正在日益改变着城市及其周边的自然环境和土地覆盖类型。不透水层被定义为由各种人工建筑，道路，机场等构成的具有不透水性的非自然地表面，不透水层率是指单位面积里不透水层所占的百分率大小。城市中大范围的不透水层会影响太阳辐射下垫面的热吸收与释放，引起“热岛”现象。在城市周边，不透水层的扩张会抑制植被对于大气中CO₂的吸收，对局部和全球范围内碳循环产生影响。因此，准确估算不透水层的分布和覆盖程度对于城市景观尺度上的生态环境研究和土地覆盖变化检测等都具有十分重要的意义。

目前，利用遥感数据提取不透水层的方法主要分为以下三种方法：

第一类是基于像素信息的分类方法，包括监督分类和非监督分类方法；

第二类是基于亚像元光谱信息分解的提取方法；

第三类是利用其他遥感信息例如植被指数NDVI与不透水层之间的相关性来进行提取的方法。

第一类基于像素的提取方法易受到地物光谱混淆问题的影响，近年来高分辨率遥感技术的发展使得局部范围内高精度的不透水层提取成为可能，但是由于波段数量较少，地物光谱混淆的问题依然存在，同时影像中高大建筑物的阴影也会对提取精度产生不利影响；第二类混合像元分解模型基于地表各组分间不存在交叉反射的假设，并不切合地表非朗伯性的实际情况，且像元分解精度很大程度上取决于端元选取的精度，因此该方法存在一定的不稳定性；第三类间接提取的方法主要针对的是区域及更大尺度上的不透水层提取，其结果往往不够精细。另外这种方法还容易受到其他外部因素的影响，例如NDVI就容易受到季节更替和土壤含水量的影响，对不透水层提取结果产生不利影响。

城市是整个自然界的特殊区域，不仅受到陆地，大气和植被的综合影响，更受到人类活动越来越频繁和持久的作用，与自然地貌相比，以不透水层为代表的人工地表具有光谱变化大，噪声强等特征，另一方面不透水层与其他自然地表覆盖类型如植被，水体和土壤的共同存在使得各种背景元素交错分布，对不透水层的提取造成了很大的干扰。

发明内容

针对以上问题，本发明提出一种基于多源遥感数据的城市不透水层提取方法并用于城市景观尺度上的不透水层分布及覆盖度估算。不同于以往的提取方法，本发明利用城市不透水层在陆地表面温度，植被指数与夜间灯光亮度数据中所反映的不同特性，以归一化遥感信息指数构建为框架，通过建立专题指数与不透水层之间的定量关系来进行不透水层的提取工作。

本发明所采用的技术方案是一种基于多源遥感数据的城市不透水层提取方法，包括以下步骤：

步骤a，根据夜间灯光亮度影像，获取各像素灯光亮度并记为OLS，计算归一化灯光亮度并基于阈值分割方法对城市地表进行区域类型分割；

步骤b，选取同一地区的多光谱影像和陆地表面温度影像分辨率，并以灯光亮度影像为基准进行配准、重采样和辐射矫正，得到真实反射率影像和各像素陆地表面温度LST；

步骤c，针对步骤b生成的真实反射率影像，建立改进的归一化水体指数MNDWI对水体进行剔除；

步骤d，针对步骤c剔除水体后的真实反射率影像，计算各像素的归一化土壤调节植被指数SAVI；

步骤e，针对步骤c剔除水体后的真实反射率影像，构建增强型归一化不透水层指数MNDISI，

$\mathrm{MNDISI}=\frac{({\mathrm{OLS}}_{\mathrm{NOR}}+{\mathrm{LST}}_{\mathrm{NOR}})-({\mathrm{SAVI}}_{\mathrm{NOR}}+\mathrm{MIR})}{({\mathrm{OLS}}_{\mathrm{NOR}}+{\mathrm{LST}}_{\mathrm{NOR}})+({\mathrm{SAVI}}_{\mathrm{NOR}}+\mathrm{MIR})}$

式中，OLS_NOR、LST_NOR和SAVI_NOR分别为步骤a所得灯光亮度OLS、步骤b所得陆地表面温度LST与步骤d所得归一化土壤调节植被指数SAVI分别进一步归一化到[0,1]范围的结果，MIR为中红外波段的反射率；

步骤f，根据同一地区的高分辨率影像获取训练样本，包括首先利用面向对象分类将每一个高分辨率影像像素标记为不透水层或者透水层，然后下采样到夜间灯光亮度影像的分辨率下并计算得到每一个像素的参考不透水层百分率数据，最后采用随机取样方法在步骤a确定的每一分割区域内选取部分像素的参考不透水层百分率数据作为不透水层率的训练样本；

步骤g，根据步骤e所得各像素的增强型归一化不透水层指数MNDISI，利用步骤f所得训练样本，拟合得到每一种类型分割区域的增强型归一化不透水层指数MNDISI与不透水层率之间的线性定量关系；

步骤h，参照步骤a所得区域类型分割结果和步骤g所得应类型的线性定量关系，得到城市不透水层率估算结果。

而且，步骤a中，归一化灯光亮度计算公式如下，

$I_j=\underset{1}{\overset{r}{\mathrm{\Σ}}}{\mathrm{DN}}_i\·\frac{n_i}{r\·N}$

式中，I_j为某一区域j的归一化灯光亮度，r为灯光亮度影像的最大亮度值，DN_i为区域j内第i等级的灰度值，n_i为区域j内第i灰度等级的像元总数，N为区域j内所有灯光像元总数。

而且，步骤a中，以归一化灯光亮度为灯光强度指标，阈值划分方法的规则如下，

如果某区域灯光强度指标大于等于预设阈值a1时，则定义该区域为中心城区；

如果某区域灯光强度指标小于预设阈值a1且大于等于预设阈值a2时，定义该区域为城乡结合区；

如果某区域灯光强度指标小于预设阈值a2，则定义该区域为附属村落。

而且，步骤c中，改进的归一化水体指数MNDWI计算公式如下，

$\mathrm{MNDWI}=\frac{G-\mathrm{MIR}}{G+\mathrm{MIR}}$

式中，G、MIR分别为绿光波段和中红外波段的反射率，当某像素计算所得归一化水体指数MDNWI大于预设阈值a3被认为是水体。

而且，步骤d中，归一化土壤调节植被指数SAVI计算公式如下，

$\mathrm{SAVI}=\frac{(\mathrm{NIR}-R)(1+l)}{(\mathrm{NIR}+R+l)}$

式中，NIR、R分别为近红外波段和红光波段的反射率，l为土壤适应性因子。

而且，步骤e中，OLS_NOR、LST_NOR和SAVI_NOR计算公式如下，

${\mathrm{OLS}}_{\mathrm{NOR}}=\frac{\mathrm{OLS}-{\mathrm{OLS}}_{\min }}{{\mathrm{OLS}}_{\max }-{\mathrm{OLS}}_{\min }},{\mathrm{LST}}_{\mathrm{NOR}}=\frac{\mathrm{LST}-{\mathrm{LST}}_{\min }}{{\mathrm{LTS}}_{\max }-{\mathrm{LST}}_{\min }},{\mathrm{SAVI}}_{\mathrm{NOR}}=\frac{\mathrm{SAVI}-{\mathrm{SAVI}}_{\min }}{{\mathrm{SAVI}}_{\max }-{\mathrm{SAVI}}_{\min }}$

其中，OLS，LST和SAVI分别为步骤a所得灯光亮度，步骤b所得陆地表面温度与步骤d所得归一化土壤调节植被指数，OLS_max、LST_max和SAVI_max分别为灯光亮度、陆地表面温度与SAVI的最大值，OLS_min、LST_min和SAVI_min分别为灯光亮度、陆地表面温度与归一化土壤调节植被指数的最小值。

而且，步骤g中，第k种类型分割区域的增强型归一化不透水层指数MNDISI与不透水层率之间的线性定量关系如下式，

k=1表示分割区域类型为中心城区，k=2表示分割区域类型为城乡结合区，k=3表示分割区域类型为乡村地区，c_k表示位于第k种类型分割区域的第c个像素点，ISP_ck表示像素c_k不透水层率拟合值，MNDISI_ck代表像素c_k的增强型归一化不透水层指数，b_k0和b_k1为第k种类型分割区域的线性拟合系数。

本发明提供的技术方案的有益效果为：提出了利用多源遥感数据指数方法提取城市不透水层的思想，较好地结合了多源遥感数据的优势，有效抑制了光谱混淆对于不透水层率估算结果的不利影响，为城市土地利用和环境变化监测提供了一种新的途径。

附图说明

图1为本发明实施例的流程图。

具体实施方式

具体实施时，本发明技术方案所提供流程可由本领域技术人员采用计算机软件技术实现自动运行。为了更好地理解本发明的技术方案，下面结合附图和实施例对本发明做进一步的详细说明。

参见图1，本发明实施例包括以下步骤：

步骤a，首先对获取的夜间灯光亮度影像进行裁剪并重投影到WGS84坐标系下。根据夜间灯光亮度影像，获取各像素灯光亮度并记为OLS；然后以行政区（例如区/县）为区域单元计算归一化灯光亮度并基于阈值分割方法对城市地表进行景观区域分割。归一化灯光亮度计算公式见公式（1）：

$I_j=\underset{1}{\overset{r}{\mathrm{\Σ}}}{\mathrm{DN}}_i\·\frac{n_i}{r\·N}---\left(1\right)$

式中I_j为某一区域j的归一化灯光亮度，r为灯光亮度影像的最大亮度值，DN_i为区域j内第i等级的灰度值，n_i为区域j内第i灰度等级的像元总数，N为区域j内所有灯光像元总数。以计算结果I_j为区域j的灯光强度指标，分别选取阈值对研究区内的区/县进行城市景观类型划分。实施例具体的阈值划分规则为：

（1）如果某区域灯光强度指标大于等于a1时，则定义该区域为中心城区；

（2）如果某区域灯光强度指标小于a1，大于等于a2时，定义该区域为城乡结合区；

（3）如果某区域灯光强度指标小于a2，则定义该区域为附属村落。

根据该规则，实施例进行阈值划分得到3种类型的分割区域，中心城区、城乡结合区和附属村落。

其中，影像数据裁剪重投影属于影像处理领域的常规方法，具体处理流程不再赘述。具体实施时，阈值a1、a2可由本领域技术人员自行预先设置为经验值，实施例中a1=0.8，a2=0.35。

步骤b，对获取的同一地区的多光谱影像和陆地表面温度影像进行裁剪并重投影到WGS84坐标系下。以灯光亮度影像的空间分辨率为基准对多光谱影像和陆地表面温度影像进行配准和重采样。然后将配准后多光谱影像和陆地表面温度影像的波段DN值进行辐射矫正转换为真实反射率和陆地表面温度，得到真实反射率影像，得到的各像素陆地表面温度记为LST。

影像数据裁剪重投影、配准、重采样、辐射矫正等属于影像处理领域的常规方法，具体处理流程不再赘述。

步骤c，针对步骤b得到的真实反射率影像，采用改进的归一化水体指数ModifiedNormalized Digital Water Index（MNDWI）剔除水体信息。实施例采用的改进的MNDWI数学表达式见公式（2）：

$\mathrm{MNDWI}=\frac{G-\mathrm{MIR}}{G+\mathrm{MIR}}---\left(2\right)$

式中G，MIR分别为绿光波段和中红外波段的反射率，当计算所得归一化水体指数MDNWI﹥a3被认为是水体像元，可通过掩膜处理剔除。具体实施时，阈值a3可由本领域技术人员自行预先设置为经验值，实施例中a3=0.1。MDNWI≤0.1的像素被保留为有效像素，具体实施时，步骤a、b、c可并行处理，也可先进行步骤c提取有效像素，步骤a中则只需对有效像素的灯光亮度进行处理。

步骤d，针对步骤c剔除水体后的真实反射率影像，根据城市地区的低植被覆盖度特点，选取归一化土壤调节植被指数Soil Adjusted Vegetation Index（SAVI）作为植被丰度信息指数。实施例中的各像素的SAVI数学表达式见公式（3）：

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

式中NIR，R分别为近红外波段和红光波段的反射率，l为土壤适应性因子，数值范围从0到1。具体实施时，参数l可由本领域技术人员自行预先设置为经验值，实施例中确定l=0.75。

步骤e，对步骤a，b，c，d得到的结果，针对步骤c剔除水体后的真实反射率影像，构建增强型归一化不透水层指数Modified Normalized Digital Impervious Surface Index（MNDISI）。实施例中MNDISI构建表达式见公式（4）：

$\mathrm{MNDISI}=\frac{({\mathrm{OLS}}_{\mathrm{NOR}}+{\mathrm{LST}}_{\mathrm{NOR}})-({\mathrm{SAVO}}_{\mathrm{NOR}}+\mathrm{MIR})}{({\mathrm{OLS}}_{\mathrm{NOR}}+{\mathrm{LST}}_{\mathrm{NOR}})+({\mathrm{SAVI}}_{\mathrm{NOR}}+\mathrm{MIR})}---\left(4\right)$

式中OLS_NOR、LST_NOR和SAVI_NOR分别为步骤a得到的灯光亮度，步骤b得到的陆地表面温度与步骤d得到的归一化土壤调节植被指数分别进一步归一化到[0,1]范围的结果，实施例中OLS_NOR、LST_NOR和SAVI_NOR的数学表达式分别见公式（5）、（6）、（7）：

${\mathrm{OLS}}_{\mathrm{NOR}}=\frac{\mathrm{OLS}-{\mathrm{OLS}}_{\min }}{{\mathrm{OLS}}_{\max }-{\mathrm{OLS}}_{\min }}---\left(5\right)$

${\mathrm{LST}}_{\mathrm{NOR}}=\frac{\mathrm{LST}-{\mathrm{LST}}_{\min }}{{\mathrm{LST}}_{\max }-{\mathrm{LST}}_{\min }}---\left(6\right)$

${\mathrm{SAVI}}_{\mathrm{NOR}}=\frac{\mathrm{SAVI}-{\mathrm{SAVI}}_{\min }}{{\mathrm{SAVI}}_{\max }-{\mathrm{SAVI}}_{\min }}---\left(7\right)$

其中，OLS，LST和SAVI分别为步骤a得到的灯光亮度，步骤b得到的陆地表面温度与步骤d得到的归一化土壤调节植被指数。OLS_max、LST_max和SAVI_max分别为灯光亮度，陆地表面温度与SAVI的最大值。OLS_min、LST_min和SAVI_min分别为灯光亮度，陆地表面温度与归一化土壤调节植被指数的最小值。

步骤f，获取同一地区的高分辨率影像并利用面向对象分类方法和像素下采样获取不透水层率计算训练样本数据。具体操作为：首先利用面向对象分类将每一个高分辨率影像像素标记为不透水层或者透水层。然后下采样到夜间灯光亮度影像的分辨率下得到每一个像素的参考不透水层百分率数据。最后采用随机取样方法在步骤a确定的每一种分割区域内从选取部分像素点（实施例取5%的像素点）的参考不透水层百分率数据作为不透水层率计算的训练样本。

影像面向对象分类和像素下采样均为现有技术，在此不再赘述具体流程。

步骤g，对步骤c得到的有效像素，利用步骤f生成的结果，以步骤a的区域分割结果分别建立MNDISI与不透水层率之间的线性定量关系。具体操作为：首先对应区域分割结果将步骤f得到的不透水层率训练样本数据分配到各种分割区域。然后根据步骤e所得各像素的MNDISI，利用最小二乘法对每一种类型分割区域的MNDISI和不透水层率训练样本数据进行线性拟合。实施例中线性拟合数学表达式见公式（8）：

k表示不同类型的分割区域，其取值可以为1（中心城区），2（城乡结合区）或3（乡村地区），c_k表示位于第k种类型分割区域的第c个像素点，ISP_ck表示该像素c_k不透水层率拟合值，MNDISI_ck代表该像素c_k的增强型归一化不透水层指数，b_k0和b_k1为第k种类型分割区域的线性拟合系数。每种类型的分割区域可得一个相应的线性定量关系。

最小二乘拟合法为现有技术，其具体流程不再赘述。

步骤h，将步骤g生成的线性定量关系对应分割区域应用于城市不透水层率估算。参照步骤a的区域类型分割结果和步骤g生成的相应类型的线性拟合公式，利用步骤e生成的MNDISI即可计算得到每个类型的不透水层率。

综上所述，本发明提出的基于多源遥感数据的城市不透水层信息提取方法，在以往注重利用不透水层的光谱信息进行提取的基础上，加入了夜间灯光亮度影像这一新数据源，能够有效规避不透水层与其他地物例如土壤的光谱混淆问题，增加了提取结果的准确度与稳定性，由于数据源便于获取，操作步骤简单，该方法更加符合实际的应用需求。

为说明本发明的实用性和效果，提供以下实验：

利用2010年美国国防气象卫星计划(DMSP)下辖的线性扫描系统（OLS）夜间稳定灯光亮度影像（DMSP/OLS）和2010年8月中分辨率成像光谱仪（MODIS）获取的地表反射率和陆地表面温度影像影像对北京地区的不透水层率进行提取，并利用同一区域的快鸟高分辨率影像获取参考不透水层率数据分别用于MNDISI—ISP线性拟合和实施例所提供方法的有效性评价，可利用步骤f选取不透水层率计算的训练样本后剩余的95%的像素点的参考不透水层百分率数。

从主观视觉效果对比和客观精度评价指标两方面对实施例所提供方法的有效性进行评价，例如对中心城区、密云、顺义三个样地的Landsat TM假彩色图像，分别得到利用亚像元光谱分解得到的不透水层提取结果、利用MNDISI得到的不透水层提取结果。从目视效果可以看出：由于主城区的地表覆盖类型主要为不透水层和植被，裸土的成分较少，因此利用两种方法估算的方法结果呈现较好的一致性。顺义地区的地表覆盖较为复杂，既有超大型的不透水层设施如机场，大型居民地等，也有大片的植被，农田和裸土分布，光谱分解方法受限于端元的选择精度，结果出现了一定的混淆性，总体的ISP水平比实际值偏高，利用MNDISI估算的效果相对较佳，大型的不透水层设施ISP值最高，其他零散分布的聚落其次，而裸土、植被和农作物覆盖的区域则ISP值很低（ISP<5%），在密云地区利用MNDISI估算ISP的优势则更加明显，由于密云水库周围既有居民地分布，同时也有大片未开发的裸土存在，二者的光谱特性较为相似，加上影像本身的大气影响，导致可供选择的端元很少，光谱分解得到的ISP分布准确率不高。利用MNDISI估算ISP则可以得到较为准确的估算结果。在密云地区，主要的不透水层分布于该地区的西南部，另外还有零散分布的一些村落区域ISP较高，其他地区的ISP水平都较低，这与实际情况是一致的。

客观评价指标选取相关系数R²，均方根误差RMSE和平均相对偏差MAE三种指标，其计算公式见公式（9）—（11）。其中

为第i个样本的ISP估算值，ISP_i为第i个样本的ISP参考真实值，ISP_mean为所有验证样本的ISP参考真实值的均值，N为样本数目。

（1）均方根误差RMSE：是一种刻画样本计算值与真实值偏差程度的一种评价指标，均方根误差RMSE越小说明方法稳健性越好。其数学表达式见公式（9）：

$\mathrm{RMSE}=\sqrt{\frac{\underset{1}{\overset{N}{\mathrm{\&Sigma;}}}{({\stackrel{\&OverBar;}{\mathrm{ISP}}}_i-{\mathrm{ISP}}_{\mathrm{mean}})}^2}{N}}---\left(9\right)$

（2）相关系数R²：又称为确定系数，用来衡量实际测量值与提取值之间的相关性强弱，其数学表达式见公式（10）：

$R^2=\frac{\underset{1}{\overset{N}{\mathrm{\&Sigma;}}}({\stackrel{\&OverBar;}{\mathrm{ISP}}}_i-{\mathrm{ISP}}_{\mathrm{mean}})}{\underset{1}{\overset{N}{\mathrm{\&Sigma;}}}{(\stackrel{\&OverBar;}{\mathrm{ISP}}-{\mathrm{ISP}}_{\mathrm{mean}})}^2}---\left(10\right)$

（3）平均相对偏差MAE：其数学表达式见公式（11）：

$\mathrm{MAE}=\frac{1}{N}\&CenterDot;\underset{1}{\overset{N}{\mathrm{\&Sigma;}}}\left|\frac{{\mathrm{ISP}}_i-{\stackrel{\&OverBar;}{\mathrm{ISP}}}_i}{{\mathrm{ISP}}_i}\right|---\left(11\right)$

按照以上提取有效性评价指标操作要求得到评价指标系数如表1所示：

表1

从表1的实验结果可知，本发明提出的基于多源数据的城市不透水层提取方法在各指标上均比传统的光谱分解方法有一定的优势，特别在城市化水平较低的城乡结合部和乡村地区，光谱混淆现象较为显著，使得光谱分解方法精度较低，而MNDISI方法则提取精度较高。

Claims (7)

1.一种基于多源遥感数据的城市不透水层提取方法，包括以下步骤：

步骤a，根据夜间灯光亮度影像，获取各像素灯光亮度并记为OLS，计算归一化灯光亮度并基于阈值分割方法对城市地表进行区域类型分割；

步骤b，选取同一地区的多光谱影像和陆地表面温度影像分辨率，并以灯光亮度影像为基准进行配准、重采样和辐射矫正，得到真实反射率影像和各像素陆地表面温度LST；

步骤c，针对步骤b生成的真实反射率影像，建立改进的归一化水体指数MNDWI对水体进行剔除；

步骤d，针对步骤c剔除水体后的真实反射率影像，计算各像素的归一化土壤调节植被指数SAVI；

步骤e，针对步骤c剔除水体后的真实反射率影像，构建增强型归一化不透水层指数MNDISI，

$\mathrm{MNDISI}=\frac{({\mathrm{OLS}}_{\mathrm{NOR}}+{\mathrm{LST}}_{\mathrm{NOR}})-({\mathrm{SAVI}}_{\mathrm{NOR}}+\mathrm{MIR})}{({\mathrm{OLS}}_{\mathrm{NOR}}+{\mathrm{LST}}_{\mathrm{NOR}})+({\mathrm{SAVI}}_{\mathrm{NOR}}+\mathrm{MIR})}$

式中，OLS_NOR、LST_NOR和SAVI_NOR分别为步骤a所得灯光亮度OLS、步骤b所得陆地表面温度LST与步骤d所得归一化土壤调节植被指数SAVI分别进一步归一化到[0,1]范围的结果，MIR为中红外波段的反射率；

步骤f，根据同一地区的高分辨率影像获取训练样本，包括首先利用面向对象分类将每一个高分辨率影像像素标记为不透水层或者透水层，然后下采样到夜间灯光亮度影像的分辨率下并计算得到每一个像素的参考不透水层百分率数据，最后采用随机取样方法在步骤a确定的每一分割区域内选取部分像素的参考不透水层百分率数据作为不透水层率的训练样本；

步骤g，根据步骤e所得各像素的增强型归一化不透水层指数MNDISI，利用步骤f所得训练样本，拟合得到每一种类型分割区域的增强型归一化不透水层指数MNDISI与不透水层率之间的线性定量关系；

步骤h，参照步骤a所得区域类型分割结果和步骤g所得应类型的线性定量关系，得到城市不透水层率估算结果。

2.根据权利要求1所述基于多源遥感数据的城市不透水层提取方法，其特征在于：步骤a中，归一化灯光亮度计算公式如下，

$I_j=\underset{1}{\overset{r}{\mathrm{\&Sigma;}}}{\mathrm{DN}}_i\&CenterDot;\frac{n_i}{r\&CenterDot;N}$

式中，I_j为某一区域j的归一化灯光亮度，r为灯光亮度影像的最大亮度值，DN_i为区域j内第i等级的灰度值，n_i为区域j内第i灰度等级的像元总数，N为区域j内所有灯光像元总数。

3.根据权利要求1所述基于多源遥感数据的城市不透水层提取方法，其特征在于：步骤a中，以归一化灯光亮度为灯光强度指标，阈值划分方法的规则如下，

如果某区域灯光强度指标大于等于预设阈值a1时，则定义该区域为中心城区；

如果某区域灯光强度指标小于预设阈值a1且大于等于预设阈值a2时，定义该区域为城乡结合区；

如果某区域灯光强度指标小于预设阈值a2，则定义该区域为附属村落。

4.根据权利要求1所述基于多源遥感数据的城市不透水层提取方法，其特征在于：步骤c中，改进的归一化水体指数MNDWI计算公式如下，

$\mathrm{MNDWI}=\frac{G-\mathrm{MIR}}{G+\mathrm{MIR}}$

式中，G、MIR分别为绿光波段和中红外波段的反射率，当某像素计算所得归一化水体指数MDNWI大于预设阈值a3被认为是水体。

5.根据权利要求1所述基于多源遥感数据的城市不透水层提取方法，其特征在于：步骤d中，归一化土壤调节植被指数SAVI计算公式如下，

$\mathrm{SAVI}=\frac{(\mathrm{NIR}-R)(1+l)}{(\mathrm{NIR}+R+l)}$

式中，NIR、R分别为近红外波段和红光波段的反射率，l为土壤适应性因子。

6.根据权利要求1所述基于多源遥感数据的城市不透水层提取方法，其特征在于：步骤e中，OLS_NOR、LST_NOR和SAVI_NOR计算公式如下，

${\mathrm{OLS}}_{\mathrm{NOR}}=\frac{\mathrm{OLS}-{\mathrm{OLS}}_{\min }}{{\mathrm{OLS}}_{\max }-{\mathrm{OLS}}_{\min }},{\mathrm{LST}}_{\mathrm{NOR}}=\frac{\mathrm{LST}-{\mathrm{LST}}_{\min }}{{\mathrm{LTS}}_{\max }-{\mathrm{LST}}_{\min }},{\mathrm{SAVI}}_{\mathrm{NOR}}=\frac{\mathrm{SAVI}-{\mathrm{SAVI}}_{\min }}{{\mathrm{SAVI}}_{\max }-{\mathrm{SAVI}}_{\min }}$

其中，OLS，LST和SAVI分别为步骤a所得灯光亮度，步骤b所得陆地表面温度与步骤d所得归一化土壤调节植被指数，OLS_max、LST_max和SAVI_max分别为灯光亮度、陆地表面温度与SAVI的最大值，OLS_min、LST_min和SAVI_min分别为灯光亮度、陆地表面温度与归一化土壤调节植被指数的最小值。

7.根据权利要求1或2或3或4或5或6所述基于多源遥感数据的城市不透水层提取方法，其特征在于：步骤g中，第k种类型分割区域的增强型归一化不透水层指数MNDISI与不透水层率之间的线性定量关系如下式，

k=1表示分割区域类型为中心城区，k=2表示分割区域类型为城乡结合区，k=3表示分割区域类型为乡村地区，c_k表示位于第k种类型分割区域的第c个像素点，ISP_ck表示像素c_k不透水层率拟合值，MNDISI_ck代表像素c_k的增强型归一化不透水层指数，b_k0和b_k1为第k种类型分割区域的线性拟合系数。

Images