利用KPCA法检测高分一号影像中的森林覆盖变化

尹凌宇， 覃先林， 孙桂芬， 刘树超， 祖笑锋， 陈小中

(1.中国林业科学研究院资源信息所,北京 100091； 2.四川省林业信息中心,成都 610081)

0 引言

森林是人类陆地生态系统的主体，在净化空气、调节气候、涵养水源、减少风沙危害等方面起着重要作用[1-2]。森林覆盖变化主要包括不同森林类型之间的转化和森林类型内部之间的变化[3]2种形式。利用遥感技术开展的森林覆盖变化检测方法主要分为3类： ①简单代数运算方法,包括影像差值、影像比值、变化矢量分析和影像回归变换方法等，这些方法通过对多时相遥感影像进行简单代数运算，根据反射率值的差异检测出影像变化区域，虽然简单易行，但没有考虑外界干扰因素的影响，因而检测结果的效果较差[4]。②监督分类后比较方法，首先对多时相影像进行监督分类，然后利用分类后的不同地物类别属性进行变化检测。这类方法虽然可以减少大气、传感器差异等外部因素对变化检测的影响提供变化类型信息； 但是在多时相遥感影像分类时，需要选择高质量和足够数量的训练样本，而且因分类不精准会造成累积误差[5]。③基于遥感数据多变量分析变换方法。如主成分分析(principal component analysis,PCA)、卡方变换(Chisquare)和穗帽变换(KT)等方法，其优点是通过数据变换减少变量间的相关性，增强影像信噪比，提高变化区域的提取精度。覃先林等[6]利用MODIS数据，分别采用红光-近红外法、共生矩阵纹理法和基于相似度的变化监测方法，对我国东北林区的森林覆盖变化监测方法进行研究，并进行了精度验证，结果表明，基于相似度的变化监测方法的正确率最高(达到了80% 以上)，而共生矩阵纹理法的漏报率最高。Platt 等[7]利用面向对象的方法对美国科罗拉多州1938—1999年的森林覆盖变化进行研究，获得了森林覆盖变化的空间分布等信息。Hese等[8]利用1989年和2000年2期TM数据，采用面向对象分类法对西伯利亚地区森林覆盖变化进行了分类比较。Collins 等[9]通过对Landsat数据进行线性变化来监测森林死亡率，研究发现PCA和多时相的KT变换比GS正交变换法更好，并且KT变换产生的湿度分量是森林变化监测中最可靠的单一指标。黄维等[10]基于PCA的变化向量分析法，用传统全局阈值法和局部最小错分概率法自动确定阈值，分别提取森林覆盖变化区域。

本文以国产高分一号(GF-1)卫星16 m分辨率多光谱宽幅(WFV)数据为数据源，利用核主成分分析(kernel principal component analysis，KPCA)方法变换构造差值影像，采用最大类间方差法(OTSU)进行自动阈值选取，提取研究区森林覆盖变化信息； 并与传统变化矢量法进行对比验证，以期探寻GF-1 WFV数据在我国森林覆盖变化检测中的应用方法。

1 研究区概况与数据源

1.1 研究区概况

以四川省甘孜州雅江县为研究区。该县位于E100°19′～101°26′，N 29°03′～30°30′之间，面积768 150 hm2； 属于青藏高原亚湿润气候区带，全年平均气温11 ℃(其中1月均温1.4 ℃，7月均温18 ℃)，年降水量650 mm，无霜期189 d，年均日照2 329 h。区内有草地35.35万hm2、林地47.9万hm2、森林20.96万hm2，森林资源覆盖率为27.29%； 主要的树种有冷杉、云杉、桦木、落叶松、柏木和栎树等。

1.2 数据源

为了减少季相性对变化检测结果的影响，本文选取2014年1月7日和2016年1月14日获取的GF-1卫星16 m分辨率WFV多光谱数据作为森林变化监测的遥感数据源。该卫星的主要参数如表1所示。

表1 GF-1卫星WFV参数Tab.1 Parameters for WFV of GF-1 satellite

