森林资源是可再生的动态资源。为解决调查间隔期(我国国有林场森林资源调查周期一般为5年，集体林区为10年)内森林资源数据存在的“一年清、两年变、三年乱”的问题，提高森林经营和林业可持续发展决策的科学性，需要准确采集森林变化信息，及时更新森林资源数据库。由于森林分布的广域性及森林变化发生地点的分散性和随机性，以往以地面调查为主的森林变化信息采集方法存在工作量大、劳动强度高、效率低、成本高(赵宪文等，2002)、精度低(李春干等，2006)等难以克服的问题，鲜见大规模长期应用的成功案例，这也是我国几十年来一直未能实现森林资源“年度出数”的主要原因，森林变化信息采集成了长期以来森林资源信息更新的瓶颈。鉴于遥感瞬间视场中地表特征随时间发生的变化引起了2个时相图像像元光谱响应的变化(马建文等，2004；周启鸣，2011；陈鑫镖，2012)，通过遥感图像变化检测，可以确定森林变化的空间位置、类型和数量(李世明等，2011)，因此，近10多年来，基于遥感图像变化检测的森林变化信息采集方法逐渐得到广泛应用(Desclée et al., 2006；Potapov et al., 2012；石军南等，2010；魏安世等，2010；王志慧等，2011；王志杰，2015)。

众多学者从不同角度基于不同应用对遥感图像变化检测的方法和理论模型进行了大量研究，提出了很多变化检测方法，如代数法、变换法、分类比较法、高级模型法、GIS集成法、面向对象法、时间序列法、可视化法等(Lu et al., 2004)。但由于传感器不同、检测目标相异、变化特征表现差异大，变化检测十分复杂和困难，没有一种方法能够适合于所有的应用问题，不同方法有时会得出不同甚至相互矛盾的结论，而且自动化程度偏低(周启鸣，2011)。因此，需要根据不同的数据源、不同的森林变化特点进行有针对性的试验。

面向对象变化检测方法通过图像分割将相邻同质的多个像元组成1个对象，并以此为单元进行分析。该方法可充分利用光谱和纹理信息，结合了目视解译的纹理分析和像元法的定量分析的优势(Desclée et al., 2006)，具有受噪声和预处理影响小的特点，并且可结合GIS专题图等其他多源信息，是目前最有潜力的一种方法(周启鸣，2011)。本研究以森林(地)覆盖变化频繁和快速、变化图斑多而小的广西壮族自治区上思县局部区域为研究区，以资源三号和高分一号卫星图像和小班专题图为数据源，采用面向对象方法，深入分析了两时相图像光谱(均值的标准差)差值的统计分布，探讨基于统计检验的面向对象高分辨率卫星遥感图像森林(地)变化的自动检测，为准确、快速、高效采集森林变化信息，及时更新森林资源数据库提供一种有效方法。

研究区位于广西壮族自治区南部的上思县东北部，是一个10.7 km×11.3 km的矩形区域，面积约12 000 hm²，中心地理位置为108°01′50″E, 22°10′17″N。南部为十万大山山脉，北部为四方岭山脉，中部为县城和平缓耕地。最高海拔750 m，最低90 m，地形陡峭，一般坡度20°~35°。地处北回归线以南，属南亚热带季风气候，全年雨热充沛且同季，夏长冬短，冬少严寒，夏少酷热，无霜期短，干湿季明显。年平均日照时数1 896.1 h，年平均气温21.7 ℃，≥10 ℃年活动积温7 700~7 800 ℃，年无霜期339天，年均降雨量1 217 mm。

常见植被为马尾松(Pinus massoniana)林(大部分为天然次生林，少部分为人工林)和桉树[主要为尾叶桉(Eucalyptus urophylla)和巨尾桉(E.grandis× E. urophylla)]人工林。据2013年林地变更调查结果，研究区内林地和非林地的面积比例分别为58.7%和41.3%。马尾松林、桉树林和一般阔叶树天然次生林面积分别占森林面积的52.1%、40.0%和4.3%。

