不透水面（impervious surface）是指由水泥、沥青等不透水面建筑材料构成的阻止地表水下渗到土壤的人工地表面，比如建筑屋顶、道路等，与具有透水性的植被相对^[1]。已有研究表明，不透水面的增加会导致城市内涝^[2]、城市热岛效应加剧^[3]等一系列环境问题。因此，不透水面分布已经成为评价城市生态环境质量的重要考核指标之一，对城市环境规划与治理及城市发展规划具有重要意义^[4]。

随着城市化进程的加速以及遥感对地观测技术的发展，国内外研究者对不透水面提取做了诸多研究，主要方法包括：线性光谱混合分析法^[5]、指数法^[6]、决策树模型法^[7]、回归模型法^[8]、基于分类的方法^[9]等。Wu等^[5]基于Landsat ETM+遥感影像，将不透水面分为高、低反照率不透水面，选择了4种端元进行混合像元分解。查勇等^[6]发现在近红外波段不透水面的反射率低于植被的反射率，而在短波红外波段不透水面的反射率高于植被的反射率，因此利用近红外波段和短波红外波段设计了归一化建筑指数（normalized difference buklt-up index，NDBI）。Deng等^[10]利用缨帽变换的3个分量分别代表不透水面、植被和水体，设计了BCI指数。岳玉娟等^[11]根据不透水面与植被的负相关关系，通过计算植被覆盖度来提取不透水面，这种方法提取的不透水面的均方根误差为0.402，精度较低。

本文基于石家庄中心城区的Landsat8遥感影像，对比分析线性光谱混合分析法（linear spectral mixture analysis，LMSA）和归一化差值不透水面指数（normalized difference impervious surface index，NDISI）在石家庄城市不透水面信息提取中的精确度，并分析石家庄中心城区不透水面的空间分布格局，对于石家庄城市规划与环境保护具有重要的参考价值和现实意义。

1 研究区与数据来源 1.1 研究区概况

本文以河北省石家庄市中心城区作为研究区，位于北纬37°27′~38°47′，东经113°30′~115°20′之间。研究区海拔高度在30-100 m，以平原为主，土地利用类型主要包括水体、不透水地表、裸地及耕地等。

早期的石家庄市由内向外成同心圆式连续扩展的发展模式，但由于这种发展模式引发了一系列问题，所以渐渐引入了城市分散组团式的发展模式，即由单一核心城市向多核心、组团方向发展。基于该发展模式的思想，在2001年，撤销了郊区，增设裕华区，并将正定县2个镇归为长安区，扩大了石家庄城市规模。2014年9年，桥东区被撤销，同时增设了3个市辖区，石家庄城市规模进一步扩张^[12]。目前石家庄市中心城区有裕华区、桥西区、长安区、新华区，总面积约为405.8 km²。研究区地理位置如图 1所示。

1.2 数据源

本文获取了2019年5月20日的Landsat8中分辨率多光谱影像数据进行城市不透水面提取技术对比研究，以同期的谷歌影像数据作为精度评估时目视解译判读的参考数据。

其中Landsat8影像数据从USGS地球资源观测和科学中心科学处理部门（http://espa.cr.usgs.gov/）获取，轨道号为P124R33、P124R34，所用影像云量低于10%，影像质量较高。

谷歌影像从91卫图获取，分辨率为0.12 m，因云、阴影等影响，2019年5月20日的数据不能满足要求，所以选择了2018年5月和2019年5月的影像数据，时间上与实验数据接近且分辨率高，可以作为精度评估的参考数据。

1.3 数据预处理

本研究使用的是经过了几何校正和辐射校正的地表反射率产品，所以还需要做的预处理工作有影像拼接和裁剪以及水体掩膜。

1.3.1 影像拼接与裁剪

由于一幅影像不能覆盖整个研究区，本研究拼接了跨越整个研究区的两幅影像，因影像是同一时段的，所以影像拼接效果较好，无色差；再用石家庄市中心城区的矢量数据裁剪拼接后的影像得到石家庄市中心城区的遥感影像。

1.3.2 水体掩膜

由于水体和阴影与低反照率地物有相似的光谱特征，从而干扰不透水面的提取，因此要对研究区的水体和阴影做掩膜处理。由于研究区为平原区且所用数据为中分辨率影像，所以可以忽略阴影的影响。因此，只需通过修正归一化水体指数^[13]（modified normalized difference water index，MNDWI）对水体进行去除，其公式如下：

$ \text{ MNDWI }=(\text{ Green }-\text{MIR})/(\text{ Green }+\text{ MIR }) $

(1)

式中，MIR为中红外波段；Green为绿光波段；分别为Landsat8 0LI影像的第7波段和第3波段。

