不透水面定义为诸如屋顶、沥青或水泥道路以及停车场等具有不透水性的地表面; 单位地表面积中不透水面的面积所占的百分比, 定义为不透水面覆盖度(impervious surface percentage, ISP)^[1]。平原河网地区的水文地貌特征较为特殊, 河道纵横交错, 河流比降小, 不透水面的扩增将会加快地表径流速度, 缩短汇流时间, 增加洪水发生的频率; 除此之外, 积蓄在不透水面上的大量污染物随地表径流进入河流, 造成河流污染, 减少健康水的供应, 成为影响平原河网地区河流生态系统的主要因素。

在众多ISP估算的模型中, 基于分类与回归树(classification and regression tree, CART)集成学习的估算方法继承了一般决策树具备的所有优点, 既可以用于分类研究, 又能进行连续变量的预测和回归, 且实现简单, 运算效率高, 成为研究不透水面的热点。然而, 决策树是一种弱学习算法, 分类精度的提高取决于树结构的合理构建与剪枝处理。Nathaniel等^[2]、马雪梅等^[3]虽然综合利用了遥感图像的多源特征, 但并未对多源特征进行选择; 廖明生等^[4]将Boosting技术引入到CART算法中, 但仅利用了遥感数据的光谱特征, 未充分利用遥感数据的纹理特征和空间变换特征。 以上两种方法在训练样区存在大量噪声的情况下, 估算精度会大大降低; YANG等^[5]通过对不同特征组合构建了多种估算模型, 对每种模型都进行试验, 选取误差最低和相关系数最高的模型作为最优模型, 这种方法虽然一定程度上避免了数据的冗余性, 但需要更多的先验知识和专家知识。

针对以上情况, 本文引入目前在知识挖掘领域采用的变精度粗糙集算法, 删除冗余特征, 求出最优的ISP估算特征集, 并采用TVDI对初始结果进行后处理。

研究区遥感图像选用2010年5月24日采集的Landsat TM图像, 全景图像质量较好, 清晰无云; 也选取部分区域ALOS图像, 用于获取ISP训练和验证数据。将ALOS多光谱10 m和2.5 m分辨率图像融合, 生成3 m分辨率数据, 成像时间均为2010年5月3日。经几何配准将TM数据和ALOS 数据统一到UTM WGS84投影坐标系下; 由于大气效应对ISP估算并没有明显影响^[6], 所以仅对其图像进行辐射定标, 将DN值转化为反射率。

研究区如图1所示, 图上右图为TM图像和样本区ALOS图像, 均由4, 3, 2波段合成。该地区地形平坦, 河网稠密, 湖群广布。随着城市化建设的迅速发展, 市区不断向郊区延伸扩展, 河湖数量和面积都在不断减少, 原有的自然水系遭到不同程度的人为破坏, 河流自然的调蓄功能萎缩。

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 研究区域Fig.1 Map of study area

首先, 利用高分变率遥感图像获取ISP的训练数据和验证数据; 然后, 获取CART模型的属性变量并对属性变量进行约简; 最后, 对中分辨率的Landsat TM遥感数据进行ISP估算和制图。

2.1 ISP训练数据和验证数据获取

采用最大似然监督分类, 将融合生成的ALOS图像分为不透水面和透水表面两大地类, 前者包括建筑和道路两个子类, 后者包括耕地、植被和水体3个子类。从该分类结果中, 计算落在以每个3 m分辨单元为中心的10像元× 10像元窗口内的不透水面像元总数百分比, 并重采样到30 m分辨率(如图2所示), 其中云区和水体不包含在百分比计算内。

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 样本区ISP分级结果Fig.2 ISP classification results of sample area

为满足各种透水面和不透水面地物的样本选取, 生成60像元× 60像元的规则格网, 如图3所示, 取格网的中心点共3 600个样本点, 对其进行排查, 修改误差较大的点, 最后将样本点分为两个相互独立的子集, 其中, 对应ISP估算模型的训练样本2 600个, 对估算结果进行验证的样本1 000个。

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 样本区规则格网Fig.3 Regular square grid of sample area