采用ENVI5.1软件对GF-1 WFV数据进行预处理，主要包括辐射定标、大气校正、正射纠正、相对几何配准及相对辐射归一化等，辐射定标时采用中国资源卫星应用中心公布的定标参数。为了更准确地得到森林变化检测结果，以2014年影像为参考，对2016年影像进行相对几何配准，配准误差控制在0.5个像元以内； 相对辐射归一化则采用迭代加权多元变化检测(iteration re-weight multivariate alteration detection，IR-MAD)法[11]，以2014年影像为参考影像，对2016年影像进行相对辐射归一化处理。

另外，还收集了森林资源调查数据、30m分辨率的ATSSLP DEM数据、行政界线矢量数据以及在该区域进行野外调查获取的相关实测数据。

2 研究方法

2.1 KPCA法

多时相遥感影像各波段间经常存在一定的相关性，这对变化检测提取变化区域增加了困难。目前，一般采用多变量的变换技术消除影像间的相关信息，降低其冗余度，如卡方变换、KT变换和KPCA法。其中，KPCA法不仅能够消除多变量影像中的二维相关信息，而且还能够消除影像间的高维相关，是一种基于核变换非线性推广的主成分分析； 其基本原理是通过引入核函数，将原始数据非线性映射到高维特征空间，然后在新特征空间中使用线性主成分分析提取样本数据的特征[12-14]。假设X=(x1，x2，…，xn)T为影像的像元矢量，N为波段数，S为协方差矩阵，若将X非线性映射到一个高维特征空间Ф，则有

(1)

对S进行特征值分解，求得特征值λi和特征向量μi，即

Sμi=λiμi；

(2)

特征向量μi可表示为Φ(x)的线性组合，即

(3)

(4)

式中a为线性组合系数，a=(a1，a2，…，an)T。

为了减少运算量，定义一个N×N的核函数K，即

Kij=K(xi,xj)=Ф(xi)TФ(xj)，

(5)

Nλia=Ka。

(6)

提取的影像数据x通过KPCA映射后的得到第t个非线性主成分pt，即

(7)

式中：μt为第t个特征向量。

由于变换前各波段之间存在很强的相关性，经过KPCA变换后的高维空间获得的主成分保留了非线性核映射最大的信息量、消除了冗余信息，从而使输出影像的各分量之间具有最小的相关性。

2.2 CVA法

为了更好地对比评价基于KPCA的变化检测方法，本文采用经典的变化矢量分析(change vector analysis，CVA)法[15-16]进行对比实验。CVA法综合利用T1和T2时相遥感影像各波段的光谱信息，通过确定从时相T1到时相T2光谱向量的变化强度和方向来检测变化信息。该方法是对影像直接差值法的扩充，通过对2时相影像各个波段的数据对应进行差值运算，得到每个像元在各个波段的变化增量，由各个波段的变化增量组成变化矢量。其中，变化的强度用变化矢量的欧氏距离表示，变化的内容用变化矢量的方向表示。

假设X=(x1，x2，…，xn)T和Y=(y1，y2，…，yn)T分别为T1时相和T2时相n个波段的遥感影像在同一位置的像元光谱值，xi和yi分别为2时相在第i波段的像元值，则变化矢量ΔC定义为

ΔC=X-Y=[(x1-y1),(x2-y2),…,(xn-yn)]T；

(8)

变化强度定义为

(9)

一般来说，当|ΔC|较大时，表示2景影像之间的差异较大，对应像元的地物类型发生变化的概率较大； 反之则较小。

2.3 差值影像构造

基于CVA的差值影像是由2个时相的影像各个波段的数据进行差值运算得到的，虽然构造差值影像结构简单，运算量较少； 但由于利用影像各个波段分别构造差值影像，影像间存在一定的相关性，对变化检测具有不利的影响，会出现明显“伪变化”和“椒盐现象”。基于KPCA方法构造差值影像时，需要先分别对每期影像进行KPCA变换，以消除影像波段间的相关信息，使变化信息集中在变换后的少数几个波段，然后再构造差值影像。本文提出的基于KPCA法构造差值影像的步骤是: 首先对影像进行KPCA 变换，选取变换后信息量最大的第一主成分量； 然后选取2时相的第一主成分量进行差值运算，得到信息量充足的差值影像。

2.4 OTSU法阈值选取

