当前，世界上很多国家都采用遥感技术从宏观尺度对森林资源动态变化进行监测，我国森林资源调查部门也将多时相多光谱遥感影像变化检测技术用于全国采伐限额、全国林地征占用检查业务，特别是近几年来，随着高分一号、高分二号、ZY-3等国产光学遥感卫星影像质量大幅度提升，对同一地区重复覆盖频率大大提高，有力支撑了国家森林覆盖动态变化监管相关业务。高分三号是2016年我国成功发射的首颗民用多模式合成孔径雷达(SAR)卫星，具有单极化、双极化和全极化等多种数据获取模式。为了尽快将高分三号多极化SAR卫星影像用于国家森林资源宏观监管业务，并为未来国内外L-波段多极化SAR的林业应用做好技术储备，研发自动化、高精度的双时相星载多极化SAR森林覆盖变化检测方法具有重要价值。

基于单极化SAR影像的自动变化检测发展较早，经典的图像二值化算法KI(kittler and illingworth)主要针对高斯分布下的光学影像，基于SAR信号服从的概率密度分布函数的阈值自动确定算法GKIT用于SAR变化检测的阈值分割具有更强的科学性(Moser et al., 2006)。马尔科夫随机场在不同时相间的极化比值影像具有较强的抗噪能力，Fisher线性变换可综合不同极化影像所反映的地物信息，依据差异影像可能服从的概率密度分布函数，采用最大期望算法(EM)和马尔科夫随机场(MRF)实现变化、未变化分类(Moser et al., 2007; 2009)。EM-MRF实际上是一种两类监督分类方法，以GKIT自动产生的变化检测结果为训练数据，可实现整个变化检测过程的自动化。国内学者对该方法进行了改进，根据原始SAR影像的分布特征推导出差异影像的概率密度函数，利用广义Gamma分布拟合出精确的SAR比值分布模型，表现出了良好的变化检测性能(王虎清，2009; 沈蕾，2010; 高丛珊等，2010; Xiong et al., 2012)。通过对多极化SAR乘性模型的改进，利用包含纹理信息和斑噪协方差矩阵的变化测度来衡量不同时相极化SAR影像的变化程度，得到不同类别间差异度更强的差异影像，并利用GKIT提取变化区域，检测精度相比传统的Wishart距离有所提高(刘萌，2013)。王纬华等(2016)以高斯分布为基础，利用马尔科夫随机场去除孤立像素，当满足影像总能量最低时停止迭代，实现了对建筑物和道路的变化检测。以往研究旨在检测地表覆盖是否发生变化，主要针对水灾、城市扩张和农业生产活动等引起的一般地表变化，并非针对森林覆盖变化的检测; 且一些地表覆盖变化检测方法仍采用高斯分布函数拟合SAR影像的概率分布(孙晓霞，2011; 张斌，2013)，并不适用以极化比值为差异影像的森林覆盖变化检测。

2010年，国内外开始了结合森林资源监测业务化应用需求的变化检测技术研究，且采用的基本上都是多时相ALOS PALSAR数据。Santoro等(2012)以对数高斯分布为比值差异影像概率密度分布函数，利用GKIT和EM-MRF进行瑞典森林皆伐(clear-cut)区检测，探讨ALOS PALSAR影像不同像元大小对检测精度的影响，结果发现以20 m空间分辨率的像元大小为基本检测单元时可获得较好的精度。董锡超(2014)采用ALOS PALSAR影像强度、纹理、ScanSAR时间序列影像进行变化检测研究，结果发现融合3种特征并采用均值滤波预处理方法可以实现较好的毁林检测。Pantze等(2014)通过辐射归一化处理不同时相的SAR影像，采用对数高斯分布模型下的广义KI算法，利用MRF融合不同通道SAR检测结果，发现基于自适应滤波的双极化SAR变化检测可达到较低的总体错误率。Mermoz等(2016)运用已有土地覆盖利用产品掩膜掉耕地、水体、城镇等非森林覆盖区域，采用EM对ALOS PALSAR的单极化数据进行森林扰动检测研究，获得2007—2010年越南、柬埔寨、老挝的年度森林损坏和再生长影像，验证HV极化的比值影像是较好的森林扰动变化特征，具有可推广性。以往研究的检测目标包括森林皆伐、森林扰动和毁林等人为或天然因素，其所导致的后果都可以归并为森林覆盖的变化。

