浑善达克沙地榆树疏林的高分辨率遥感识别方法

薛传平， 高志海， 孙斌， 李长龙， 王燕， 张媛媛

(中国林业科学研究院资源信息研究所，北京 100091)

0 引言

浑善达克沙地是内蒙古4大沙地之一，是中国京津冀经济区重要的生态屏障。沙地榆树疏林是一种极具区域特色的植被类型，由于它主要分布于温带地区，并且与分布于热带、亚热带地区的稀树草原(savanna)非常相似，一些研究就称之为温带savanna[1-3]。榆树疏林是浑善达克沙地生态系统的重要组成部分，它的退化会直接导致沙地生态系统的退化和荒漠化的扩张，因此保护好榆树疏林对于防治沙地荒漠化扩张很有必要[4]。准确的榆树疏林覆盖数据是其他科学合理研究的前提条件，其准确的监测结果是管理决策者重要参考依据。由于常规的野外树木调查入库花费人力、财力很大，寻找一种替代的收集榆树疏林覆盖数据的方法很有必要。近年来，遥感技术发展十分迅速，影像空间分辨率有了很大提高，高空间分辨率遥感影像的可选择性越来越高，遥感已成为一种分析大范围林地生态系统的重要工具。用高空间分辨率遥感影像提取树冠轮廓，可有效提高调查效率，减少外业工作强度，对摸清榆树疏林的规模和分布状况，从而制定合理的保护管理对策有实际意义。但如何有效提取树冠仍是一个难题，也是近几年来国内外遥感领域研究的热点。

目前，国内对榆树疏林的研究多局限于生态方面[5-6]，运用遥感技术对榆树疏林树冠提取的研究较少。面向地理对象影像分析技术(geographic object based image analysis，GEOBIA)是将影像分割为有意义的地理对象，利用对象的光谱、几何和纹理等信息，发展与应用自动化或半自动化遥感影像解译分析的理论、方法和工具，提高影像解译准确度及效率，节省人力、物力和财力，被认为是处理遥感中各种复杂特征识别问题的一种有效方法[7]。GEOBIA已在非洲、澳大利亚和美国等稀树草原树冠提取上得到很好应用。Laliberte等使用QuickBird数据，采用多尺度分层分割和CART决策树方法提取干旱区灌木[8]；Gibbes等使用IKONOS数据，基于面向对象分类的方法提取稀疏树木[9]；Boggs使用SPOT-5和QuickBird数据，利用基于像素的归一化植被指数(normalized difference vegetation index，NDVI)阈值方法和面向对象分类方法分别提取研究区稀疏树冠，结果表明面向对象方法要比阈值提取结果精度普遍高[10]；施敏燕等基于QuickBird数据，采用面向对象多层次分割方法提取额济纳胡杨林树冠信息，结果表明该方法可准确示出植株的分布情况[11]。

国外关于Savanna的研究可为榆树疏林的研究提供有益的参考价值，本文探讨使用国产高空间分辨率数据利用遥感技术提取榆树覆盖的方法。该方法将很大程度降低外业调查的工作量，拓宽国产高空间分辨率影像在疏林草原生态系统上的应用前景。

1 研究区概况及数据源

1.1 研究区概况

研究区位于内蒙古自治区正蓝旗，该区为天然沙地榆集中分布区，以“蓝旗榆”著称，其北部为浑善达克山地腹地，南部为低山丘陵，该区地理坐标为：E115°00′～116°42′，N41°56′～43°11′。正蓝旗气候属于干旱、半干旱大陆性季风气候，年平均气温为1.7 ℃，年平均降水量为350 mm左右，年平均蒸发量为2 700 mm左右，呈现“收入少，支出多”的特点。研究区属于欧亚草原区，植物物种丰富，是一个含有散生沙地榆树林的温带型干旱地区，榆树疏林群落结构和组成均单一，很少有其他树种。根据锡林郭勒盟森林资源统计，1976年正蓝旗内榆树疏林面积达2 482 hm2[12]。沙地榆树多呈丛状分布，每丛3～5株，最多达20余株，根据树龄不同可分为5 a以下幼苗、5～30 a幼树、30～70 a成熟树及70 a以上老龄树，最大高度约9～10 m，最大胸径约90 cm，最大冠幅约9～12 m[13]。沙地榆树在每年5月初开始萌动，于10月底落叶，另外考虑到8月底—9月初是该地区打草季节，为减少草地对榆树信息提取的影响，本研究首选影像成像时间是9月中旬—10月底或5月底—7月初。研究区地理位置如图1所示。

图1 研究区地理位置示意图

1.2 数据源