本文采用最大类间方差法(OTSU)[17-18]自动提取阈值。OTSU法是由日本学者大津展之(Nobuyuki Otsu)于1979年提出的处理影像分割二值化的一种方法(亦称最大类间差法或大津算法)，是一种非参数、自适应的阈值确定方法。该方法基于聚类的思想，把影像的反射率值按反射率级分成2个部分，使得2个部分反射率值之间的差异最大，每个部分内部反射率值之间的差异最小，通过计算其方差来寻找一个合适的划分级别。因此，在区分变化和未变化信息时，可采用OTSU算法自动选取阈值、进行二值化。

对于一景影像，假设对发生变化区域与未发生变化区域的分割阈值为t，其中发生变化区域的像元个数占该影像像元总数目的比例为w0，均值为u0； 未发生变化区域的像元个数占该影像像元总数目的比例为w1，均值为u1； 则整景影像的均值为

u=w0u0+w1u1，

(10)

建立函数

g(t)=w0(u0-u)2+w1(u1-u)2，

(11)

g(t)就是当分割阈值为t时的类间方差表达式。OTSU算法使函数g(t)取得最大值时，把所求出的t称为最佳阈值。

从L个反射率级遍历t，使得t为某个值的时候，发生变化区域和未发生变化区域的方差最大，则这个t值便是所要求得的阈值。

3 结果与分析

3.1 相对辐射归一化结果分析

采用IR-MAD法提取出2时相未变化像元后，就可以利用这些像元点建立线性回归方程，然后对目标影像进行相对辐射归一化处理。以2014年GF-1影像为参考，对2016年GF-1影像进行相对辐射归一化处理，建立的回归方程如图1(反射率被拉伸至0～10000)所示。为了定量检验线性回归效果，分别计算出每个波段的决定系数R2。

(a) B1 (b) B

(c) B3 (d) B4

图1各波段未变化像元散点图及回归方程

(目标影像为2016年GF-1影像，参考影像为2014年GF-1影像)

Fig.1Scatterplotsandregressionequationsofnochangedpixelsineachband

图1显示，利用 IR-MAD法提取的2期GF-1影像的未变化像元点，各波段的R2都超过0.99，拟合效果非常好。为了定量检验相对辐射归一化效果，对参考影像和相对辐射归一化前后的目标影像进行统计，提取出的未变化像元各波段反射率的均值如表2所示。

表2 相对辐射归一化前后目标影像和参考影像各波段反射率均值Tab.2 Average reflectivity of different bands of target image before and after radiation normalization and reference image

从表2可以看出，对2期GF-1影像相对辐射归一化处理后，目标影像和参考影像各波段的反射率均值相近，差异明显减小。

3.2 基于KPCA的变换结果分析

本文选择高斯核函数对2期影像分别进行KPCA变换，分别得到4个主成分(PC)变量。2期影像的KPCA变换的定量统计结果见表3，统计信息包括特征值、方差贡献率和累计方差贡献率。其中，方差贡献率是第k个KPCA主成分的方差在全部方差中所占的比重，反映了该主成分包含信息量的大小； 累积方差贡献率表示前k个KPCA主成分的共有信息量，为k个主成分的方差和在全部方差中所占的比重。从表3可以看出，2期影像的KPCA第一主成分的方差贡献率均超过了91%，说明该变量包含了最多的变化信息； 其余各变量的特征值依次减小，说明后面几个KPCA变量中包含的有效变化信息越来越少，所包含的噪声可能会很多，降低了变化信息提取精度。因此，本文利用2期影像第一主成分变量构造差值影像，进行变化检测。

表3 GF-1影像KPCA主成分方差分布Tab.3 Variance distribution of principal components of GF-1 image KPCA

3.3 OTSU阈值选取结果

采用OTSU法对已得到的差值影像提取最佳阈值。选取雅江县东部的1块GF-1局部影像(图2)进行放大和对比分析。

(a) 2014-01-07 (b) 2016-01-14

图2雅江县东部局部影像

Fig.2LocalimagesofeastpartinYajiangCounty

对该区域采用KPCA和CVA这2种方法所得到差值影像，使用阈值确定方法得到2个变化区域的二值影像(图3)。

(a) KPCA法 (b) CVA法

图3变化区域二值影像

Fig.3Binaryimagesofchangedareas

从图3可以看出，对2种方法构造的差值影像，利用阈值法都能较准确地提取出主要的变化区域。

3.4变化信息提取精度检验