从遥感图像变化检测的角度分析，研究区森林(地)变化具有如下特点：1) 林地变化频繁。桉树人工林为高度集约经营，林木生长迅速，轮伐期为4~6年，因此，森林-无林地(采伐迹地)变化频繁。2) 森林恢复快速。桉树人工造林6个月林分平均高达3.8 m (黄荣林，2006)，2.25~2.67年时平均高为13.1~17.3 m (陈健波等，2011)，3年生年均蓄积生长量为49.8 m³·hm^-2a^-1(张健军等，2012)，最高可达70.3 m³·hm^-2a^-1(项东云等，2008)。森林采伐主要为桉树林和马尾松林，均为皆伐。桉树林采伐后，通过萌芽方式更新，马尾松林采伐后，当年春夏或次年春季通过植苗造林方式更新为桉树林。据调查，2013年上思县森林采伐面积为8 952.3 hm²，至当年底，除24.7%仍保留为采伐迹地外，其余采伐迹地经更新造林后均转化为新造林地。采伐迹地-幼林的转化极为迅速，有研究表明，当成像时间与林分采伐时间间隔大于9个月时，资源一号02C图像和天绘一号图像NDVI差值的标准差的差值很小，难以通过遥感图像进行检测(付强等，2013)。3) 变化图斑多。2013年全县变化图斑多达4 142个，其中造林更新图斑2 255个，采伐(含毁林开垦)图斑1 758个。4) 变化图斑面积小。全县变化图斑平均面积为3.9 hm²，面积小于1.0、2.0和3.0 hm²的变化图斑数量分别占全部变化图斑的46%、65%和75%。

卫星遥感数据为2013年12月7日成像的资源三号(ZY-3) 图像(全色波段空间分辨率为2.1 m，波长为0.50~0.80 μm；多光谱波段空间分辨率为5.8 m，波长范围：蓝波段为0.45~0.52 μm，绿波段为0.52~0.59 μm，红波段为0.63~0.69 μm，近红外波段为0.77~0.89 μm)和2015年1月15日成像的高分一号(GF-1) 图像(全色波段空间分辨率为2 m，波长为0.45~0.90 μm；多光谱波段空间分辨率为8 m，波长范围与ZY-3图像相同)。两景图像分别采用1:10 000正射航空图像(DOM)和分辨率为30 m的DEM进行几何精度校正(正射)，总误差不大于1个像元，并进行分辨率融合，两景图像均重采样至空间分辨率为2 m。以GF-1图像为参考图像对ZY-3图像进行直方图匹配。

小班专题图为2013年森林资源变更调查成果图(1 244个小班，ArcGIS shp格式)。该调查于2014年通过对2013年ZY-3等图像变化检测提取变化图斑边界，经实地逐一核对、修正其边界后对2008年森林资源数据库进行更新而得。小班专题图和遥感图像具有相同的投影系统。

鉴于本研究针对的是森林(地)的变化检测，故检测目标不包含非林地(农用地、建设用地和水体等)。根据研究区森林(地)变化情况，从遥感图像变化检测的角度分析，森林(地)包括3种类型：1) 成林地，在ZY-3和GF-1图像上，林龄较大的林分纹理较为粗糙，林龄较小的林分纹理较为光滑；2) 无林地，包括采伐后未更新或更新后幼林未郁闭成林、完全或部分裸露的林地和已经占用征收的林地，在假彩色图像上表现色彩基本一致，几无纹理，可见林区道路，与周围的成林地区分明显；3) 幼林地，指采伐后造林更新6个月至1年林木尚未完全郁闭的林地，在彩色图像上色彩与成林地相近，纹理不明显，运材林道清晰可见。见图 1。

图 1 成林地、幼林地和采伐迹地在GF-1图像(RGB=红-近红外-蓝)上的表征 Fig.1 An illustration of old, young, and newly-cleared forest area on the GF-1 image (bands 3-red, 4-near infrared, and 1-blue)

将不变作为1个变化类型，则3种林地类型共有9种变化类型：1) 成林地-成林地；2) 成林地-无林地；3) 成林地-幼林地；4) 无林地-无林地；5) 无林地-成林地；6) 无林地-幼林地；7) 幼林地-幼林地；8) 幼林地-成林地；9) 幼林地-无林地。从变化检测的角度看，1)、4)、7) 3种类型为不变类型，其余类型均属变化类型。

