城市化水平是衡量一个国家和地区经济、社会、文化、科技发展水平的重要标志, 而城市空间范围的扩展是城市化的一个重要衡量指标。尽管城市区域只占地球表面很小的面积(约3%), 然而其扩展变化对局地气候、生态环境以及全球变化影响显著^[1]。

遥感是获取城市范围及其扩展信息的重要手段。国内外学者在利用遥感数据进行城市范围提取和扩展研究方面开展了大量工作, 并且提出不少方法^[2-4]。各种不同分辨率的图像数据用于城市建成区的提取和变化分析^[5-7], 其中, 美国Landsat系列卫星数据(MSS, TM, ETM+和OLI)和美国军事气象卫星计划(Defense Meteorological Satellite Program, DMSP)搭载的OLS (Operational Linescan System)传感器提供的夜间灯光数据, 因获取时间早(均自1972年开始获取数据)、历史数据丰富, 并免费开放获取, 广泛应用于城市范围提取以及城市扩展研究^[8-9]。

Landsat系列卫星数据因其较高的空间分辨率(30 m)、较大的空间覆盖范围以及丰富的光谱信息和空间信息, 适合单个城市或区域城市范围的提取和扩展分析^[10]。但是, 利用Landsat系列数据进行城市建成区提取和扩展分析仍然存在许多问题。在30 m分辨率尺度上, 城市内部不同地物类别(植被和人造地物)之间混合现象普遍, 导致混合像元大量存在。由于不同地物类别在像元尺度混合的比例在空间上变化显著, 因此城市类别的光谱特征差异性很大^[11]; 同时城市类别与周围的地物类别(如裸土)之间存在显著的光谱特征相似性。这些因素给利用Landsat系列数据进行城市建成区的提取与变化分析带来困难。因此, 只利用光谱信息提取城市建成区及探测城市扩展的精度有限。近年来, 不少研究者尝试加入空间信息(如图像纹理)以及综合利用不同传感器的数据来改进城市建成区及其扩展信息的提取精度^[12-14]。

DMSP/OLS夜间灯光数据具有独特的夜间低光探测能力和反映人类活动信息(尤其是城市区域的人类活动信息)的特点, 可以提供大范围的城市用地信息, 因此广泛应用于区域或全球尺度的城市范围提取和变化分析研究^[15-16]。但是, 夜间灯光数据的空间分辨率比较低(1 km), 并且存在溢出效应(overglow)^[17-19]。

尽管Landsat TM/ETM+数据和DMSP/OLS夜间灯光数据提供了互补的信息, 但由于DMSP/OLS夜间灯光数据的分辨率较低, 很少有研究者将夜间灯光数据与Landsat TM/ETM+数据结合起来, 用于提取城市建成区范围及扩展信息。许多研究利用DMSP/OLS夜间灯光数据提取区域尺度的城市范围, 并利用Landsat TM/ETM+来验证DMSP/OLS夜间灯光数据的提取结果^{[15, 20]}。最近, 有研究利用DMSP/OLS夜间灯光数据对Landsat数据的城市提取结果进行后处理和优化, 进行区域城市扩展的提取和分析^[21], 但是, 后处理和优化的过程依赖于较高分辨率的Landsat TM/ETM+数据的城市提取结果。本文的目的是结合DMSP/OLS夜间灯光数据和Landsat TM/ETM+数据以及提取的多波段图像纹理^[12], 提取城市扩展信息, 并进行比较分析。城市地区主要由具有不透水表面的建成区土地所组成。本文中, 城市建成区代表被特定的物理结构直接覆盖的区域, 包括建筑物、街道和不透水层等一切人造结构^{[5, 22]}。

1 研究区及数据 1.1 研究区

本文以我国北方的超大城市天津市为研究区(图 1)。天津市位于华北平原的东北部, 面积为11946.88 km², 是中国第三大城市和直辖市, 区位优势明显。天津市是华北、西北广大地区的出海口, 重要的国际港口城市以及生态城市, 是中国北方对内、对外开放两个扇面的轴心, 是中国北方的经济中心、环渤海地区经济中心, 具有极其重要的社会、经济、生态价值。由于经济的快速发展, 天津市的城市范围也显著扩展。

1.2 数据及预处理

本文采用Landsat TM/ETM+数据和DMSP/ OLS夜间灯光数据, 提取和分析近20年来天津市及其周边的城市扩展信息。选取1992年7月、2001年9月和2006年7月的Landsat TM/ETM+多光谱数据。上述3个时相的图像数据的获取时间均为夏季, 季相基本一致。同时, 选取天津地区上述3个年份的DMSP/OLS夜间灯光数据。此外, 本文还收集了天津市行政区图和由《天津市统计年鉴》得到的3个年份的建成区面积统计数据。