综上可以看出，非针对森林覆盖的多极化SAR地表变化检查方法是以SAR信号统计理论为基础，可充分利用影像上下文信息的高精度、自动化方法; 而针对森林覆盖的变化检测技术研发重点在于对国际上最新发展的地表覆盖自动化检测算法进行改进和集成应用，大都采用已有土地覆盖/利用图对多极化SAR地表覆盖变化检测结果进行掩膜处理，以消除非森林覆盖区易出现大量虚警的问题。国外在这方面已开展了很好的应用探索，但国内类似研究工作尚十分缺乏，不利于国产卫星SAR在国家森林资源监测中的业务化应用。虽然与高分三号最近似的国外卫星是Radarsat-2，但由于该卫星是商业卫星，历史存档数据很少，很难获取到适合我国林区森林覆盖变化检测的多时相数据。ALOS PALSAR是一颗L-波段多极化、全极化SAR卫星，2007—2010年，每年实现1次全球覆盖，并将定量化处理得到的25 m双极化数据产品和森林-非森林制图产品实现全球免费共享(Shimada et al., 2014)。鉴于此，本研究采用ALOS PALSAR数据，发展一种半自动化的双极化SAR森林覆盖变化检测方法，以期为我国高分三号及未来L-波段多极化SAR卫星的国家森林资源监测业务化应用提供技术支撑。

1 研究区概况与数据 1.1 研究区概况

研究区位于黑龙江省逊克县西北部，地理位置127.83°—128.02°E，49.37°—49.49°N，面积142 km²，海拔170~290 m，年均气温0.5 ℃，年降水量650 mm，植被生长期较短。图 1a灰色填充多边形为研究区所在逊克县边界，图 1b给出了本研究采用的整景PALSAR影像对逊克县的覆盖情况。研究区主要地表土地利用覆盖类型包括林地、耕地、城乡居民用地和水域等。主要树种有蒙古栎(Quercus mongolica)、白桦(Betula platyphylla)、黑桦(Betula dahurica)、山杨(Populus davidiana)、榆(Ulmus pumila)和兴安落叶松(Larix gmelinii)，主要森林类型包括蒙古栎林、杨桦林(包括白桦林、黑桦林)及阔叶混交林等。研究区自2000年实施天然林保护工程以来，很少发生大面积的人为毁林事件，森林覆盖变化主要由森林火灾、合法采伐和退耕还林等引起，变化面积一般较小，无疑增加了森林覆盖变化遥感自动或半自动检测的难度。本研究森林覆盖是指郁闭度大于或等于0.20的乔木林地和灌木林地。

1.2 双极化SAR遥感数据

ALOS卫星是日本对地观测卫星，其有效载荷之一为L-波段相控阵合成孔径雷达(简称PALSAR)，具有高分辨率、扫描和极化3种观测模式。本研究采用其双极化(HH、HV)FBD模式Level 1.1级单视复数(SLC)数据产品，距离向像元大小为9.37 m，方位向像元大小为3.17 m，影像中心入射角为38.7°。2期影像的成像时间分别为2007年6月22日和2009年6月27日，覆盖范围东西向70 km、南北向58 km。为了更好地对检测算法进行定量评价，选取图 1b方框所示范围的PALSAR数据开展研究，对应的前一期SAR影像如图 1c所示。

1.3 变化检测精度评价参考图

获取覆盖研究区的2期Landsat-5 TM影像，成像时间分别为2007年6月11日和2009年5月14日，与2期ALOS PALSAR成像时间尽量保持一致。通过对2期Landsat-5 TM影像的目视解译，人工勾绘森林覆盖变化面积相对较大的区域。对一些森林覆盖面积变化较小、在Landsat-5 TM影像上变化特征不太明显的区域，采用Google Earth存档的高空间分辨率遥感影像进行仔细确认并解译勾绘。目前Google Earth已成为中低分辨率遥感土地覆盖、变化检测结果验证的重要参考数据来源(Mermoz et al., 2016)。