2.2 CART模型属性变量获取

综合考虑平原河网地区流域特点和植被-不透水面-土壤模型^[7], 由于平原河网地区水体和耕地面积比例较大, 裸露土壤较少, 因此将每一个像元视为不透水面、植被、水体3种代表性土地覆盖类型的线性组合。国内外学者针对这3种地物已研究发展了几十种不同的指数模型, 本文在此基础上, 选取了比较常用的指数模型。针对不透水面选取了归一化建设用地指数(normalized difference built-up index, NDBI), 城市指数(urban index, UI), 基于指数的建筑用地指数(index based built-up index, IBI); 针对水体提取了归一化水体指数(normalized difference water index, NDWI), 改进的归一化水体指数(modified normalized difference water index, MNDWI), 混合水体指数(combined index of NDVI and MIR for water body identification, CIWI); 针对植被提取了比值植被指数(ratio vegetation index, RVI), 垂直植被指数(perpendicular Vegetation index, PVI), 归一化差值植被指数(normalized difference vegetation index, NDVI)和土壤调节植被指数(soil adjust vegetation index, SAVI)。

纹理也是遥感图像的重要信息, 本文选择常用的灰度共生矩阵法, 对TM图像主成分变换生成的第一主成分, 衍生8个纹理特征值, 分别包括均值(mean)、方差(variance)、协同性(homogeneity)、对比度(contrast)、非相似度(dissimilarity)、熵(entropy)、角二阶矩(second moment)和相关度(correlation), 其中纹理特征的窗口大小选择3像元× 3像元, 移动步长取1, 移动方向为0。除此之外, 还选择了最小噪声分离变换(minimum noise fraction rotation, MNF)的前4个波段信息; 以及穗帽变换的前3个分量: 亮度分量(KT-b)、绿度分量(KT-g)和湿度分量(KT-w)。

综上所述, 本文获取CART模型的属性变量共31个, 如表1所示。

表1

Tab.1

表1(Tab.1)

表1 CART模型属性变量 Tab.1 Attribute variable of CART model analysis 代码 类型 代码 类型 B1 TM1 B12 MNDWI B2 TM2 B13 CIWI B3 TM3 B14 NDBI B4 TM4 B15 UI B5 TM5 B16 IBI B6 TM7 B17 MNF1 B7 RVI B18 MNF2 B8 PVI B19 MNF3 B9 NDVI B20 MNF4 B10 SAVI B21 KT-b B11 NDWI B22 KT-g B23 KT-w B28 非相似度 B24 对比度 B29 熵 B25 方差 B30 角二阶矩 B26 协同性 B31 相关度 B27 对比度

表1 CART模型属性变量 Tab.1 Attribute variable of CART model analysis

2.3 CART模型属性变量约简

以上生成的31个属性变量间的相关性是未知的, 这一缺点将导致决策树中子树的重复, 降低分类的效率。随着分类数的增加, 叶子数和节点数将增多, 使决策树变得复杂, 掩盖有用的规则信息。Ziarko^[8]等人提出的变精度粗糙集理论与决策树具有很强的优势互补性, 可以有效地消除数据冗余; 与标准粗糙集相比, 变精度粗糙集通过设置阈值参数β , 放松了标准粗糙集对近似边界的严格定义, 0.5< β ≤ 1, 随着β 的增大, 近似边界变窄, 即变精度粗糙集意义下的不确定区域变小, 对数据的不一致性有一定的容忍度, 容忍度的大小可由阈值β 确定。本文在参考文献^[9]的基础上以变精度粗糙集的分类质量的量度对条件属性进行选择, 属性约简基本流程如下:

1)将遥感图像视为一个信息系统, 将生成的2 600个样本点作为对象, 样本的31个特征作为条件属性, 对应的ISP值作为决策属性, 构成一个信息决策表;