在一般面向对象的图像分析中，图像分割均针对一景图像(Flanders et al., 2003；Wulde et al., 2004)，为使分割后的图像对象能准确描述变化图斑的边界，故两时相图像中全部8个波段同时参与分割。图像分割在Defines^®Developer软件平台中进行，为寻求最能准确表现变化图斑边界的图层，进行由大尺度至小尺度的多尺度分割，分割参数见表 1。第1层分割的对象为全部像元，小班专题图参与分割，目的是将2013年小班调查的土地类型属性赋予图像对象，为后续尺度更小的分割作掩膜用途，第2、3层仅对林地范围进行分割。

表 1 图像多尺度分割参数 Tab.1 Parameters for multi-resolution image segmentation

表 1 图像多尺度分割参数 Tab.1 Parameters for multi-resolution image segmentation

采用目视评价方法，确定第3层图斑边界与图像表征最吻合，因此后续图像分析和检测均以该层为基础进行。第3层数据中，共有林地对象(图斑)12 861个，删除面积小于500 m²的碎小斑块后，剩余12 609个，平均每个图斑包含1 405个像元，最少包含125个像元，最多包含34 804个像元。

数据集中共有8个波段的图像数据，其中ZY-3和GF-1图像各有4个波段，分别用ZY-3_B1，…，ZY-3_B4和GF-1_B1，…，GF-1_B4表示。由于每个图斑包含一定数量的像元，因此可以计算第i个图斑在第j个波段的灰度均值(m_ij)：，式中DN_ijk为第i个图斑在第j个波段中第k个像元的灰度值，n为第i个图斑包含的像元数。同理，可以计算各个图斑在各个波段灰度值的标准差(s_ij)。此外，还可计算各个图斑在各个波段的一系列纹理指标、形状指标等。因此，每个对象可以提取大量的特征因子，也即在面向对象的图像分析中，每个图斑可用多个特征组成的多维特征向量进行描述。由于一些对象特征之间具有不同程度的相关性，因此，为确保算法的正常运行、减少运算量，根据应用目的，有针对性地提取相关性小、对分析目标作用明显的对象特征。本研究提取的是每个图斑的灰度均值和标准差，共16个对象特征。将由T1时相(ZY-3图像)、T2时相(GF-1图像)提取的特征向量分别记为X_T1和X_T2，则第i个图斑的特征向量可表示为：

(1)

由于本研究的目的是变化信息提取，故更关注两时相图像的差异情况。两时相图像相应波段具有相同的光谱范围，因此可计算各个图斑4个波段的灰度均值和标准差的差值(可视为差值图像中各图斑的均值和标准差)，两时相特征向量的变化向量(ΔX_i)为：

(2)

显然，ΔX_i反映了第i个图斑对应的地面区域2个年度间光谱响应的变化。一般而言，当林地覆盖变化不剧烈时，如森林不采伐、采伐迹地未更新等，Δm_ij(i, j=1, 2, …, 4) 较小，反之Δm_ij较大，故可用ΔX_i进行森林变化检测。

将各个图斑的X_T1、X_T2和ΔX_i从Defines^® Developer软件中导出。全部图斑在差值图像中的灰度均值用其标准差进行频率分析，结果见图 2。由图 2可知，无论是均值(图 2a)还是标准差(图 2b)，均近似呈正态分布。

图 2 差值图像中4个波段灰度均值和标准差的频率分布 Fig.2 Frequencies of object mean and standard deviation of four bands of the different images

进一步分析发现，在各波段灰度均值和标准差的二维平面空间上，不变图斑和变化图斑的分布具有明显的规律性：不变图斑集中分布于中心区域，而变化图斑分布于外围区域。图 3为变化图斑和不变图斑在近红外波段(第4波段)灰度均值和标准差二维平面空间上的分布。

图 3 2类图斑在差值图像近红外波段灰度均值和标准差二维平面空间上的分布 Fig.3 A scatter plot of changed/unchanged objects on the axes of object mean and standard deviation of NIR spectral band of the different images

由图 2，可以假定差值图像的均值和标准差均呈正态分布，即

(3)