采用该方法解译出森林覆盖变化区域13个，包含641个变化像元; 其他区域为森林覆盖未变化区域。图 2给出了13个森林覆盖变化区域的中心位置，背景为2期Landsat-5 TM影像(RGB组合方式为近红外、红色、绿色波段)。图 2中2期影像及变化区域中心位置已采用SAR影像预处理阶段建立的SAR影像地形校正地理编码(GTC)模型转换到多视强度影像(MLI)坐标空间，图 2右侧图为解译得到的森林覆盖变化参考图，将用于评价本研究所提出方法的有效性。

2 研究方法

森林覆盖变化检测技术流程如图 3所示。首先需要对PALSAR数据进行定量化预处理; 后续的Fisher特征变换和EM-MRF分类都属于两类别(变化、未变化)有监督分类算法，因此需要先生成一个初始化的森林覆盖变化检测图，仅含变化、未变化2类，为Fisher特征变换和EM-MRF分类提供训练数据。在农田、草地、灌木等非森林覆盖类型分布区，尽管2期SAR影像成像时间在同一月份，但不同年份植被物候期也可能会有所变化，且地表土壤含水量和粗糙度等发生变化的可能性也较大，再加上SAR影像固有的斑点噪声现象，若直接采用Moser等(2009)提出的完全自动化地表覆盖变化检测方法，则会出现大量变化检测虚警。为解决该问题，本研究分别对2期SAR多视强度影像(MLI₁、MLI₂)进行森林-非森林分类，得到初始森林覆盖变化检测图，进而用于Fisher特征变换和EM-MRF分类。最后，以森林覆盖变化参考图为地面实况数据对森林覆盖变化检测精度进行定量评价。

2.1 定量化预处理

首先采用干涉SAR主、辅影像自动配准算法，建立SLC₂与SLC₁的空间匹配模型，将SLC₂重采样到SLC₁影像坐标空间; 其次按照PALSAR辐射定标公式分别对SLC₁、SLC₂进行辐射定标处理，但不进行分贝化处理，以便直接用于后续多视化、比值处理; 然后进行多视化处理，结合后续变化检测试验，确定按距离向1视、方位向5视进行多视化处理可得到最佳变化检测结果，处理结果为MLI₁、MLI₂，各包含HH、HV 2个极化数据，数据类型为32位浮点; 再利用ASTER DEM基于SAR影像模拟法，通过MLI₁的HH极化影像建立地形校正地理编码模型(M_GTC)，将森林覆盖变化参考图变换到MLI₁所在的影像坐标空间。这样就实现了2期SAR数据和森林覆盖变化参考图在MLI坐标空间的严格配准。利用自编写的脚本文件驱动GAMMA软件的相关处理模块，实现以上整个定量化处理过程的自动化。

2.2 森林-非森林分类

只要对MLI₁、MLI₂分别进行森林-非森林分类，获得2期森林-非森林分布图，就可以通过简单的空间叠加分析计算得到森林覆盖变化检测图，这是一种经典的变化检测方法，即所谓的“2期分别分类森林覆盖变化检测法”。为了提高该处理步骤的自动化程度，本研究采用阈值分割法进行森林-非森林分类。首先在MLI₁影像上分别选择若干森林、非森林地表覆盖类型的感兴趣区域(AOI)，通过2个类别HV极化像元值的统计直方图分析，确定森林-非森林分割阈值，采用阈值分割法将每个像元分为森林、非森林；然后采用同样的分析和处理方法得到MLI₂的森林-非森林分类结果。将2个森林-非森林分布图叠加分析得到初始森林覆盖变化图。

2.3 Fisher特征变换

为综合利用HH、HV极化所包含的地类变化特征信息，本研究将HH比值、HV比值的对数进行Fisher特征变换，计算得到一个新的综合变量用于EM-MRF分类。

设u₁=HH₁/HH₂、u₂=HV₁/HV₂分别代表HH极化比值和HV极化比值。对u₁、u₂做如下变换：