为了验证影像变化信息提取的准确性，对所采用的2种变化检测方法的变化信息提取结果分别进行了精度评价。

利用对2014年1月29日发生在雅江县八角楼乡的森林火灾形成的火烧迹地的GPS采点调查结果，结合通过目视解译，在GoogleEarth 平台上分别从2013 年11月26日和2016年3月7日获取的GF-1影像中选取变化与未变化样本点(图4)，其中变化样本点131个，未变化样点289个； 计算2种检测方法变化信息提取的错分误差、漏分误差、总体精度等指标，得到基于KPCA和CVA差值法的变化检测精度(表4)。

(a) KPCA法 (b) CVA法

图4变化区域提取结果

Fig.4Resultextractedfromchangedareas

表4 2种变化检测方法精度Tab.4 Accuracy of two methods for change detection

从表4可以看出，本文采用的2种变化检测算法的总体精度都超过了80%，其中CVA法的提取精度为83.33%，KPCA法比其精度提高了5.94%，达到了89.27%。对于KPCA法，未变化区的用户精度最高，达93.88%； 变化区用户精度最低，仅为80.28%，说明变化区误分最为严重； 未变化区生产者精度为90.32%，变化区生产者精度为87.02%，说明变化区部分存在漏分现象。作为对比的CVA法的类别提取精度均低于KPCA法； 对变化区而言，生产者精度远低于KPCA法，仅为70.23%，说明变化区漏分现象很严重。这是因为CVA算法混合了很多噪声信息，掩盖了部分变化信息，因而出现“伪变化”和“椒盐现象”； 而KPCA算法通过数据变换后，减少了变量间的相关性，增强了影像的信噪比，提高了变化区域的提取精度。

4 结论

本文探讨了基于不同年份的2期GF-1WFV多光谱遥感数据进行森林类型覆盖变化信息检测的方法，得出如下结论：

1)采用迭代加权多元变化检测(IR-MAD)法的相对辐射归一化法，使作为参考影像和目标影像的2期GF-1影像各波段间的的均值差异明显减小，克服了2期影像成像质量的差异，各波段的R2都超过0.99，有效消除了不同时相间的辐射差异。

2)采用KPCA方法对2期影像的有效信息进行提取，得到第一主成分的信息量均达到了91%，使差值变化信息更集中，有利于提高变化监测精度。

3)经精度验证，采用KPCA方法和CVA方法的总体检测精度达到80%以上，都能检测出森林火灾引起的森林覆盖变化； 但由于火灾燃烧程度不同，造成火烧迹地影像的灰度值不同，因而造成火烧迹地、尤其是火烧迹地的边界都存在被误判或漏判的现象。

本文研究的KPCA方法能提取GF-1 WFV影像中的有效变化信息，用变换后的第一主分量构造差值影像，虽然可以有效抑制噪声，但在一定程度上舍弃了部分影像信息，可能会遗漏部分森林变化区域。下一步研究将结合遥感影像中的纹理、形状、颜色等特征信息，综合利用多种特征进行变化检测，以提高森林覆盖变化检测的精度。

[1] Fazakas Z,Nilsson M,Olsson H.Regional forest biomass and wood volume estimation using satellite data and ancillary data[J].Agricultural and Forest Meteorology,1999,98-99:417-425.

[2] 李俊明,邢艳秋,杨 超,等.基于环境与灾害监测预报小卫星的树种识别[J].东北林业大学学报,2013,41(11):41-45,50.

Li J M,Xing Y Q,Yang C,et al.Forest type identification based on hyperspectral remote sensing image of environment and disaster monito-ring satellite[J].Journal of Northeast Forestry University,2013,41(11):41-45,50.

[3] 吴雪琼,覃先林,周汝良,等.森林覆盖变化遥感监测方法研究进展[J].林业资源管理,2010(4):82-87.

Wu X Q,Qin X L,Zhou R L,et al.Progress of study on forest cover change detection by using remote sensing technique[J].Forest Resources Management,2010(4):82-87.

[4] 钟家强,王润生.基于独立成分分析的多时相遥感图像变化检测[J].电子与信息学报,2006,28(6):994-998.

Zhong J Q,Wang R S.Multitemporal remote sensing images change detection based on ICA[J].Journal of Electronics & Information Technology,2006,28(6):994-998.