首先对3个年份的Landsat TM/ETM+数据进行预处理。由于2001年的图像有少量云雾覆盖, 对该图像进行去雾霾处理(haze reduction), 去除一部分薄云。对3个时相的TM/ETM+图像进行几何配准, 配准的均方根误差(RMS)在0.4个像元以内。将3幅配准后的图像均裁剪成为4435×2997像元大小, 裁剪后的图像如图 2所示。该地区的主要地物类别包括建成区、农田、裸地、陆地水体、海洋、林地等6类。对每个时相的Landsat图像分类所采用的训练样本和精度评价所用的检验样本见表 1。

由于DMSP/OLS传感器没有进行辐射校正, 不同卫星的OLS灯光数据不能直接对比, 因此, 采用Elvidge等^[23]提出的相互校正(intercalibration)方法, 对不同卫星的OLS灯光数据进行校正。为了便于对两类图像数据所得到的信息进行叠加与对比分析, 将夜间灯光数据重采样为与Landsat TM/ ETM+数据相同的30 m像元大小, 并转换为UTM椭球坐标系统。天津地区3个年份的OLS夜间灯光图像如图 3所示。

2 方法

本文的总体思路可以概括为:将Landsat TM/ ETM+数据与DMSP/OLS数据结合起来提取多年的城市建成区范围, 得到不同时间段的城市扩展情况。为了减少因裸土和城市建成区的光谱相似性导致的误分, 首先利用DMSP/OLS数据得到城市建成区的大致范围, 并且利用这个范围内的光谱数据以及图像纹理信息, 提取城市建成区范围。通过与只利用Landsat TM/ETM+数据光谱信息和利用由DMSP/OLS数据提取的建成区范围和城市扩展信息以及统计数据中相应的建成区面积对比, 评价和验证由Landsat TM/ETM+数据得到的城市建成区扩展结果。总体流程如图 4所示。

2.1 城市范围提取与扩展检测

本文采用分类后比较方法^[24], 从多时相数据提取城市扩展, 即首先分别对不同时相的图像数据进行分类, 然后比较分类结果, 得到城市扩展结果。

尽管Landsat TM/ETM+数据广泛应用于城市范围提取及扩展分析^{[10, 13]}, 但由于城市建成区光谱特征的复杂性^[12]及与周围地物类别(如裸土)的相似性, 只利用光谱信息来提取城市建成区范围及扩展信息的精度有限, 因此将光谱信息与其他信息相结合是改进提取精度的重要策略^[12-13]。为了改进基于Landsat TM/ETM+数据的分类精度, 本文利用DMSP/OLS夜间灯光数据和纹理信息。利用灯光数据的目标是消除灯光区外裸土与城市建成区的混淆, 加入纹理信息是为了减少灯光区内裸土与建成区的混淆。

首先利用DMSP/OLS夜间灯光数据得到城市建成区的大致范围, 然后对这个范围内的Landsat TM/ETM+图像进行进一步分类。由于DMSP/OLS夜间灯光数据反映的是人类活动的特征, 城市建成区和大型居民点等在DMSP/OLS夜间灯光图像上具有较高的灰度值, 在图像上表现为较亮的特征; 裸土和植被等区域在DMSP/OLS夜间灯光图像上的灰度值为零或很低, 在图像上较暗。由于夜间灯光数据的分辨率不高, 范围较小的村庄在图像上基本上不显示, 所以可以通过灰度的差异把城市区域与裸土和植被等区域区分开。因此, 采用一个较小的阈值, 从DMSP/OLS图像上提取大致的城市区域。尽管由于DMSP/OLS灯光数据具有溢出效应, 导致图像上探测到的灯光区面积往往大于实际城市建成区范围, 但对利用灯光数据提取的城市区域范围内的Landsat TM/ETM+图像进行分类, 可以消除灯光区域以外地区裸土与城市建成区的混淆, 在一定程度上改善分类精度。

为了减少灯光区范围内的城市建成区与裸土的混淆, 将图像的纹理特征与原始的光谱信息相结合进行地表覆盖类别分类。由于城市建成区的光谱特征变异性大, 纹理特征比较显著, 而裸土的光谱特征较均匀, 纹理特征并不明显, 因此, 加入纹理信息可减少两类地物间的混淆。