$ v=u_{1}^{{{\alpha }_{1}}}\times u_{2}^{{{\alpha }_{2}}}。$

(1)

对上式两边取对数得到线性方程：

$ \ln v={{\alpha }_{1}}\ln {{u}_{1}}+{{\alpha }_{2}}\ln {{u}_{2}}。$

(2)

以初始森林覆盖变化图为训练数据，分别计算变化、未变化类别像元值对数的类内散布矩阵S_i[式(3)]、总类内散布矩阵S_w[式(4)]和总类间散布矩阵S_b[式(5)]：

$ {{S}_{\text{i}}}=\sum\limits_{x\in {{D}_{i}}}{\left(x-{{m}_{i}} \right)}{{\left(x-{{m}_{i}} \right)}^{\text{T}}}; $

(3)

$ S_{\mathrm{w}}=S_{1}+S_{2}; $

(4)

$ {{S}_{\text{b}}}=\left({{m}_{1}}-{{m}_{2}} \right)\left({{m}_{1}}-{{m}_{2}} \right)^{\text{T}}。$

(5)

式中：i=1，2，代表 2个类别；D_i代表 2个不同类别的像元；x=[lnu₁, lnu₂]^T；m_i为x的平均值，上标T表示向量转置运算。

Fisher特征变换就是要找到使得准则函数J(α)[式(6)]最大化时的直线变换方程[式(2)]的参数向量(α)。α的计算公式如式(7)所示：

$ J(\alpha)=\frac{{{\alpha }^{\text{T}}}{{S}_{\text{b}}}\alpha }{{{\alpha }^{\text{T}}}{{S}_{\text{w}}}\alpha }; $

(6)

$ \alpha=S_{w}^{-1}\left(m_{1}-m_{2}\right)。$

(7)

2.4 EM-MRF分类

EM算法实现分为E步和M步。E步根据初始森林覆盖变化图推导出的数据特征参数计算对数似然函数的期望，M步用于选择期望最大值时的参数，从而更新类别信息和特征参数，代入到E步中，如此反复，直到算法收敛。

由于EM算法需要根据差异影像中不同类别像元值的概率密度函数求得条件概率，用于判断像元所属类别，因此选择并高精度拟合变化、未变化类别像元值的条件概率密度函数是EM算法的关键，拟合原则为保证真实概率分布(TPD)和拟合概率函数(FPD)的分布趋势尽可能接近。本研究利用森林覆盖变化参考图，在Fisher特征变换后的差异影像中分别对变化、未变化类别的像元值进行对数高斯概率密度函数[式(8)]拟合：

$ P\left({{v}_{n}}|{{\delta }_{i}} \right)=\frac{\exp \left[ -\frac{{{\left(\ln {{v}_{n}}-{{k}_{1i}} \right)}^{2}}}{2k_{2i}^{2}} \right]}{{{k}_{2i}}{{v}_{n}}\sqrt{2\pi }}, n\in N, i=1, {{2}_{{} }}。$

(8)

式中：N为像元数；v_n为像元灰度值；δ_i为分布函数的参数集合，包括对数均值k_1i和对数标准差k_2i。

图 4展示了2个类别的TPD和FPD拟合效果，说明变化、未变化类别像元值的概率密度函数可用对数高斯函数进行很好拟合。

考虑到马尔科夫随机场和Gibbs随机场的等价性，可以通过最小化表示影像的能量函数获得全局的最优类别信息(张斌，2013)。在实际应用中，通常会按照与中心像元的距离大小将周围像元划分成不同阶数的邻域系统，根据不同的局部概率定义方式，常见MRF模型有Ising模型、Potts模型、MLL模型和多尺度MRF模型等。高阶的邻域系统计算压力较大，不具普及性，因此应权衡邻域系统的阶数和计算复杂度，在获得较高精度的同时具有较快的计算速度。本研究采用二阶邻域系统上的Potts模型，在保证较好计算精度的同时减轻计算压力。Potts模型在二元势函数的定义下对影像进行二值化处理，设x_a是中心像元类别属性，x_b是邻域系统Z内某个像元类别属性，将势函数的参数定义为β，本研究的势函数定义见式(9)：