1)遥感数据及预处理。本研究选取了一景GF-2数据，获取时间为2016年9月16日。数据的预处理包括辐射校正、几何纠正和影像融合等。辐射校正是在ENVI5.3软件中FLAASH模型下进行，以消除因辐射误差而引起的影像畸变；几何纠正是在ERDAS9.2软件下进行，参考影像是空间分辨率为2 m高精度正射校正后ZY-3数据，几何纠正的绝对误差小于1个像元。对几何纠正后影像进行投影转换为WGS84 UTM Zone 50N。影像融合算法为NNDiffuse Pan Sharpening，该算法能很好地保留影像的色彩、纹理和光谱信息。GF-2卫星传感器具体参数如表1所示。

表1 GF-2卫星传感器参数

2)地面调查数据。自2011年开始连续6 a对浑善达克沙地及其周边开展野外调查，获得了大量野外样地实测资料。其中，2014—2016年间重点对正蓝旗境内榆树疏林进行3次地面调查，设置大小为50 m×50 m样地，调查内容包括样地内每木经纬度、胸径、冠幅(包括东西、南北2个方向)和树高，同时拍摄样地照片作为补充。本研究所用的GF-2影像区域内含8个实测榆树样地，140棵单木。为弥补GPS定位误差，结合外业调查数据与目视解译判读结果，在融合后影像上目视手动勾绘150个图斑作为参考数据，用于榆树区粗提取过程阈值及条件的设置，其中有24个榆树群和126棵榆树单木。

2 研究方法

为提高识别精度和数据处理速度，本研究以融合后的GF-2影像为数据源，先建立NDVI阈值初步判断，再采用GEOBIA方法深度分析研究区的典型沙地榆树，这对北方干旱、半干旱地区典型树种识别具有重要参考意义；同时对评价该地区的生态环境也具有重要的现实意义。

2.1 榆树区粗提取

由于研究区榆树大多分布在沙质土壤上，并且分布稀疏，可用NDVI初步掩模非榆树区。NDVI能够很好地区分植被区与非植被区，但阈值太小会保留大面积林下植被，阈值过大会排除一些小榆树。树冠反射率受光合物质叶子和非光合物质枝干的共同影响，与其他植被类型相比，尤其是草地，它们有较明显的光谱变化。因此，基于影像对象的空间和光谱信息提取榆树，可通过对NDVI设置阈值生成对象，使用对象的NIR标准差作为区分树冠的一个指标。为防止大面积背景参与分类，结合统计数据，设置面积阈值大于20 m2且小于3 000 m2的限制条件。单棵榆树或小面积榆树群的形状接近椭圆形，利用其圆度和椭圆度几何特征可以很好地将其识别。NDVI计算公式为

(1)

式中：RNIR为近红外波段的反射率；RRED为红光波段的反射率。

人工设置NDVI阈值及条件A，粗提取榆树分布区。NDVI阈值从最小值0.1开始，以0.05为逐步增量，直到NDVI达最大值0.4为止。只有满足条件A，即面积介于20～3 000 m2之间、NIR标准差大于450，且椭圆度大于0.4的对象为潜在榆树分布区P。否则，小于NDVI最小值的对象为其他地类，这里统一归并为背景B，此过程是迭代进行的。图2中B为背景，P为潜在榆树区，条件A为面积、NIR标准差和椭圆度3个限制条件的交集。具体榆树区粗提取过程如图2所示。

图2 榆树粗提取过程

根据统计目视勾绘的参考样本特征值设置NDVI阈值及条件A。本文参考数据特征值统计结果见表2。

表2 参考数据特征值统计

①：b为NDVI均值；n为对象内像素数量；bi为NDVI值；σ为NIR标准差；c为NIR均值；ci为NIR值；a为面积；p为像素大小；r为圆度；s为最小包围对象椭圆半径；l为对象内最大封闭椭圆半径。

2.2 基于GEOBIA方法的榆树区准确提取

本研究使用GEOBIA方法对粗提取结果进一步准确提取。此方法从提出到现在10余a内，已广泛应用于城市监测、灾害监测等很多领域，但在干旱区植被提取方面应用很少[14-15]，其优势是以对象为处理单元克服像素级分类的椒盐现象，且能够充分利用影像对象的光谱、几何和纹理等特征，但特征选择和规则集构建的不确定性是制约其发展的关键因素。因此，本研究基于分离阈值(separability and thresholds，SEaTH)算法优选特征及自动构建规则集。该算法由Nussbaum等人提出，是一种基于J-M(Jeffries-Matudita)距离评价2个类别在某特征上的关联程度和根据高斯概率分布模型计算特征阈值的方法[16]。