本文采用一种基于多元变差函数(multivariate variogram)的多波段纹理(multivariate texture)特征^{[12, 25]}。利用多波段纹理的一个优点是不需要进行波段选择^{[12, 25]}, 即可以利用多光谱图像的所有波段来计算图像的纹理, 反映的是图像灰度值的多元空间变化; 而利用经典的单波段图像纹理(如经典的基于灰度共现矩阵, gray-level co-occurrence matrix, GLCM)^[26]需要选择一个合适的波段计算图像纹理。多波段纹理特征成功地应用于地表覆盖分类, 尤其是城市地区的地表覆盖分类^{[12, 25]}以及岩性识别和提取^[27]。基于多元变差函数的多波段纹理计算公式为

$ \begin{array}{l} {\gamma _m}(\mathit{\boldsymbol{h}}) = \frac{1}{2}E[{\bf{DN}}(x) - {\bf{DN}}(x + \mathit{\boldsymbol{h}}) \cdot \\ \;\;\;\;\;\;\;\;\;\;\mathit{\boldsymbol{M}}{({\bf{DN}}(x) - {\bf{DN}}(x + \mathit{\boldsymbol{h}}))^{\rm{T}}}], \end{array} $

(1)

DN(x)是像元位置为x处的像元矢量, DN(x+h)是像元位置为x+h处的像元矢量, h是距离矢量; E是数学期望; M为距离函数, 如欧氏距离、马氏距离和光谱角距离等。本文采用光谱角距离来计算多元变差函数纹理特征, 因为该纹理在城市地区的地表覆盖分类中的性能最好^[12]。基于光谱角距离的多元变差函数纹理计算公式为

$ \begin{array}{l} {\gamma _{{\rm{m\_Sa}}}}(\mathit{\boldsymbol{h}}) = \frac{1}{{2N(\mathit{\boldsymbol{h}})}}\sum\limits_{i = 1}^{N(h)} {{\rm{arc}}\cos } \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\left\{ {\frac{{\sum\limits_{k = 1}^p {{\rm{d}}{{\rm{n}}_k}({x_i}) \cdot {\rm{d}}{{\rm{n}}_k}({x_i} + \mathit{\boldsymbol{h}})} }}{{{{\left[ {\sum\limits_{k = 1}^p {{{\left( {{\rm{d}}{{\rm{n}}_k}({x_i})} \right)}^2}} } \right]}^{1/2}} \cdot {{\left[ {\left[ {\sum\limits_{k = 1}^p {{{\left( {{\rm{d}}{{\rm{n}}_k}({x_i} + \mathit{\boldsymbol{h}})} \right)}^2}} } \right]} \right]}^{1/2}}}}} \right\}, \end{array} $

(2)

其中, dn_k(x)是像元x的第k波段的灰度值, dn_k(x+h)是像元x+h的第k波段的灰度值, p是图像的波段数; N(h)是距离矢量为h的像元对数。多波段纹理计算中主要有两个参数:窗口大小和距离矢量长度(步长)^[12], 本文通过试验选出最佳参数。

为了有效地利用所得到的多波段纹理信息和光谱信息进行地表覆盖分类, 本文采用支持向量机(SVM)作为分类算法, 得到不同年份的地表覆盖分类图。支持向量机是一种统计学习方法, 主要通过解算最优化问题, 在高维特征空问中寻找最优分类超平面解决分类问题, 可有效地处理小样本和高维数据^[28]。研究表明, 支持向量机方法能得到比最大似然法和决策树等传统分类方法更高的分类精度^[28]。分类后, 将农田、裸地、陆地水体、海洋、林地等类别合并为“非城市”类, 这样每个年份最终的分类结果只包括城市和非城市两个类别。通过比较两个不同年份的城市建成区提取结果, 得到相应时间段的城市建成区扩展信息。

2.2 结果评价与验证方法

为了全面评价所提出的方法, 本文采用以下几种方法。