$ C\left( {{x}_{a}},{{x}_{b}} \right)=\left\{ \begin{align} & 0\left( {{x}_{a}}={{x}_{b}} \right); \\ & \beta \left( {{x}_{a}}\ne {{x}_{b}} \right)。\\ \end{align} \right. $

(9)

推导出能量函数[式(10)]及条件概率[式(11)]：

$ U\left(x_{a}\right)=\sum C\left(x_{a}, x_{b}\right)=\beta f_{a}\left(x_{a}\right); $

(10)

$ P\left({{x}_{a}}|{{x}_{{{Z}_{a}}}} \right)=\frac{\exp \left[ -U\left({{x}_{a}} \right) \right]}{\sum\limits_{{{x}_{a}}\in {{D}_{i}}}{\exp }\left[ -U\left({{x}_{a}} \right) \right]}=\frac{\exp \left[ -\beta {{f}_{a}}\left({{x}_{a}} \right) \right]}{\sum\limits_{{{x}_{a}}\in {{D}_{i}}}{\exp }\left[ -\beta {{f}_{a}}\left({{x}_{a}} \right) \right]}。$

(11)

式中：f_a为Z中与x_a不同相的元素数量。

通过以上势函数的描述，本研究将势函数的参数β定义为1，将Potts模型与式(8)中的条件概率相结合，得出MRF模型新的能量函数：

$ U\left(D_{i} | v_{n}, c_{n}, \theta\right)=\frac{\left(\ln v_{n}-k_{1 i}\right)^{2}-k_{2} \ln k_{2 i}}{2 k_{2 i}}-\lambda f_{i n}。$

(12)

式中：c_n表示公式中的参数集合，包括k_1i、k_2i和λ；θ表示Potts模型的信息；系数λ用于平衡二者之间的权重关系。

将式(12)计算的能量函数代入式(11)，得出新的类别条件概率公式：

$ P\left({{L}_{n}}={{D}_{i}}|{{v}_{n}}, {{c}_{n}}, \theta \right)=\frac{\exp \left[ -U\left({{D}_{i}}|{{v}_{n}}, {{c}_{n}}, \theta \right) \right]}{\sum\limits_{i=0}^{1}{\exp }\left[ -U\left({{D}_{i}}|{{v}_{n}}, {{c}_{n}}, \theta \right) \right]}。$

(13)

式中：L为每个像元的类别标签。

根据以上基本理论和算法描述，确立本研究的具体算法步骤如下：

1) 利用初始森林覆盖变化检测图的类别信息，采用最大似然法估计变化、未变化类别的概率密度分布参数，从而拟合差异影像的概率密度分布。

2) 将最大期望条件概率与影像空间上下文信息相结合，计算像元的类别条件概率，根据最大后验概率准则，更新像元的类别标签。

3) 采用式(14)~(16)更新变化、未变化类别的概率分布参数：

$ k_{1 i}^{t+1}=\frac{\sum\limits_{n=1}^{N} w_{i n}^{t} \ln v_{n}}{\sum\limits_{n=1}^{N} w_{i n}^{t}}, i=1, 2; $

(14)

$ k_{2 i}^{t+1}=\frac{\sum\limits_{n=1}^{N} w_{i n}^{t}\left(\ln v_{n}-k_{1 i}^{t+1}\right)^{2}}{\sum\limits_{n=1}^{N} w_{i n}^{t}}, i=1, 2 ; $

(15)

$ \lambda^{t+1}=\operatorname{argmax} \sum\limits_{n=1}^{N}\left[\lambda \sum\limits_{i=1}^{2} w_{i n}^{t} f_{i n}^{t}-\ln \sum\limits_{i=1}^{2} \exp \left(\lambda f_{i n}^{t}\right)\right], \lambda>0。$

(16)

式中：t为迭代次数；w_in^t=P(L_n=D_i|v_n, c_n^t, θ^t)；λ初始值设置为0.3，在实际应用中可采用牛顿-拉夫逊方法计算λ的近似值。