除采伐(皆伐)、造林更新外，研究区森林经营活动主要为1~3年生桉树林的除草、施肥，均在林下进行，无间伐、择伐，不会导致林地覆盖明显变化，同一区域的地物反射特征保持稳定，在两时相图像上的灰度值相差较小。经图像分割取平均值后，不变图斑的对象特征在两时相图像上的差异进一步削弱。因此，可以假定不变图斑的对象特征(均值和标准差)服从正态分布。对于采伐、造林更新活动而导致变化的林分，由于林地覆盖状况发生了剧烈变化，其太阳辐射的反射特征发生了较大变化，在两时相图像上表现为对象特征值相差较大。当将全部图斑在差值图像上的对象特征作频率分布分析时，不变图斑位于0值附近，变化图斑位于分布曲线的前端和后端(图 2)，作灰度均值和标准差二维平面分布分析时，不变图斑位于中心，变化图斑位于不变图斑的外围(图 3)。

由以上分析，对ΔX_i通过x_i=(x_i－x_min)/(x_max－x_min)作中心化处理后，可以假定在差值图像上不变图斑的向量X_i服从均值为m、协方差矩阵为∑的正态分布，则向量X_i的平方和

(4)

为一个新的随机变量，显然，C_i服从自由度为n (n=8) 的χ²分布，即C_i~χ²(n)，则不变图斑的C_i小于χ_1－α²(n)的概率为1－α：

(5)

若C_i>χ_1－α²(n)，则以1－α的置信度确定第i个图斑为变化图斑(Desclée et al., 2006)。

显然，若确知不变图斑的均值m和协方差矩阵∑的真值，则可通过式(4)、(5) 将变化图斑一次性全部检测出来；然而，即使样本足够大，也只能得到不变图斑m和∑的近似估计值和。为此，可通过反复迭代的方式逐步剔除变化图斑，直至数据集中不含变化图斑为止。采用式(4)、(5) 进行变化检测的流程(Desclée et al., 2006)如下：

1) 采用全部数据集，计算差值图像8个向量的和，根据式(4) 计算C_i，根据式(5) 第1次剔除变化图斑；

2) 采用剔除变化图斑后的剩余数据重新计算、和C_i，再一次剔除变化图斑；

3) 重复2)，直至无变化图斑剔除为止。

在ArcGIS软件平台支持下，通过验证样地进行检测精度评估。具体方法如下：1) 采用Hawthorne Analysis Tools for ArcGIS软件(Beyer, 2015)，在试验区范围内随机生成1 000块面积为30 m×30 m的地面验证样地，得到一个精度检验样地文件(ArcGIS shp格式)；2) 将精度检验样地文件与遥感图像、2013年小班分布图叠合，通过目视解译方法逐一确定各个样地所属的变化类型，若样地中包含2个或以上变化类型，以所占面积大的类型为样地变化类型，若2个或以上类型的面积相等，以外延面积大的类型为样地变化类型；3) 将精度检验样地文件分别与检测结果(ArcGIS shp格式)进行空间叠置分析，得到各个样地所属的检测类型；4) 采用混淆矩阵计算相关精度指标。

编写算法的计算机程序，分别以危险率α=0.05、0.02、0.01和0.001(相应的置信度分别为0.95、0.98、0.99和0.999) 进行独立检测(全部运行时间约2 h)，4块检测分别进行了25、23、20和15次迭代，各次检测到的变化图斑数量变化情况见图 4a，α=0.01时各次迭代检测到的变化图斑在差值图像近红外波段均值和标准差二维平面空间的分布见图 4b。α=0.01时的检测结果见图 5。

图 4 4个危险率检测中各次迭代检测到的变化图斑数量变化曲线(a)和α=0.01时各次迭代检测到的变化图斑在差值图像近红外波段灰度均值(d_b4) 和标准差(d_stdb4) 二维平面空间的分布(b) Fig.4 Number changes of change objects detected in every iteration of four detections based on confidence level of 0.95, 0.98, 0.99 and 0.999(a) and their distributions detected on the confidence level of 0.99 on the axes of the difference of DN means and standard deviations of NIR spectral band

