森林生物量是森林质量的重要体现，全面掌握森林生物量动态变化，对了解森林的生长过程、完整分析植被生态系统自然现象的变化、实现生态系统的有效修复具有重要意义。鉴于森林生物量、特别是区域范围内森林生物量准确计算的重要意义，近年来遥感技术被广泛应用于森林生物量的估算中，并且成为目前遥感和林业应用领域的研究热点之一^[1-3]。由于遥感技术仅能反映森林地上生物量（aboveground biomass，AGB）的变化，因此，基于遥感手段的森林生物量反演均指森林AGB。现今使用的遥感手段中，合成孔径雷达（synthetic aperture radar，SAR）由于克服了传统光学遥感受多云、多雾天气影响的局限，并具有获取森林冠层水平、垂直结构等特征的能力，成为当下森林AGB反演最具潜力的遥感数据源之一。随着SAR数据源的丰富，SAR数据可以提供的信息类型越来越多。不同波长的SAR数据对森林冠层的穿透性不同，有可以反映冠层散射的X、C波段，也有能够穿透冠层与树干、树枝发生作用的L、P波段^[4-8]；协同多频段SAR数据森林监测的优势，可实现多维度森林AGB的监测，不仅可以较全面地了解森林的SAR散射特征，还可有效提高森林AGB遥感反演的精度和饱和点。

森林AGB反演中常用的参数模型包括经验模型、物理模型和半经验半物理模型。Attema和Ulaby提出的水云模型（water cloud model，WCM）最初常用于农作物生物量、叶面积指数的反演^[9]。WCM属于经典的半经验半物理模型，既具有经验模型通过简单的线性或非线性回归来反演植被参数的优势，又可通过物理散射机制还原植被的垂直结构。WCM因其兼具模型简单和物理可解释性的优势，被广泛用于农作物长势监测参数反演^[10-12]。Santoro等^[13]将WCM进行改进，提出了干涉水云模型，并采用L波段SAR数据对瑞典北方森林蓄积量进行了反演试验，取得了较高的反演精度，在森林蓄积量高达350 m³/hm²的森林覆盖区，小班水平的误差仅为10 m³/hm²。戴玉芳等^[14]构建了我国东北针叶林AGB的Yamaguchi极化分解散射分量直方图，发现使用L波段反演我国东北针叶林AGB满足水云模型使用的前提假设，为后续利用WCM对该地区森林进行AGB反演奠定了基础。

受限于SAR数据获取的困难性，目前使用SAR信息及水云模型进行森林AGB的研究还多集中在单极化或单频段SAR数据的后向散射信息，对于全极化信息、不同频段极化信息的利用及其反演结果的对比研究还开展较少。此外，WCM多用于农作物参数反演，用于森林AGB反演的研究目前在国内鲜有报道，已有研究多集中在国外的研究区，且研究模型仅使用了单极化的后向散射信息或干涉信息^[15-20]。不同频段SAR数据对森林AGB水平、垂直分布表征能力差异明显，全极化SAR数据对森林散射体的形状、方向敏感，因此采用WCM、并将极化信息引入WCM，对比不同频段、极化信息在WCM中的AGB反演潜力，可为未来将多频SAR数据有效融合于WCM提供科学依据。因此，本文利用X、C、L、P波段全极化数据，首先使用WCM模型进行了AGB反演，对比了极化和波长对WCM反演结果的影响；然后将极化信息采用地体散射比的形式引入WCM，发展了极化水云模型（PolWCM），并探索了地体散射比参数在X、C、L、P波段进行AGB反演的可行性。

1.   研究区概况、数据获取及预处理

1.1   研究区概况

本文研究区为内蒙古大兴安岭国家级森林生态系统定位站，该生态站（121°20′ ~ 121°36′E，50°50′ ~ 50°59′N，图1）位于内蒙古呼伦贝尔盟的根河市。研究区地势相对平缓，区域内80%的坡度小于15°，海拔高度分布在800 ~ 1 200 m。研究区森林覆盖率大于75%，森林AGB平均值在50 t/hm²左右。研究区为典型的天然兴安落叶松（Larix gmelinii）林，其主要树种为兴安落叶松，有零星的白桦（Betula platyphylla）、樟子松（Pinus sylvestris var. mongolica）以及山杨（Populus davidiana）等。

图  1  研究区地理位置

Figure  1.  Geographic location of the study area

下载: 全尺寸图片 幻灯片

1.2   X、C、L、P波段极化SAR数据获取及预处理

1.2.1   SAR数据获取