1)分割方法。本文采用的分割方法为Multiresolution segmentation 多尺度算法。它是一种自下而上的区域生长算法，通过合并相邻像元或者对象，保证对象内部像元之间同质性最大、异质性最小。多尺度分割主要设置参数包括分割尺度、波段权重及匀质性因子。分割尺度的设置至关重要，如何选择最佳分割尺度也是很多学者研究的热点，目前已出现了一些确定最佳分割尺度的方法或工具，如最大面积法、目标函数法和面积比均值法等。但很多研究仍采用多次试验比较分割结果的方法来确定最佳分割尺度，因为随着分割尺度增大分割对象也增大，但是两者之间并没有直接关系，因此，仍需要多次尝试来确定分割尺度。

2)特征筛选。目视查看潜在榆树区发现，除榆树外还包括灌木、草地、草本湿地、人工杨树林和沙地5种类别。根据样本的代表性和分布均匀性，本文选取470个分割对象作为样本，其中灌木47个、草地52个、草本湿地78个、榆树152个、沙地108个和人工杨树林33个。分别统计其特征值，估计其概率分布，计算2种类别C1和C2之间的J-M距离来评价其可分离度[17]。J-M距离取值范围为[0，2]，其值为0时表示2个类别在某特征上完全不能区分，其值趋向于2时表示2个类别在某特征上的分离度较大，因此选择J-M距离较大的特征组成最优特征集。具体计算公式为

J=2(1-e-B)

(2)

(3)

式中：B为巴氏距离；m1和m2表示类别的某特征值；σ1和σ2表示类别的某特征标准差。

影像对象常用特征包括光谱特征、几何特征和纹理特征。本文考虑的光谱特征包括各波段及NDVI的均值、标准差、比率、最大方差及亮度等；几何特征包括面积、长度、宽度、长宽比、圆度、椭圆适合性和矩形适合性等；纹理特征包括各个波段图像灰度共生矩阵(gray-level co-occurrence matrix，GLCM)的匀质性、反差、熵、均值、标准差、相关性、相性和角二矩阵等。

(4)

(5)

式中n1和n2是类别C1和C2的样本数目。由式(4)和(5)计算出2个阈值T1和T2，选取位于m1和m2之间的阈值作为最佳阈值T。但分割时会存在一些类别之间因J-M距离值偏低而不易区分，若此时直接使用T作为最佳阈值分割结果会不准确。针对这个问题，Marpu等[18]提出阈值调整规则，T′为调整后阈值，J为J-M距离值，调整规则如下：

如果J<0.5，那么忽略该特征；

如果0.5≤J≤1.25，那么T′=m2；

如果1.25

如果J≥1.75，那么T′=T。

4)分类识别方法。根据SEaTH算法筛选的特征及阈值结果对榆树类别进行类描述，使用eCogni-tion中规则分类的classification算法进行榆树分类识别，将其他地类归为背景; 最终将影像分为榆树和背景2个类别。

3 结果分析

3.1 榆树区粗提取

局部榆树粗提取结果如图3所示，其中，蓝色为背景，红色为粗提取榆树分布区，可以发现大多数非榆树区被归为背景类中，榆树区被初步识别出来。目视查看榆树粗提取效果，发现冠幅小的榆树单木能够被很好地识别，并且提取的冠幅接近其真实冠幅，而聚集分布的多株榆树除了榆树之外，其林下植被、沙地和灌丛等也会被粗提取出来。此外，粗提取结果还包括一些草本湿地、人工杨树林等非榆树区，仍需要进一步准确提取。

(a) GF-2融合影像 (b) 粗提取结果

3.2 多尺度分割

利用NDVI阈值粗提取榆树分布区只是一个初步的提取过程，提取结果仍包括一些灌丛、草地、草本湿地、人工杨树林及沙地等非榆树区，需进一步对粗提取结果进行准确提取，首先是影像分割。本文选取多尺度分割方法，分割尺度及匀质性因子是通过多次试验比较分割结果来确定。由于分割过程中主要利用影像的光谱因子，所以光谱因子相对更为重要。波段权重的设置参考Ren等[19]提出的最佳指数因子(optimum index factor，OIF)，该方法是一种简单有效筛选波段的方法，通过计算各波段的标准差和波段间的相关系数的比值选择最佳的波段组合，其具体计算公式为

(6)

式中:Si为第i波段的标准差；Rj为任意2个波段的相关系数。OIF值越大表示其波段组合信息量越大。

通过计算GF-2影像波段间最佳指数因子筛选出最佳波段组合，进而设置多尺度分割的波段权重。本文GF-2影像各波段间相关系数及OIF值分别见表3和表4。

表3 GF-2各波段的标准差及波段间相关系数

①：B4为近红外波段，B3为红光波段，B2为绿光波段，B1为蓝光波段。

表4 GF-2各波段组合OIF值

需要注意分割尺度不易设置过小，防止单棵树冠被过度分割。经过多次尝试，最终确定分割尺度为40，各波段权重分别为0(B1)、1(B2)、2(B3)、1(B4)、1(NDVI)，光谱因子为0.7，形状因子为0.3，光滑度和紧密度为0.5。局部分割效果如图4所示。