使用波段运算工具计算水体指数，再通过反复试验结合目视判读，当结果与真实情况最为吻合时确定阈值提取出水体，最后通过掩膜工具将水体去除。

2 研究方法 2.1 线性光谱混合分析法

空间分辨率为30 m的Landsat8 OLI影像中存在大量混合像元。混合像元是指那些包含多种地物类型的像元，所以为获取某类地物的分布图，必须确定组成混合像元的各种地物及其比例。以往研究者多用光谱混合分析法解决此类问题以获取某种地物的纯净光谱，即端元光谱。光谱混合分析法中比较成熟的方法是LSMM。该模型定义为像元在某一波段的反射率是该像元内各个端元的反射率的线性组合，而系数则是其基本组分在像元中的面积比例^[14]。其模型如下：

$ {{R}_{i}}=\sum\limits_{k=1}^{\text{N}}{{{f}_{k}}}{{R}_{ki}}+{{e}_{i}} $

(2)

式中，R_i是波段i的反射率；f_k是端元k所占的比例；R_ki是端元k在波段i的反射率；N是像元内端元的总数；e_i是残差。其中，f_k需满足两个约束条件：①端元k所占比例非负，即f_k ≥ 0；②各端元所占比例之和为1，即$\sum\limits_{k=1}^{N}{{{f}_{k}}}=1$。通过提取端元波谱，利用带约束的线性光谱分解模型得到各个端元在每个像元中比例可以避免传统混合像元分解过程中丰度图出现负值或大于1情况的发生^[15]。

线性光谱混合分解主要需要经过以下几个步骤：确定端元的数目、提取端元、线性光谱混合分解、精度评估^[16]。

2.1.1 确定端元数目

根据线性光谱混合模型的定义可知，确定端元的个数及各端元在各波段的反射率是线性光谱混合分解的关键。通过对石家庄中心城区地物的分析，并结合Ridd^[17]提出的植被-不透水面-土壤模型（V-IS模型），本研究选择高反照率不透水面、植被、土壤、低反照率不透水面四种端元。

2.1.2 基于纯净像元指数提取端元

由于从散点图上直接选取端元有一定困难，而且受人为因素影响大，无法保证结果的准确性，所以本研究使用基于纯净像元指数的方法来提取端元。

原始遥感影像中存在噪声，同时各个波段间有不同程度的相关性，导致数据冗余，所以计算纯净像元指数前需要首先对遥感影像进行最小噪声分离变换（minimum noise fraction，MNF），去除图像中的噪声以及不同波段间的相关性，以保证结果的精度和信息的丰富程度。

为了提高端元选择的准确性，需要在高纯度的像元内选择端元。所以本研究利用最小噪声分离变换后的前4个波段计算纯净像元指数，得到“纯净像元指数图像”。该图像中每个像元的DN值与像元被标记为极值的次数相对应，DN值越高其属于纯净像元的概率越大，所以通过反复实验设置10为阈值提取出纯净像元，缩小了端元选择的范围，从而将端元选择的范围集中于高纯度的像元内。

为进一步提高精度，减少肉眼识别造成的误差，本研究将选择的纯净像元投影到ENVI软件中的N维可视化工具中，以实现像元的可视化操作。利用N维可视化工具，可同时使用2个以上的特征波段进行交互分析，且相同特征的纯净像元自动聚集在一起，因此通过N维可视化工具选择端元可以进一步提高精度。

本研究将纯净像元指数大于10的像元输入到N维可视化空间中，选择最小噪声分离变换的前三个波段构成三维散点图选择符合条件的像元所形成的端元集合，绘制其平均波谱曲线。

2.1.3 线性光谱混合分解

根据基于纯净像元指数提取的端元，使用线性光谱混合分解，对石家庄市中心城区的Landsat 8 OLI影像进行混合像元分解，得到石家庄中心城区植被、土壤、高反照率不透水面和低反照率不透水面4种端元的丰度图和RMSE（root mean square error）分布图像。

线性光谱混合分析法的关键即端元数目的确定与端元光谱的提取。本研究依据Ridd提出的V-I-S模型选择高反照率不透水面、低反照率不透水面、植被和土壤四个端元能够有效代表石家庄中心城区各地物类型；同时利用纯净像元指数结合n维可视化工具提取端元光谱能够大大降低肉眼识别的误差，从而得到更为准确的端元。

2.2 归一化差值不透水面指数法