本研究共获取到机载或星载的X、C、L、P波段4种频率全极化SAR数据、机载的LiDAR（light setection and ranging）数据和SRTM DEM（shuttle radar topography mission digital elevation model）辅助数据（各数据的覆盖范围见图2a）。其中X、C和L波段SAR数据为星载数据，分别为TerraSAR-X（图2b）、RADARSAT-2（图2c）和ALOS-2 PALSAR-2（图2d）数据；P波段SAR数据由国产机载SAR系统获取（图2e）。

图  2  X、C、L、P 4个波段覆盖范围概况图及地理编码后Pauli RGB

HH和HV分别为HH和HV两个极化后向散射系数。HH and HV are polarization backscattering coefficients of HH and HV, respectively.

Figure  2.  Overview of the coverage of SAR data and their geocoded images displayed as Pauli RGB at X-, C-, L- and P-band

下载: 全尺寸图片 幻灯片

1.2.2   SAR数据预处理

4个波段SAR数据均为全极化数据，其荷载方式及星座、所属国家、成像时间、中心入射角、距离向和方位向采样等参数见表1。X、C、L和P波段SAR数据预处理的共同步骤包括：单视复数图像转换为S2矩阵、辐射定标、滤波处理和地理编码。处理过程使用的SAR数据处理软件包括Gamma、PolSARpro等，同时结合部分IDL、Python等编码语言。由于各频段SAR数据获取的传感器略有差异，且搭载传感器的平台又有星载和机载之分，预处理过程中会略有差异，如RADARSAT-2 C波段数据在Gamma读入时已做定标处理，其他处理过程中的多视视数、定标数值、滤波方法等见图3。

表  1  根河试验区X、C、L、P波段数据详细参数

Table  1.  X -, C -, L -, P-band data parameters in Genhe Experimental Area

荷载及星座 Load and constellation	所属国家 Nation	获取时间 Acquisition time	中心入射角 Central incidence angle/（°）	距离向分辨率 Range resolution/m	方位向分辨率 Azimuth resolution/m
TerraSAR-X X-band	德国 Germany	2015−08−19	41.41	0.91	2.68
RADARSAT-2 C-band	加拿大 Canada	2013−08−20	46.09	4.96	4.73
ALOS-2 PALSAR-2 L-Band	日本 Japan	2014−08−29	36.52	2.86	2.64
Airborne CAMSAR P-Band	中国 China	2013−09−13	50.06	0.63	0.67

下载: 导出CSV 

| 显示表格

图  3  多频SAR数据预处理技术路线

“――”代表数据提供方已作辐射定标，本文不需再作定标。“――” means that the data provider has made radiometric calibration and no further calibration is required in this paper.

Figure  3.  Flowchart of multifrequency SAR data preprocessing technology

下载: 全尺寸图片 幻灯片

1.3   机载LiDAR数据获取与预处理

本文采用LiDAR提取的DSM（digital surface model）、DEM数据用于各频段极化SAR数据的地理编码，由LiDAR数据生成的森林AGB产品用于多频SAR数据WCM和PolWCM的训练及反演结果的验证。本研究获取的机载LiDAR数据是由中国林业科学研究院资源信息所提供，DEM（图4a）、飞行区域森林AGB（图4b）和研究区森林AGB（图4c）数据经插值后分辨率均为2 m，DEM、AGB数据的详细生成方法见文献[21]和文献[22]。

图  4  LiDAR提取研究区相关数据

Figure  4.  LiDAR extracts relevant data of the study area

下载: 全尺寸图片 幻灯片

本文以研究区LiDAR森林AGB图为基础，按照750 m的空间采样间隔，在ArcGIS软件中运用ArctoolBox模块中随机空间采样功能，并在此基础上增加人工干预（去除道路及裸露地物）共选取相应样点113个（图4c）。113个森林AGB样本的平均AGB为46.70 t/hm²；0 ~ 10 t/hm²的采样点有3个；10 ~ 20 t/hm²的采样点有14个；20 ~ 30 t/hm²的采样点有17个；30 ~ 40 t/hm²的采样点有17个；40 ~ 50 t/hm²的采样点有18个；50 ~ 60 t/hm²的采样点有14个；60 ~ 70 t/hm²的采样点有10个；70 ~ 80 t/hm²的采样点有8个；80 ~ 90 t/hm²的采样点有7个；> 100 t/hm²的采样点5个。并以此用于本文中各频段WCM、PolWCM模型训练及反演结果验证。其中建模的样本为60%，验证的样本数量为40%。

2.   研究方法