4) 返回第2步迭代运算，直到相邻2次迭代条件概率的差值(ΔP)小于设定阈值，或迭代次数达到设定次数(N_I)为止。

2.5 精度检验方法及对比试验方案

根据以上流程得到森林覆盖变化图，将其与参考图进行比较，采用如下指标定量评价变化检测效果：1)正确率(TDA)，检测到的正确变化像元和正确未变化像元总和与全部像元的百分比; 2)虚警率(FA)，被错分到变化类别的像元个数与参考影像中未变化像元个数的百分比; 3)漏警率(MA)，没有检测到的变化像元个数与参考影像中变化区域像元个数的百分比。

将以上描述方法与如下2种方法进行比较：

1) 2期分别分类森林覆盖变化检测法(CBFC)：采用该方法获得2期森林-非森林分布图，通过对2期分布图的空间叠加分析得到森林覆盖变化图(FCCM)。

2) CBFC与EM-MRF变化检测结果直接掩膜法(CBFC-EM-MRF)：采用HV比值影像通过GKIT自动阈值分割获得土地覆盖变化图作为EM-MRF分类的输入，得到土地覆盖变化图(LCCM); 然后利用产生的FCCM从LCCM中掩膜掉森林覆盖未变化区域，得到森林覆盖变化检测图。

2.6 森林覆盖变化检测图的地理编码处理

以上处理得到的森林覆盖变化检测图仍然位于MLI坐标空间，采用已建立的M_RTC模型将其地理编码到大地坐标空间，得到最终的森林覆盖变化检测结果(该步骤未在图 3上画出)。

3 结果与分析 3.1 差异影像特征比较分析

图 5展示了HH、HV极化幅度比值差异影像和经Fisher特征变换后得到的差异影像，与图 2中的森林覆盖变化参考图进行对比分析，表 1为3幅差异影像中森林覆盖变化、未变化的类间方差、类间距、类内信息熵、总体信息熵和等效视数。

虽然HH幅度比值影像的总体信息熵和类内信息熵最低，但是类间方差和类间距远远小于HV幅度比值影像和Fisher特征变换后影像，因此HH比值影像难以实现较高精度的森林覆盖变化检测。与HV幅度比值影像相比，Fisher特征变换后的影像类间方差和类间距变大，类内信息熵和总体信息熵均变小，且差异影像的等效视数由0.31(HH)、0.56(HV)升高到1.66，影像噪声得到了一定程度抑制。以上数据说明以Fisher特征变换后差异影像作为EM-MRF的输入更有利于森林覆盖的变化检测。

3.2 EM-MRF分类算法收敛性分析

利用EM-MRF模型对Fisher特征变换后的差异影像进行迭代分类，设置最大迭代次数N_I=50，相邻2次迭代条件概率差值ΔP < 0.001，满足其中一个条件即停止迭代，输出变化检测结果。图 6所示为相邻2次迭代之间的条件概率差值随迭代次数的变化趋势，随着迭代次数增加，差值逐渐变小，当迭代次数达到10次左右时，差值变化趋于稳定，显示本研究提出的变化检测方法具有良好的收敛性。

3.3 森林覆盖变化检测效果定性分析

图 7所示为3种检测方法的森林覆盖变化检测结果，白色表示森林覆盖变化区域，灰色表示森林覆盖未变化区域。可以看出，CBFC方法存在严重的检测虚警，主要原因是为了提高自动化程度，CBFC采用比较简单的阈值分割法进行森林-非森林分类，而且没有对分类结果进行去类别噪声处理。虽然CBFC变化检测效果不好，但可以满足CBFC-EM-MRF及本研究方法对初始森林覆盖变化图的要求，这2种方法都可以检测到主要的森林覆盖变化、未变化区域，本研究方法的检测图与参考图最为接近，效果最好。

为进一步评价本研究方法的有效性，图 8展示了2块较大的变化区域与对应2期Landsat-5影像的对照。可以看出，本研究方法检测到的森林覆盖变化区域面积和轮廓与Landsat-5影像保持了最好的一致性，而其他2种方法都存在不同程度的虚警或漏警现象，进一步说明了本研究方法的有效性。

3.4 森林覆盖变化检测精度评价