1) 比较本文提出的方法与只利用Landsat TM/ ETM+光谱数据的城市范围提取与变化检测精度。通过选取检验样本, 采用混淆矩阵, 计算分类精度及变化检测精度。混淆矩阵通过逐像元比较分类结果与地表真实数据得到, 从中可计算总体分类精度、Kappa系数和各类精度(生产者精度和用户精度^[29]。

2) 采用阈值法从DMSP/OLS夜间灯光数据提取城市建成区范围以及城市扩展信息, 验证本文提出的结合Landsat TM/ETM+数据与DMSP/OLS数据的城市范围与扩展检测方法。阈值法简单, 容易操作, 并可在一定程度减少灯光的溢出问题^[30-31], 是广泛采用的城市建成区范围提取方法^{[17, 30-32]}。通过对灯光影像设定不同的阈值, 将得到的建成区范围与参考图像中的城市范围进行对比, 得到的建成区范围与参考图像最接近时的灯光数据阈值为最佳阈值^[17]。本文采用从较高分辨率的Landsat TM/ ETM+图像提取的城市范围作为参考数据。与已有研究^[33]相同, 对多年的DMSP/OLS数据采用同样的阈值提取城市范围。在利用阈值法得到每个年份的城市建成区范围之后, 通过比较不同年份的城市建成区范围, 得到不同时间段的城市扩展情况。与Landsat TM/ETM+数据和DMSP/OLS数据相结合得到的城市范围与扩展检测结果进行对比, 比较1992-2001年和2001-2006年两个时间段的城市扩展区域, 分析城市随着时间推移的空间分布变化信息。

3) 计算并比较不同方法得到的城市区域面积的年增长率。城市面积年增长率(annual growth rate, AGR)^[34]表达城市在某一时段的面积变化情况, 计算公式为

$ {\rm{AGR }} = {\rm{ }}100 \times \left( {{{\left( {{\rm{U}}{{\rm{E}}_{{\rm{end}}}}/{\rm{U}}{{\rm{E}}_{{\rm{start}}}}} \right)}^{(1/d)}} - {\rm{ }}1} \right), $

(3)

式中, UE_start为在初始年份的城市面积范围, UE_end为在末尾年份的城市面积范围, d为以年为单位的研究时间跨度。本文分别对Landsat TM/ETM+数据与灯光数据相结合所得结果、DMSP/OLS夜间灯光数据所得结果及建成区面积统计数据计算AGR, 并做定量比较。

3 结果与讨论 3.1 结合Landsat TM/ETM+数据和DMSP/OLS数据的城市扩展结果

首先选取阈值为15, 对DMSP/OLS数据提取3个时相天津市建成区的大致范围。在结合纹理信息的分类中, 步长选为1, 分别计算窗口大小为3×3, 5×5和7×7的多波段纹理。通过试验, 发现窗口大小为5×5时分类效果最好。对得到的城市建成区大致范围内不同时间的Landsat TM/ETM+图像和相应的多波段纹理进行分类并合并后, 得到城市和非城市两类。3个时相图像的建成区提取精度见表 2。

从表 2可以看出, 与灯光数据结合的3个时相图像分类总体精度均在93%以上, Kappa系数均在0.86以上, 显著高于只利用Landsat TM/ETM+数据光谱信息的分类精度。例如, 在1992年的结果中, 总体精度(OA)从87.11%提升到94.51%, Kappa系数从74.35%提升到89.03%。从各类别的分类精度来看, 与只利用光谱数据的分类结果相比, 本文所提方法的结果中, 城市类别的用户精度(UA)显著提高, 生产者精度(PA)基本上保持不变; 与此相对应, 非城市类别的生产者精度显著提高, 用户精度基本不变。例如, 在1992年的分类结果中, 城市类别的用户精度从79.82%提高到91.00%, 生产者精度均为98.48%;非城市类别的生产者精度从76.26%提高到90.71%, 用户精度基本上不变。这表明, 结合DMSP/OLS夜间灯光数据和多波段纹理, 大大减少了非城市类别(裸土)与城市类别的混淆, 显著提高了分类精度。

图 5为只利用Landsat TM/ETM+数据光谱信息的分类结果与结合灯光数据和Landsat TM/ETM+数据的分类结果的对比(以1992年图像局部为例)。从图 5可以看出, 只利用Landsat TM/ETM+数据的光谱信息, 存在显著的混淆(图 5(b)); 然而在结合了纹理信息和灯光数据的结果中(图 5(c)), 灯光区域外裸土与建成区的混淆被消除。在灯光区内, 结合纹理信息的分类结果对裸土与建成区的混淆现象也有所减少(如图 5(a)中红圈区域, 在图 5(b)中误分为建成区的裸土区域, 在图 5(c)中正确分为裸土), 建成区范围提取结果更加准确。其他两个年份的结果与此相似。上述结果验证了方法的有效性。