2.1   X、C、L、P各波段SAR特征提取

分别提取各波段6个SAR特征，具体包括Freeman 分解3分量，即森林冠层的体散射分量（F_Vol）、地面与树干的二次散射分量（F_Dbl）、地表散射分量（F_Odd）；另外3个特征为HH、HV、VV后向散射系数。

2.2   水云模型

WCM是以农作物为研究对象提出的植被微波散射模型。该模型对植被覆盖层的散射机制进行了简化，假定：“云”代表植被层，由类似水分子的相同大小、均匀分布在整个植被空间的各向同性颗粒组成；该模型仅考虑植被冠层与地表直接散射，忽略多次散射；植被层变量包括和“云”的水分含量成比例的高度和密度两个变量。本文参考Santoro等^[13]提出的用于瑞典北方森林蓄积量反演的水云模型，将农作物替换为森林，即将森林覆盖区的后向散射能量由冠层直接散射σ_forest和地表直接散射σ_ground两个部分通过式（1）来表达，其中β为经验参数， V_forest表示森林AGB。通过式（2）到式（6）的推导得出由微波后向散射理论值作为自变量、森林AGB作为因变量的森林AGB WCM反演模型。本文分别将不同极化通道的后向散射系数带入式（6）表示的WCM中的σ来反演森林AGB。式（6）中的β₀、β₁和β₂参数通过最小二乘法拟合获得。

2.3   将极化信息引入水云模型

由于在传统的WCM中没有考虑极化信息，随着全极化数据的丰富，其在具体应用中受到限制，为了拓展WCM的应用，本文对经典WCM进行了改进，将极化信息通过地体散射比的方式引入WCM，并将其用于森林AGB的反演。具体通过极化分解参数构建地表散射与体散射分量比值μ，并将其代入式（1），来发展传统的WCM，称其为PolWCM，其推导过程见式（7） ~ （9）。

2.4   模型验证方法

本章采用两种模型（WCM和PolWCM）对研究区的森林AGB进行了反演，反演结果精度的定量评价通过反演结果与真值之间的R²和RMSE表征。这2个参数的计算公式如下：

式中： ${y_i}$为第$i$个实测值，${\hat y_i}$为第$i$个预测值，$\bar y$为实测值的平均值，n为样本总数。

3.   结果与分析

3.1   WCM反演森林AGB

分别采用各波段HH、HV、VV极化的后向散射系数来替代WCM中的参数σ，使用最小二乘法进行拟合，获取WCM中$\; {\beta _0}$、$\; {\beta _1}$和$\;{\beta _2}$，研究区森林AGB的反演结果见表2和图5。除X波段外，以HV通道的后向散射系数作为总后向散射输入WCM，反演的森林AGB精度最佳，其中L波段反演结果的R²值最高，RMSE值最低，分别为0.46和18.00 t/hm²。其次是P波段HV极化的反演，R² 为0.43，RMSE为21.18 t/hm²。该研究结果表明WCM更适用于采用波长较长波段（L和P）的HV极化后向散射系数进行森林AGB反演，对于波长较短的X波段则反演结果精度较差。值得注意的是，采用L波段的HH和VV的后向散射系数代入WCM进行森林AGB反演，反演结果精度仍高于其他波段HV极化后向散射系数代入WCM的反演结果，这可能是由于研究区的森林AGB平均值约为50 t/hm²，略低于L波段的饱和点，因此L波段更能表征研究区森林的微波散射特征。

表  2  X、C、L、P波段WCM反演结果

Table  2.  Inversion results of WCM for X-, C-, L- and P-band

波段 Band	极化通道 Polarization channel	R2	RMSE/(t·hm−2) RMSE/(t·ha−1)
X	HH	0.03	26.12
	HV	0.01	26.49
	VV	0.05	25.70
C	HH	0.19	24.32
	HV	0.33	23.32
	VV	0.30	22.76
L	HH	0.24	21.27
	HV	0.46	18.00
	VV	0.30	20.34
P	HH	0.12	27.40
	HV	0.43	21.18
	VV	0.09	27.04

下载: 导出CSV 

| 显示表格

图  5  X、C、L、P波段WCM反演森林AGB

散点图中的红色实体线为森林AGB估测值拟合线，黑色点为训练样本，蓝色点为验证样本。图中RMSE单位为t/hm²。下同。The red solid line in the scatter diagram is the fitting line of forest AGB estimation values, the black points are the training samples, and the blue points are the verification samples. The unit of RMSE is t/ha. The same below.

Figure  5.  X-, C-, L-, P-band WCM inversion results of forest AGB