3种森林覆盖变化检测方法的检测精度见表 2。由于试验区森林覆盖发生变化的区域很少、很小，森林覆盖变化参考图中变化像元只占总像元数的0.14%，剩下的99.86%为未变化像元，因此2个类别参考样本存在严重数量不均衡现象，参考检测正确率的定义，会导致3种检测方法的变化检测正确率都偏高，即使漏警率都在10.00%以上，正确率也都能达到85.00%以上。这种情况下，虚警率、漏警率更能说明变化检测方法的精度，尤其漏警率是用户最希望尽量降低的指标。CBFC只利用了HV极化后向散射特征，受斑点噪声影响较严重，也没有利用到影像空间上下文信息，因此检测正确率较低，虚警率和漏警率高。虽然CBFC-EM-MRF与本研究提出的EM-MRF方法相比正确率较高、虚警率较低，但差别很小，不到1.00%，而漏警率却比本研究方法高出36.91%，且本研究检测结果中，变化区域的“空洞”现象和个别孤立像元大大减少，降低了斑点噪声引起的“椒盐”现象，充分说明本研究方法具有更高的森林覆盖变化检测能力。

为更精细地评价本研究方法的可靠性，以参考数据中的13个森林覆盖变化区域(RFCN)为基准，分别统计每个变化区域内参考数据的森林覆盖变化像元个数(FCN)和3种变化检测结果中检测到的变化像元个数，评价3种算法的正确率和平均正确率，如表 3所示。可以发现，CBFC和CBFC-EM-MRF方法仅在第1、11、12个变化区域具有较好的检测精度，在剩下的10个变化区域中正确率均较低，大量的森林覆盖变化区域未检测到，与表 2中的较高漏警率相一致；而EM-MRF方法仅未检测到第10、11个变化区域，这2个区域内的像元个数较少，即使未检测到但对整体的森林覆盖变化检测影响并不大，平均正确率均高于CBFC和CBFC-EM-MRF，显示出本研究方法的优越性。

4 结论

本研究提出一种将“2期分别分类森林覆盖变化检测法”(CBFC)与“贝叶斯最大期望-马尔科夫随机场(EM-MRF)变化检测法”相结合的森林覆盖变化检测方法，采用阈值分割法分别对2期多极化SAR影像进行森林-非森林分类得到初始森林覆盖变化图，作为后续对多极化比值影像进行Fisher特征变换和EM-MRF分类的训练数据。以黑龙江省逊克县为试验区，对该方法的有效性进行评价，主要结论如下：

1) 2个极化比值差异特征经Fisher特征变换后得到新的变化检测差异特征可有效增强森林覆盖变化、未变化类别的对比度，且2个类别的像元值服从对数高斯概率密度分布，为EM-MRF分类取得稳定、有效的结果奠定了理论基础。

2) 由于SAR影像受斑点噪声的影响，只采用阈值分割法难以取得较好的变化检测效果，而通过马尔可夫随机场算法加入影像空间上下文信息以后，不仅提高了所提取森林覆盖变化区域的空间连贯性，而且大大降低了虚警率和漏警率。

3) 定量评价结果表明，本研究森林覆盖变化检测方法虚警率为1.58%，漏警率为11.87%，正确率为98.36%，检测精度明显优于CBFC和CBFC-EM-MRF。

本研究方法在对2期多极化SAR影像进行森林-非森林分类时，需要用户通过统计分析确定影像分割阈值，除了该步骤外，其他处理过程都可实现自动化，所以可认为是一种半自动化方法。该半自动化多极化星载SAR森林覆盖变化检测方法算法具有收敛性好、检测可信度高、需要用户交互较少等特点，对我国高分三号及未来其他相关多极化SAR卫星的森林资源监测业务应用具有重要参考价值。下一步应进一步分析该方法对初始化变化检测图的敏感性，提高森林-非森林分类步骤的自动化程度; 另外，本研究试验区覆盖范围较小，发生变化的森林地块较少，变化面积也较小，所提出方法对大区域森林覆盖变化检测的适用性还不明确，尚需更多试验区、更大覆盖范围的试验研究。