根据不透水面区别于其他地类的光谱特性的差异构建遥感指数，可以快速高效地提取大范围不透水面信息。由于不透水面包括建筑屋顶、道路等多种地物，其组成复杂，不容易找出合适的波段范围来得到不透水面，因此徐涵秋等^[18]研究了不透水面的共性，发现不透水面在热红外波段的反射率高于其他地类，在近红外波段、可见光波段和中红外波段的反射率均低于其他地类的反射率，所以将近红外波段、可见光波段和中红外波段作为不透水面的弱反射组，将热红外波段作为其高反射组，创建了归一化差值不透水面指数（NDISI）来突出不透水面信息。这个指数通过多波段的组合来提取不透水面，能够降低提取不透水面的指数中存在的错分误差。NDISI的指数表达为：

$ \text{NDISI}=\frac{\operatorname{TIR}-\frac{(\text{VIS}1+\text{NIR}+\text{MIR}1)}{3}}{\operatorname{TIR}+\frac{(\text{VIS}1+\text{NIR}+\text{MIR}1)}{3}} $

(3)

式中，NIR是近红外波段；MIR1是中红外中的某一波段；TIR是热红外波段；VIS1是可见光中某一波段；水体指数MNDWI可以代替可见光波段。由于本实验预先掩膜水体，所以使用可见光中某一波段进行计算。通过多次实验，本研究选择Landsat8中的波段11、波段3、波段5、波段7计算指数。

由于热红外波段影像分辨率为100 m，所以在计算指数前需要在ENVI软件中先对其进行重采样，使其分辨率与多光谱影像分辨率一致。

使用ENVI软件中的Band Math工具计算NDISI得到指数影像后，再利用通过阈值选择感兴趣区功能结合目视判读确定阈值，通过Band Math工具对指数影像进行二值化，得到不透水面分布图。

2.3 精度评估

在ENVI软件中随机生成80个验证区域，窗口大小为90 m×90 m（3×3像元），以降低不同影像几何误差的影响。分别统计每个验证区域两种方法提取的不透水面的比例。另外利用生成的90 m×90 m的随机样本下载同期高分辨（0.12 m）的谷歌影像，在ArcGIS软件中通过目视解译结合数字化计算其不透水面比例作为地面真实参考数据。采用均方根误差RMSE作为精度评价的指标，计算公式如下：

$ \text{RMSE}=\sqrt{\frac{\sum\limits_{k=1}^{\text{N}}{{{\left( {{X}_{k}}-{{x}_{k}} \right)}^{2}}}}{~\text{N}}} $

(4)

式中，X_k是Landsat影像中不透水面丰度的估计值；x_k是高分辨率谷歌影像中不透水面丰度的“真实值”；N为样本数目。

3 结果与分析 3.1 基于LSMA方法的不透水面提取 3.1.1 端元提取结果

通过PPI结合N维可视化工具获取的端元光谱如图 2所示。由图可知，高反照率不透水面的反射率高于其他地类，可以有效避免传统线性光谱混合分析法中高反照率不透水面与植被的混淆；在近红外波段，低反照率不透水面的反射率低于其他地类；在短波红外波段，低反照率不透水面的反射率高于植被。因此，该实验选择的高、低反照率不透水面端元能够较容易与其他地类区分。

3.1.2 基于LSMA方法的不透水面提取结果

线性光谱分解结果的均方根误差为0.036，精度满足要求。将线性光谱混合分析法得到的高反照率不透水面比例图和低反照率不透水面比例图相加得到石家庄中心城区不透水面比例图如图 3所示。在ArcGIS软件中对石家庄中心城区高分辨率Google影像进行数字化，并统计得到基于LSMA方法的提取结果的RMSE为0.246。

由图 3可以看出，长安区和新华区北部的周边范围有部分不透水面丰度偏低的区域；裕华区东南角及内部也有部分不透水面丰度偏低的区域，而其余地区不透水面丰度比较高。

3.2 基于NDISI的不透水面提取结果

将用NDISI指数计算得到的图像进行二值化处理，通过反复试验，选择0.32作为阈值，得到基于NDISI的石家庄不透水面分布图如图 4所示。

由图 4可以看出，该指数提取的不透水面分布信息与线性光谱混合分析法提取的不透水面分布信息近似，均是长安区和新华区北部的周边范围为透水面，裕华区东南角及内部部分区域为透水面，而其余地区几乎全是不透水面。

NDISI指数的提取结果用ArcGIS软件分析可得：石家庄市中心城区不透水面总面积为296.08 km²，占石家庄市中心城区总面积的72.96%。其中，长安区不透水面面积为97.70 km²，占长安区总面积的69.49%；裕华区不透水面面积为78.78 km²，占裕华区总面积的73.48%；新华区不透水面面积为65.98 km²，占新华区总面积的72.75%；桥西区不透水面面积为53.62 km²，占桥西区总面积的79.67%。

3.3 精度验证