下载: 全尺寸图片 幻灯片

从图5中可以看出X波段各极化的后向散射系数分布较松散，没有明显的规律性，因此反演的结果也较差。相比X波段，C波段样点的聚集性明显增强，但是饱和点较低。L波段各个极化的后向散射系数随森林AGB动态变化呈一定的规律性，并且饱和点也明显提高，特别是HV极化，样本点的分布非常接近拟合的模型曲线。P波段除HV极化外，HH和VV极化后向散射随森林AGB的变化也未有明显的规律性，尽管总后向散射值的动态变化范围相比X、C波段有显著的增加，但也不能较好地表示森林AGB的动态变化。

3.2   PolWCM算法反演森林AGB

基于样本数据，将$\mu$采用Freeman-Durden分解的地表散射分量（F_Odd）与体散射分量（F_Vol）的比值来代替，通过最小二乘法拟合，将比值代入极化水云模型进行森林AGB反演。其反演结果见表3和图6。表3总结了采用PolWCM在4个波段进行森林AGB反演的结果。PolWCM的反演结果略优于WCM，特别是在X、L和P波段，RMSE的值都明显降低。但是R²的值除P波段外，都明显降低。 图6描述了将极化地体散射比引入WCM进行森林AGB反演的结果精度，在X和C波段，采样点对应的4个频段SAR数据的极化地体散射比的比值分布不均匀，在X波段中$\mu$值多分布在0.20 ~ 0.25之间；而在C波段，该值则多分布在0.10 ~ 0.20之间，同等AGB动态变化范围内，在C波段地体散射比的动态变化范围略大于X波段，因此，C波段的反演结果精度略高于X波段，但极化地体散射比与森林AGB均表现出较弱的相关性。相比波长较短的X和C波段，L和P波段的动态变化范围则显著增宽，在L波段$\; \mu$值的范围约在0.01 ~ 0.55之间，在P波段则略宽，范围约为0.05 ~ 0.65，因此，反演结果精度与WCM反演结果类似，较长波段反演结果精度较优。采用$\; \mu$描述的PolWCM在X和C波段出现明显的饱和现象，且饱和点均较低；但在波长较长的L和P波段，则饱和现象不明显。

表  3  X、C、L、P 波段PolWCM反演森林AGB结果

Table  3.  Inversion forest AGB results of PolWCM for X-, C-, L- and P-band

参数 Parameter	波段 Band	R2	RMSE/(t·hm−2) RMSE/(t·ha-1)
FOdd/FVol	X	0.00	24.90
	C	0.01	24.71
	L	0.40	17.70
	P	0.51	18.08
注：FOdd和FVol分别为Freeman极化分解三分量中的地表散射与体散射分量。Note: FOdd and FVol are the surface scattering and volume scattering components of Freeman polarization decomposition.

下载: 导出CSV 

| 显示表格

图  6  X、C、L、P波段PolWCM反演森林AGB结果

F_Odd和F_Vol分别为Freeman极化分解三分量中的地表散射与体散射分量。F_Odd and F_Vol are the surface scattering and volume scattering components of Freeman polarization decomposition.

Figure  6.  Forest AGB retrieved by PolWCM at X-, C-, L- and P-band

下载: 全尺寸图片 幻灯片

综合WCM和PolWCM的反演结果可知，WCM在森林AGB反演中具有极化和波段依赖性。在所有用于反演的不同极化通道的后向散射系数中，HV极化后向散射系数引入WCM可获得最优的反演结果，且在C、L和P 3个波段中均表现出较好的反演能力；其次，WCM在波长较短（X和C）的波段，反演效果较差，而在长波长（L和P）中反演效果较好。PolWCM同样具有波长依赖性，在波长较长的L和P波段中反演结果精度较高；同时，其在L和P波段的反演结果优于同等条件下WCM的反演结果，R²分别为0.40和0.51，RMSE分别为17.70和18.08 t/hm²。

4.   讨论与结论

4.1   讨　论