图 5 图斑边界与两时相原始图像(a、b，RGB=红-近红外-蓝)叠合效果及α=0.01时的检测结果(c) Fig.5 False color composite subsets (RGB=NIR-Red-Green) of ZY-3 and GF-1 image overlaid by the two period data segmentation results, and the hatched regions on the output map corresponding to changed objects detected with α=0.01

采用混淆矩阵进行检测精度评价，结果见表 2。

表 2 4种置信度(P_1-α=0.95、P_1-α=0.98、P_1-α=0.99和P_1-α=0.999) 检测的混淆矩阵 Tab.2 Error matrix of detection at four confidence levels (P_1-α=0.95, P_1-α=0.98, P_1-α=0.99 and P_1-α=0.999)

表 2 4种置信度(P1-α=0.95、P1-α=0.98、P1-α=0.99和P1-α=0.999) 检测的混淆矩阵 Tab.2 Error matrix of detection at four confidence levels (P1-α=0.95, P1-α=0.98, P1-α=0.99 and P1-α=0.999)

由表 2可知，随着置信度提高，漏检率升高、误检率降低，当置信度为0.99时，总体精度和Kappa系数均最大，说明在这一置信水平上结果最好。

当置信度为0.99时，检测错误主要是将不变的成林地-成林地检测为变化(误检)、将变化的成林地-幼林地检测为不变(漏检)以及将不变的无林地-无林地检测为变化(误检)，此外，无林地-幼林地、幼林地-成林地的漏检率分别达到40%和36%。

漏检错误主要出现在与幼林地相关的变化类型上，包括成林地-幼林地、无林地-幼林地和幼林地-成林地。这3种类型虽属变化类型，但其中一些图斑的地表特征变化不够明显，两时相图像的光谱差异较小，因此，容易造成漏检。无林地-无林地为不变类型，但在3个可靠性水平的检测中，其误检率分别达到75%、65%和45%，表面上不好解释。通过图像观察发现，随着时间的推移，部分林区道路图斑受路旁林木生长和采伐的影响，图像表征发生了较大变化，故误检为变化图斑；一些无林地(荒地)图斑经占用征收后变成了建设用地，地表特征也发生了较大变化，故误检为变化图斑。进一步分析发现，“无林地-无林地”在两时相图像4个波段灰度值差值的标准差远大于其相应的均值，其变异系数高达3.43~42.34，说明一些无林地尽管没有发生类型转换，但其光谱变化很大，因此，在迭代过程中将相当一部分光谱变化较大的“无林地-无林地”检测为变化图斑，造成了误检。

通过前面分析可知，本研究讨论的变化检测方法集图像分割、图像差表达、统计检验和变化图斑识别于一体。

面向对象的图像分析方法综合考虑了遥感图像的光谱特征和空间特征，通过基于光谱异质性标准和空间异质性标准的多尺度分割(Walter, 2004)后，得到了大量光谱高度同质且与相邻区域明显异质的图像对象(图斑)。当两时相图像共8个波段的数据同时参与分割并根据各自的标准差设置其权重时，图像分割过程中考虑了各个局部区域两时相图像的光谱差异，可有效将变化区域与不变化区域分离，得到准确的变化图斑边界(图 5)。此外，在图像分割过程中，小班专题图参与分割，可确保分割后的图斑边界不超出小班边界，并将小班内部光谱异质区域(即变化区域)切割出来，运行检测程序后可将变化区域检测出来，经适当编辑后即可得到用于小班专题图更新的变化图斑图层，为实现以小班为单元的变化信息采集、更新小班空间信息奠定了基础。

在前期试验中发现，不但变化图斑的灰度值差值大于不变图斑，而且其标准差差值也大于不变图斑，因此，在变化检测中，标准差的差值是十分有用的信息，同时利用两时相图像灰度值和标准差的差值信息，有助于提高检测精度。

本研究采用的方法直接对分割图像进行运算，通过迭代方法逐次将变化图斑检测出来，不需任何训练样本和外部信息，具有非监督分类的特性(盛辉等，2004)。

基于统计检验的面向对象高分辨率遥感图像森林变化检测方法具有严密的统计学基础，直接对分割图像进行计算，不需建立训练样本，不依靠任何外来信息，不需人为干预，整个检测过程自动完成，并且检测效果良好，在森林(地)变化检测中具有良好的应用前景。