2)利用Rosetta软件中的NaiveScaler算法^[10]对信息决策表的条件属性进行离散化;

3)计算条件属性和决策属性的等价集, 以及下β 近似集和上β 近似集;

4)决策表属性约简, 用数据集U初始化根节点R, 用根节点R初始化队列Q, 取出队列中的第一个节点N, 如果N不纯, 估计每个节点上的β 分类质量的量度, 选出β 分类质量的量度最大的属性, 根据属性的不同取值将N分裂为N₁, N₂, …, N_|N|;

5)将根据选为结点的属性所分裂的N₁, N₂, …, N_|N|输入队列, 转向步骤4);

6)检查决策表的一致性, 即从决策表中将属性集的某个属性删除后, 若决策表中不出现新的不一致则该属性可以约去, 否则该属性不能约去。重复步骤4), 5), 直至没有属性冗余。

运用上述约简算法对离散化后的数据进行求解, 具体步骤见参考文献^[11-12], 分别取β 的经验值0.667和0.75, 得到2组不同的属性约简值, 作为CART树的独立变量对研究区进行不透水率估算, 决策树信息和分类精度如表2所示。

表2

Tab.2

表2(Tab.2)

表2 β 约简信息表 Tab.2 β reduction information β 取值范围 属性约简集合 分类质量 决策树的结点数 ISP估算精度(R2) [0.5, 0.667] {B1, B7, B12, B14, B15, B16, B18, B19, B22, B24, B25, B27, B31} 0.84 41 0.89 [0.5, 0.75] {B1, B4, B7, B11, B12, B14, B17, B19, B20, B21, B22, B24, B25, B28, B30, B31} 0.64 65 0.68

表2 β 约简信息表 Tab.2 β reduction information

如表2所示, β 取0.667时, 31个属性特征被约简为13个, 并取得较高的样本分类质量和ISP估算精度, 决策树的结点数仅为41; 当β 取0. 75时, 分类质量和估算精度均降低, 决策树的结点数反而增多。而未经变精度粗糙集属性约简的31个属性变量构造的决策树结点数达到426个, 可见当β 取0.667时, 得到的13个属性构建的决策树的结构简单, 大大降低了决策树的复杂度, 并获得较高的估算精度, 因此选择此阈值下得到的属性组合为最优特征集。

将利用ALOS图像获得的2 600个ISP训练样本作为CART模型中的目标变量, 经变精度粗糙集约简的13个属性特征作为属性变量, 估算整个实验区的ISP, 为验证变精度粗糙集属性约简的有效性, 采用相同的训练样本对单一CART树算法即未经属性约简的31个特征进行了ISP估算实验。采用上文形成的1 000个验证样本作为ISP参考值, 分别对2种方法的估算结果进行验证, 验证指标选择了ISP估算值和ISP参考值的均方差(MAE)、均方根误差(RMSE)和决定系数(R²)来评价, 如表3所示。

表3

Tab.3

表3(Tab.3)

表3 ISP估算精度比较 Tab.3 Accuracy Comparation of ISP estimation 估算方法 MAE/% RMSE/% R2 单一CART算法 18.9 25.1 0.67 集成CART算法 9.5 11.9 0.89

表3 ISP估算精度比较 Tab.3 Accuracy Comparation of ISP estimation

由表3可知, CART集成算法的精度明显优于单一CART算法, 均方差和均方根误差分别降低了9.4%和13.2%, R²提高了0.22。图4给出了两种方法估算的ISP与参考值间线性拟合结果。

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 ISP估算值与参考值散点图Fig.4 Scatterplot of estimating and reference ISP