WCM于20世纪90年代初被用于森林生物量或蓄积量的反演，使用的SAR数据包括机载和星载数据，频段主要为X、C和L波段。Pulliainen等^[23]首次基于机载高分辨率SAR数据构建了WCM用于森林AGB反演，通过将X和C波段的HH和VV双极化后向散射系数代入WCM用于森林AGB的反演的研究，结果表明：C波段VV极化后向散射系数与森林AGB相关性最高，适用于WCM反演森林AGB的使用条件^[23]；该研究结果与本文X和C波段的HH、VV极化后向散射系数反演森林AGB的结果高度一致，反演结果中C波段优于X波段，VV极化优于HH极化。Santoro等^[13]、Askne等^[15]、Cartus等^[16]考虑森林的垂直和水平结构变化，将郁闭度因子、干涉信息引入最初的WCM，并将其用于C波段、L波段星载SAR数据进行森林蓄积量的反演，结果指出该模型在样地水平反演精度较高、在像元水平反演精度较低。而Cartus等^[16]在L波段的研究表明，HV后向散射的反演结果优于HH极化，这与本文研究结果一致。戴玉芳等^[14]构建了我国东北针叶林AGB的Yamaguchi极化分解散射分量直方图，得出L波段以体散射分量为主导分量，而二次散射极化散射分量可以忽略不计，满足WCM反演我国东北针叶林AGB的前提需求；本文采用L波段数据构建WCM反演森林AGB精度均较高的结果也说明了这一点。Santoro等^[17]综述了L波段WCM在瑞典森林蓄积量反演中的研究，指出了WCM森林蓄积量反演中的极化依赖性；本文研究同样表明了WCM的极化依赖性。Cartus等^[24]采用光学数据作为辅助数据训练WCM获得模型中的地表散射和植被散射分量，并将其用于X、C、L和P波段区域尺度森林AGB反演，结果与本文相似，即：P > L > C > X；该研究同时指出即使组合多频信息精度也未高于P波段信息单独反演的结果；另外，该研究中各频段的极化依赖性受环境因子、待估森林结构等影响明显。

将除后向散射信息以外的其他SAR信息引入WCM，可以有效提高森林AGB的反演精度。已有研究中，包括将干涉信息、时间序列信息等引入WCM。Santoro等^[13]、Askne等^[15]、Cartus等^[16]将干涉相干性信息引入WCM，发展了干涉水云模型，森林蓄积量和生物量的反演精度和饱和点明显提高。以上国内外研究中均未使用极化信息改进WCM，这也表明了本文将4个波段的极化信息以地体散射比的方式引入WCM形成PolWCM的必要性，从试验结果来看，将WCM引入极化信息保证了模型的散射特征和简单化，在总体上未出现明显饱和点的条件下，可明显提高森林AGB反演精度。同时，本文研究结果也印证了戴玉芳等^[14]提出的L波段较适宜使用水云模型反演中国北方针叶林AGB的假设。与Santoro等^[24]和Kumar等^[25]的研究结果相比，本文反演结果R²虽未达到其精度，但RMSE值总体上均明显高于这些研究结果。

4.2   结　论

本文通过将X、C、L和P 4个波段HH、HV和VV 3个极化方式的后向散射系数分别代入WCM进行研究区的森林AGB反演，同时以极化分解参数为基础，构建地体散射比，并将其应用于WCM来反演森林AGB，探索了极化方式、波长、极化信息引入方式等对基于WCM进行森林AGB反演的影响方式及程度。研究得出以下结论：

（1）WCM进行森林AGB反演具有极化依赖性和波长依赖性。

在森林AGB平均水平较低时，例如本文研究区，平均值约为50 t/hm²时，在所有波段中，将HV极化后向散射系数代入WCM进行森林AGB反演，反演结果多优于同波段的其他极化通道后向散射系数的反演结果；在森林AGB大于50 t/hm²时，HV或VV极化后向散射系数反演结果较优；在所有森林AGB水平上，HH极化后向散射系数反演结果精度均最差。WCM的反演精度在后向散射系数所使用的波长较长时反演结果精度较优，在低森林AGB水平时，L波段优于P波段；但是生物量水平较高且研究尺度较大时，P波段优于L波段。

（2）将极化信息以地体散射比的形式引入WCM进行森林AGB反演可明显提高反演精度。

将极化信息采用地体散射比形式（F_Odd/F_Vol）引入WCM推导形成PolWCM，用于X、C、L和P 波段森林AGB的反演，R²值有明显提高，RMSE值有明显降低，总体反演精度明显提高。

本文中WCM模型的参数通过最小二乘方式获得，未考虑其物理意义，后续研究中，可以引入辅助数据，考虑模型中各参数的物理意义，进而探索其在增强物理可解释性的基础上提高反演精度的可能性。其次，本文中P波段为机载极化SAR数据，在与X、C和L波段对比时，并未考虑机载数据与星载数据的差异，由于中心入射角变化对机载和星载SAR数据的影响差异明显，未来研究中需要进一步考虑机载和星载数据差异对森林AGB反演结果的影响。