[5] 李世明,王志慧,韩学文,等.森林资源变化遥感监测技术研究进展[J].北京林业大学学报,2011,33(3):132-138.

Li S M,Wang Z H,Han X W,et al.Overview of forest resources change detection methods using remote sensing techniques[J].Journal of Beijing Forestry University,2011,33(3):132-138.

[6] 覃先林,陈尔学,李增元,等.基于MODIS数据的森林覆盖变化监测方法研究[J].遥感技术与应用,2006,21(3):178-183.

Qin X L,Chen E X,LI Z Y,et al.Forest cover change monitoring using MODIS data[J].Remote Sensing Technology and Application,2006,21(3):178-183.

[7] Platt R V,Schoennagel T.An object-oriented approach to assessing changes in tree cover in the Colorado front range 1938―1999[J].Forest Ecology and Management,2009,258(7):1342-1349.

[8] Hese S,Schmullius C.Forest cover change in Siberia-results from the Siberia-II project[EB/OL].(2011-01-11)[2016-08-01].http://www.Ecognition.com/document/forest-cover-changesiberia-results-siberia-ii-project.

[9] Collins J B,Woodcock C E.An assessment of several linear change detection techniques for mapping forest mortality using multitemporal Landsat TM data[J].Remote Sensing of Environment,1996,56(1):66-77.

[10] 黄 维,黄进良,王立辉,等.基于PCA的变化向量分析法遥感影像变化检测[J].国土资源遥感,2016,28(1):22-27.doi:10.6046/gtzyyg.2016.01.04.

Huang W,Huang J L,Wang L H,et al.Remote sensing image change detection based on change vector analysis of PCA component[J].Remote Sensing for Land and Resources,2016,28(1):22-27.doi:10.6046/gtzyyg.2016.01.04.

[11] 刚慧龙.基于IR-MAD算法的多元遥感影像变化检测方法研究[D].青岛:山东科技大学,2014.

Gang H L.The Study of Change Detection Method for Multivariate Remote Sensing Data based on the IR-MAD Algorithm[D].Qingdao:Shandong University of Science and Technology,2014.

[12] 沈照庆,陶建斌.基于模糊核主成分分析的高光谱遥感影像特征提取研究[J].国土资源遥感,2009,21(3):41-44,99.doi:10.6046/gtzyyg.2009.03.08.

Shen Z Q,Tao J B.Hyperspectral remote sensing image feature extraction based on fuzzy kernel principal component analysis[J].Remote Sensing for Land and Resources,2009,21(3):41-44,99.doi:10.6046/gtzyyg.2009.03.08.

[13] 杨 帆,马 洁.基于KPCA数据重构的参数优化方法[J].北京信息科技大学学报,2014,29(6):46-50.

Yang F,Ma J.Parameter optimization method based on KPCA data reconstruction[J].Journal of Beijing Information Science and Technology University,2014,29(6):46-50.

[14] 王国胜.核函数的性质及其构造方法[J].计算机科学,2006,33(6):172-174,178.

Wang G S.Properties and construction methods of kernel in support vector machine[J].Computer Science,2006,33(6):172-174,178.

[15] Malila W A.Change vector analysis:An approach for detecting forest changes with Landsat[C]//Proceedings of Sixth International Symposium on Machine Processing of Remotely Sensed Data.West Lafayette,Indiana:Purdue University,1980.

[16] 林克正,班守峰,张玉双.一种基于变化向量分析的变化检测方法[J].哈尔滨理工大学学报,2008,13(4):47-49.

Lin K Z,Ban S F,Zhang Y S.A change detection based on change vector analysis[J].Journal of Harbin University of Science and Technology,2008,13(4):47-49.

[17] 魏伟一.基于小波域灰度拉伸的Ostu图像分割[J].西北师范大学学报(自然科学版),2009,45(6):46-48.

Wei W Y.New Ostu image segmentation based on intensity stretching on DWT fields[J].Journal of Northwest Normal University(Natural Science),2009,45(6):46-48.

[18] 蔡燕伟,崔红霞.基于彩色空间的二维OSTU的图像分割算法[J].电子设计工程,2014,22(5):131-133.

Cai Y W,Cui H X.Image segmentation algorithm based on two-dimensional OSTU of color space[J].Electronic Design Engineering,2014,22(5):131-133.