从图4可以看出, 单一CART算法的估算结果与参考值散点图的点分布较分散(图4(a)), 拟合直线较多地偏离了直线y=x(黑色虚线), 其斜率和截距分别为0. 53和43.03, ISP高估现象较为严重, 大多数点聚集在y=x直线的上方; 而CART集成算法的ISP估算结果与参考值的线性拟合要优于单一CART算法, 散点图(图4(b))的点基本都聚集在对角线附近, 其斜率高达0.86, 截距仅为4.62, ISP低值高估的现象得到改善, 但在高值区还存在低估现象, 拟合线较多地偏离了直线y=x。通过实验分析, TVDI^[13]与ISP存在较好的相关性, 以此关系寻找后处理规则消除ISP高估现象。由于高值区主要出现在市区, 本文分别以苏、锡、常3个市中心点向外10 000 m缓冲区, 在缓冲区内随机获取600个像素点的ISP和TVDI, 并将ISP统一到0~1内, 与TVDI做散点图, 如图5所示。

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 TVDI与ISP散点图Fig.5 Scatter plot of TVDI and ISP

从图5可以看出, TVDI和ISP存在较好的相关性, 随着ISP的增大, 拟合线逐渐偏移到直线y=x的上方, 说明随着ISP的增大, TVDI值略高于ISP值。由此可得到后处理规则为: 在缓冲区内, 当ISP≤ 0.7且TVDI> 0.7时, 该像元的ISP值为以该点为中心5像元× 5像元窗口内的ISP≥ 0.7所有像元的均值。经此规则后处理后的ISP结果再与1 000个验证样本做精度评价, 散点图如图6所示。

图6

Fig.6

	Figure Option ViewDownloadNew Window
	图6 处理后ISP估算值与参考值散点图Fig.6 Scatter plot of estimating and reference ISP after postprocessing

由图6可知, 虽然处理后R²没有变化, 但斜率增大到0.96, 均方根误差减小到10.01%, 拟合线几乎与直线y=x重合, ISP高值低估现象得到明显的改善。研究区ISP估算结果城区(无锡市)放大图见图7。

图7

Fig.7

	Figure Option ViewDownloadNew Window
	图7 研究区ISP估算结果城区放大图Fig.7 Estimating ISP results of the urban area

从图7可以直观地看出, 无锡市城市空间分布呈现出不均匀特征, 主要包括集中连片分布和轴向分布两种典型的分布模式。苏— 锡— 常城市中心线自西北向东南, 由惠山区到北塘区、崇安区、南长区和无锡新区, 经过了无锡国家高新区和无锡新加坡工业园, 该地区靠近京杭运河和沪宁铁路干线, 并且312国道从中穿过, 水陆交通极为发达, 经济发展迅速, ISP较高, 均在50%以上; 城市内的绿地以及城乡结合部的农田区和植被区, 虽然植被覆盖较好, 不透水率有所降低, 但由于受到周围不透水地物的阻隔, 斑块面积小而破碎, 不透水性仍要高于其他农田和植被区的ISP, 基本上处于30%~50%之间, 如图中黄色斑块; 最低的是植被覆盖度较好且周围没有大面积不透水地物分布的农田区和植被区, ISP基本上都在30%以下。

本文针对利用多源特征进行ISP估算存在属性冗余性的不足, 将变精度粗糙集理论引入到CART分析中, 最终获取相关性小且能满足分类精度的属性变量, 并在此基础上利用TVDI阈值对初始分类结果进行后处理。

1)通过精度验证可知, 该方法有效可靠, ISP估算精度明显优于单一CART算法, 均方根误差降低了13.2%, R²提高了0.22, 有效地改善了单一CART算法ISP估算结果中存在的低值高估现象。

2)针对CART集成算法估算结果中存在的高值低估现象, 利用TVDI构建后处理规则对其进行改善, 使ISP估算值与参考值拟合线的斜率由0.86提高到0.96, 截距由4.62降低为2.45, 进一步提高了ISP的估算精度。

3)空地如农田中的土壤和城市内的未利用地存在一定的光谱混淆, 两者较难区分, 空地ISP的估算可能会出现虚高现象; 以及CART集成算法对样本数量的敏感性, 还需进一步的研究解决。