图 6为与灯光数据结合的分类结果。从图 6可看出, 天津市的范围随着时间变化显著扩展: 1992年时, 中心城区和滨海新区的面积均比较小, 且二者之间没有明显的连接区域; 2001年分类结果显示中心城区呈蔓延式扩展, 滨海新区已经显著发展起来, 与中心城区显示出互相连接的趋势; 2006年, 中心城区进一步向周边地区扩展, 滨海新区也有所扩大, 已经与中心城区沿着二者之间的海河连接带明显相连。

表 3是只利用Landsat TM/ETM+数据光谱信息与本文方法的城市扩展信息检测精度评价对比。可以看出, 与只利用光谱数据的城市扩展检测结果相比, 本文方法城市扩展检测精度更高, 如1992-2001年的城市扩展检测总体精度(OA)从84.02%提高到88.40%, Kappa系数从57.61%提高到67.33%;各类别的变化检测精度也均有所提高。本文提出的方法在两个时间间隔的城市变化检测精度均为87%以上, Kappa系数均大于67%。

图 7为结合3个时相的Landsat TM/ETM+数据和灯光数据得到的不同时间段天津的城市扩展结果。可以看出, 1992-2001年扩展的部分主要分布在天津市中心城区和滨海新区, 以二者为中心呈辐射式扩展, 并沿二者连接部分有一定扩展; 2001-2006年中心城区和滨海新区进一步扩展, 呈现内部紧密扩展与外围稀疏扩展相结合的趋势。

3.2 结果对比

采用阈值法从DMSP/OLS数据提取不同时间的天津市城市范围。通过试验, 发现阈值为60的结果最好。图 8是由DMSP/OLS夜间灯光数据得到的不同时间段天津城市扩展情况。从图 8可以看出, 1992-2001年天津市主要围绕中心城区和滨海新区扩展, 2001-2006年, 在中心城区和滨海新区持续向外扩展的同时, 二者间的连接区域也逐渐加宽, 与结合灯光数据的Landsat TM/ETM+数据城市扩展提取结果一致。

表 4是由Landsat TM/ETM+数据结合灯光数据得到的建成区面积、由DMSP/OLS夜间灯光数据得到的建成区面积和建成区面积统计数据, 以及相应的年增长率。从年增长率来看, 利用本文提出方法提取的建成区面积年增长率与统计数据年增长率基本上一致, 由DMSP/OLS夜间灯光数据得到的城市年增长率则均偏大, 这可能是由DMSP/OLS夜间灯光数据溢出效应导致的。我们可以看到, 由Landsat TM/ETM+数据结合灯光数据得到的建成区面积、DMSP/OLS夜间灯光数据得到的建成区面积和建成区面积统计数据三者反映的城市扩展趋势是一致的, 即2001-2006年的城市年增长率均高于1992-2001年, 表明2001-2006年城市的加速扩展。这些结果说明结合Landsat TM/ETM+数据和灯光数据所得到的1992-2006年天津市城市范围变化是可信的。从建成区面积来看, 本文提出的方法以及由灯光数据得到的城市建成区面积均高于建成区面积统计数据, 原因可能为建成区面积统计数据统计的只是市区面积, 本文方法提取的建成区以及夜间灯光数据提取的建成区包括天津市辖区范围内的县和镇, 即大的居民点。

前面的分析验证了本文提出方法的有效性。因此, 虽然本文以天津市为研究区, 但结合Landsat TM/ETM+数据与DMSP/OLS夜间灯光数据的方法可以应用于其他大型城市建成区范围提取和扩展检测的研究中。

4 结论

本文将多时相的Landsat TM/ETM+数据与DMSP/OLS夜间灯光数据相结合, 提取天津市的城市范围与扩展信息, 并对结果进行了对比。结果表明, 与只利用Landsat TM/ETM+的光谱数据相比, 利用夜间灯光数据并加入多波段的图像纹理, 可减少裸土与建成区的混淆情况, 有效提高城市建成区提取精度和扩展检测的精度。通过与只利用Land-sat TM/ETM+数据光谱信息和利用DMSP/OLS夜间灯光数据提取结果及建成区面积统计数据对比, 验证了本文所提出方法的有效性。对天津地区来说, 1992-2006年, 城市建成区面积逐渐增加, 且在大致2000年之后扩展速率明显加快。从空间分布来看, 城市建成区的扩展主要分布在中心城区、滨海新区以及二者之间的连接地区。因此, 综合利用Landsat TM/ETM+数据和DMSP/OLS夜间灯光数据提取和分析城市扩展, 可得到可靠的城市建成区扩展趋势, 是一种大型城市建成区范围提取和扩展信息研究的有效方法。