本文选取的80个样本区域分布如图 5所示。精度评估结果表明，线性光谱混合分析法提取结果的RMSE为0.246，归一化差值不透水面提取结果的RMSE为0.332，所以对于石家庄中心城区405.8 km²的范围内不透水面的提取而言，线性光谱混合分析法获取的不透水面的结果较好。

通过目视解译，本研究发现归一化差值不透水面指数（NDISI）提取的结果中，在石家庄中心城区周边有部分裸地被误分为不透水面从而导致不透水面被高估，同时存在部分平房被漏分从而导致不透水面被低估；在石家庄中心城区内各个区边缘部分则存在少部分植被被误分为不透水面从而导致不透水面高估的现象。

通过分析，其原因主要为：一方面，不透水面建筑材料的主要成分中多为砂土石，其或多或少地保留了其原材质的光谱特征，使得不透水面与裸土的光谱特征有一定程度的相似性从而产生误分^[19]，导致石家庄中心城区周边范围存在裸地和房屋被误分的现象；另一方面，它作为一种传统的“硬”分类方法，忽略了混合像元的问题，将某一像元全部归为不透水面或透水面，使得不透水面面积估算不够准确。如一个30 m×30 m的像元，经过线性光谱混合分解得到其不透水面丰度为0.6，而用指数法则其不透水面丰度为1，这将导致不透水面面积的高估；而若经线性光谱混合分解得到的不透水面丰度为0.2，用指数法则不透水面丰度多为0，这将导致不透水面面积的低估。因此石家庄中心城区存在少部分植被被误分为不透水面的现象。另外，热红外波段影像的分辨率比较低，可能会加剧中分辨率影像的混合像元问题，而指数法又忽略了混合像元的问题，因此给不透水面提取结果带来了误差^[18]。

3.4 不透水面空间分布分析

根据基于线性光谱混合分析法得到的不透水面的比例统计分析得到：石家庄市中心城区不透水面总面积为316.00 km²，占石家庄市中心城区总面积的77.81%。各个区的统计情况如图 6所示，长安区不透水面面积为105.39 km²，占长安区总面积的74.96%；裕华区不透水面面积为86.26 km²，占裕华区总面积的80.45%；新华区不透水面面积为68.84 km²，占新华区总面积的75.90%；桥西区不透水面面积为56.85 km²，占桥西区总面积的84.47%。

根据上述分析结果，可以看出石家庄中心各城区不透水面所占比例相近，表明石家庄中心城区各个区发展比较均衡。

为了进一步分析石家庄中心城区不透水面的空间分布特征，根据参考文献^[20]的研究成果，采用等间距法将石家庄区域分为极高覆盖度不透水面（ISC≥90%）、高覆盖度不透水面（70%≤ISC＜90%）、中覆盖度不透水面（50%≤ISC < 70%）、低覆盖度不透水面（30%≤ISC < 50%）、极低覆盖度不透水面（10%≤ISC < 30%）、自然地表（ISC < 10%），如图 7所示。从图中可以看出，不透水面主要分布在市中心，长安区和新华区北部不透水面丰度偏低。这表明，石家庄中心城区南部发展较快，石家庄可能进一步向南部扩张。

4 结束语

基于石家庄中心城区Landsat8影像，本文使用线性光谱混合分析法（LMSA）和归一化差值不透水面指数（NDISI）并结合水体掩膜运算，提取了石家庄市中心城区的不透水面，并分析了石家庄市中心城区的不透水面空间分布格局，主要结论如下：

1）利用Landsat8 OLI影像，使用基于纯净像元指数的端元提取方法选择端元，并应用线性光谱混合模型进行混合像元分解，能够有效解决中分辨率影像中的混合像元问题，相比于指数法得到了更高的精度。

2）基于LSMA方法得到的不透水面比例计算得到石家庄中心城区不透水面面积为316.00 km²，占中心城区总面积的77.81%；其中，桥西区、裕华区、新华区和长安区不透水面面积占比分别为84.47%、80.45%、75.90%、74.96%，空间分布相对分散。总体而言，石家庄中心城区各个区不透水面面积比例接近，其中裕华区和桥西区不透水面面积比例最大，表明石家庄中心城区南部发展较快，石家庄有可能进一步向南部扩张。

随着遥感技术的快速发展，针对混合像元分解中不透水面易与土壤混淆的问题，可以利用多源数据融合的方法来解决，比如引入LiDAR的高程数据来区分土壤与建筑物。

随着石家庄城市化进程的加速发展以及进一步推进“海绵城市”建设的需求，城市不透水地面分布信息可为城市规划决策提供重要基础数据。因此，本研究结果可为石家庄城市布局及规划建设提供基础数据，对于石家庄城市规划与环境保护具有重要的参考价值和现实意义。