(a) 原始影像 (b) 分割效果

3.3 基于SEaTH算法的榆树疏林提取规则集构建

SEaTH算法可根据统计的样本特征值，自动筛选出区分两两类别之间的最优特征及其对应的最佳特征阈值，然后根据筛选结果构建提取某地类的规则。本文主要目的是提取榆树的空间分布，不考虑其他地类，因此重点是构建提取榆树的规则。表5为榆树与其他地类区分的规则，其中，榆树与草地、人工杨树林、沙地之间的J-M距离值均大于1.75，尤其是与沙地间J-M距离值大于1.9，表示榆树与这些类别之间容易区分，特征阈值直接取阈值T。而榆树与灌木、草本湿地之间J-M距离值约1.6，表示榆树与灌木、草本湿地之间易混淆，此时，根据Marpu等[18]提出的阈值调整规则，特征阈值取(T+m2)/2。

表5 榆树的提取规则

①：RatioG为绿光波段比值;StandarddeviationNIR为近红外波段标准差;GLCMEntropy(90)为灰度共生矩阵熵特征;RatioR为红光波段比值;HSITransformationIntensity(R=r，G=g，B=b)为HSI变换强度。

3.4 榆树分类识别结果及精度评价

根据表5采用SEaTH算法构建的榆树准确提取规则，进一步掩模非榆树区，将能够同时满足其中5条规则的对象被识别为榆树，其他地类则为背景。

图5为研究区局部地区榆树覆盖信息识别提取结果，选取2种不同生境(沙丘、塔拉)榆树种群并对比其提取结果，其中蓝色为背景，红色为最终提取的榆树。

(a) 沙丘生境榆树影像 (b) 识别结果

(c) 塔拉生境榆树影像 (d) 识别结果

本研究使用混淆矩阵对分类识别结果进行精度评价。在榆树识别结果上对榆树与背景2个类别分别生成随机点各300个，然后在融合后的GF-2影像上目视判读各个点的类别作为参考，通过对比榆树分类识别结果，逐一判断这600个随机点的分类正误，构建混淆矩阵进行评价，结果见表6。

表6 榆树识别精度评价

本方法提取的总体精度达88.17%，Kappa系数为0.76，其中，背景的用户精度达99.33%，榆树的制图精度达99.14%。研究区榆树种群复杂，不同的发育阶段和生境对榆树疏林识别有不同的影响。其中，冠幅较小的幼树容易被忽略。由于与周围背景有明显差异，沙丘生境中榆树(图5(a))很容易被提取，并且准确率较高，而塔拉(草原)生境中榆树(图5(b))由于受背景的影响，其提取精度明显不如沙丘生境中的榆树。此外，高大灌木与榆树易混淆，选取的300个榆树类别样本中有69个为误分，其中46个样本实际地类为灌木，其高度约0.5～2.5 m，主要包括榆树幼苗、红柳、锦鸡儿和柴桦等。

4 结论与讨论

目前对沙地榆树疏林的研究多是采用样地调查、实地走访和查阅历史文献等传统方法，并且研究内容多局限于生态方面，运用遥感技术去探讨榆树疏林空间分布的研究很少。本文尝试通过遥感手段使用GF-2数据采用基于GEOBIA技术提取研究区榆树，为沙地榆树信息提取提供有效决策支撑，获得主要结论如下：

1)NDVI阈值可快速掩模大多数非榆树区，准确提取榆树的潜在分布区，减少后续榆树准确提取的工作量，适用于榆树疏林区的粗提取。

2)本文使用GF-2基于GEOBIA方法进一步识别榆树，具有一定的可行性。经验证，其总体精度为88.17%，榆树提取的制图精度达99.14%，基本可以满足榆树疏林区识别的要求。

3)采用单景GF-2影像，结合研究区榆树疏林的分布特性，提出先使用NDVI粗提取榆树分布区，与直接采用GEOBIA方法提取榆树分布区相比，可提高识别精度并加快数据处理速度。

4)考虑到特征选择和规则集构建的不确定性对GEOBIA方法的制约，提出基于SEaTH算法优选特征及自动计算特征阈值，减少人为的主观干扰，取得了较好的效果，使浑善达克沙地榆树疏林信息提取进入了新的阶段，对于该地区榆树疏林管理有一定的帮助。但是浑善达克沙地面积广阔，情况复杂多样，本研究提出的方法需要进一步验证。

5)本文方法也存在一些不足，仍有待进一步改善。例如塔拉生境中榆树由于受背景影响，使用本文方法提取效果比较一般。下一步研究可考虑结合多个时相的遥感数据来提高塔拉生境中榆树的识别准确率，改进本研究方法的不足。