森林生态系统是森林群落及其环境在功能流作用下形成的一定结构、功能和自调控的自然综合体，是陆地生态系统中面积最大、最重要的自然生态系统，在全球碳循环、控制温室气体以及维持地球生物圈稳定等方面具有非常重要的地位 (Houghton，2005)。当前，温室效应日益严峻，世界各国正加强国家之间、洲际之间的合作，扩大森林面积、促进森林可持续经营增加碳汇以应对全球气候变化。为寻求工业发展与生态平衡的折中点，精确地估算森林碳储量和生物量显得十分必要。近年来，随着遥感技术的不断进步, 运用遥感技术进行森林生物量估算已成为一种新兴的研究手段。与传统的实地勘察测量方法不同，遥感技术可以实时、准确、无破坏地对森林生物量进行估算，进而推算碳储量，且易于从小尺度推广到大区域，尤其对一些地面测量难以开展的偏远天然林来说，遥感估算的优势更加凸显。

目前，在对森林生物量进行估测的遥感方法中，主要使用微波遥感、激光雷达以及光学遥感等技术。微波遥感具有一定的穿透能力，如合成孔径雷达 (synthetic aperture radar，SAR) 可与森林生物量的主体即树干和树枝发生作用 (陈尔学等，2007)，但受地形起伏影响较大，且当生物量达到一定水平时后向散射系数趋于饱和。基于极化干涉SAR的极化相干层析技术通过提取森林相对反射率垂直分布，将拟合得到的层析测量树高建立生物量模型，有助于克服高值生物量估测的信号饱和问题 (李文梅等，2014)，但目前由于数据源的局限，尚无法大面积应用。激光雷达技术 (light detection and ranging，LiDAR) 在林木高度测量和森林垂直结构信息获取方面能实现更高效、更精准的监测，具有其他遥感技术无可比拟的优势 (庞勇等，2005)；同时，激光雷达反演的高精度森林结构参数还可用作其他遥感手段的训练数据和验证数据 (庞勇等，2005；2012；郭庆华等，2014)。

相对于微波遥感、激光雷达而言，光学遥感技术发展时间更长、数据源更丰富。随着国内外各类星载与机载传感器技术平台的不断更新与发展，高分辨率光学影像的获取更加便捷。高分辨率遥感影像可以充分利用丰富的光谱信息、植被指数信息、纹理信息以及空间信息等特征反映更详细的森林遥感信息 (李德仁等，2012)。基于光学遥感影像估算生物量通常有2种方式：一是直接以遥感波段作为自变量，使用单波段或多波段与生物量建立估算模型；二是将遥感影像的各波段经过数学运算组合成不同形式的遥感指数，以这些遥感指数直接或间接参与建模 (徐新良等，2006)。纹理信息作为遥感影像的派生特征，可揭示图像内容的水平结构信息并反映图像内容的空间变化规律，逐渐被加入到生物量估算模型中参与反演。Lu等 (2005)对巴西亚马逊流域热带森林的TM影像数据与灰度共生矩阵 (gray level co-occurrence matrix，GLCM) 的8个纹理参数进行相关性分析，结果显示GLCM中的均值与相关性因子对森林生物量预测贡献度较高。Sarker等 (2011)利用ANVIR-2数据对中国香港地区森林生物量进行预测，并比较5种参量 (波段反射率、波段比、植被指数、纹理指数和纹理指数比) 对生物量的预测效果，结果显示纹理指数比方法对香港地区亚热带、热带森林的生物量估算效果最优。尽管纹理信息在预测生物量上很有潜力，但是实际应用中却不超过60%，甚至在对结构简单的温带森林的应用上也很少 (Sarker et al., 2011)。法国学者Couteron等 (2005)在估算法属圭亚那地区热带雨林生物量时提出了FOTO (Fourier-based textural ordination) 纹理方法，该方法是一种基于傅里叶变换思想的新型纹理提取技术，研究显示其在对高密度热带植被的生物量估算中没有出现饱和现象。Proisy等 (2007a)利用FOTO方法对法属圭亚那地区的红树林生物量进行了估算，从影像上提取出森林结构参数与垂直参数并构建了红树林的生物量模型 (Proisy et al., 2011)。Singh (2012)在估测马来西亚婆罗洲森林生物量时指出，由于大气和信号饱和的影响，用树冠顶层的反射值对生物量进行预测其有较大不确定性，利用纹理信息则更加有效，并发现FOTO纹理因子对生物量的估算精度优于GLCM纹理因子。

FOTO纹理参数提取研究在国外的典型热带林区陆续得到了不同程度的应用，而在国内尚未对该方法展开相应的研究与探讨。本研究以黑龙江省小兴安岭地区凉水国家自然保护区为研究区，开展FOTO纹理参数估算温带森林生物量的尝试，以期为提取新型纹理参数估算森林生物量提供新的参考途径。

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

黑龙江省凉水国家自然保护区 (47°11′N，128°53′E) 地处我国东北东部山地小兴安岭南坡达里带岭支脉的东坡，黑龙江省伊春市带岭区境内 (图 1)。保护区内地形复杂，以山地为主，森林覆盖率98%，海拔280~707 m，平均坡度10°~15°，为典型的低山丘陵地貌，具有明显的温带大陆性季风气候特征。年平均气温-0.3℃，年平均降水量676 mm，年平均相对湿度78%，无霜期100~120天。境内地带性土壤为山地暗棕壤。气候特点：春季风大；夏季气温较高，降雨集中；秋季降温急剧；冬季漫长，寒冷干燥而多风雪 (刘丽娟等，2013)。保护区自然资源丰富，植被群落类型复杂多样，既有处于演替顶极状态的原始阔叶红松 (Pinus koraiensis) 林、鱼鳞云杉 (Picea jezoensis var. microsperma) 林和兴安落叶松 (Larix gmelinii) 林，又有处于不同演替阶段的次生林。保护区内有大片较原始的红松针阔混交林，是我国目前保存下来最为典型和完整的原生红松针阔混交林分布区，也是我国和亚洲东北部很具代表性的温带原始红松针阔混交林区 (庞勇等，2012)。

1.2 遥感数据

遥感数据采用2009年9月4日和13日拍摄的机载航空影像 (图 1)，搭载国产运-5小型多用途飞机，携带激光扫描仪 (LMS-Q560) 和CCD相机 (DigiCAM-H/22) 进行同步数据采集。获取数据总覆盖面积约200 km²。CCD相机拍摄的南北航向航片经过正射校正和拼接处理后，得到完整覆盖研究区的高分辨率遥感影像。遥感数据空间分辨率为0.5 m，含红 (700 nm)、绿 (540 nm) 和蓝 (480 nm)3个可见光波段。

1.3 生物量数据

高精度机载激光雷达反演的生物量数据可作为其他遥感数据反演时的验证数据。本研究所用的森林地上生物量参考数据是由与航片同步获取的机载激光雷达数据与地面样地调查数据联合反演得到的 (图 2)，生物量影像空间分辨率为10 m，反演精度高于90%，因此可以作为可靠的参考生物量。机载CCD影像与激光雷达反演的生物量数据是经过同步飞行得到的，二者经过实验室系统检校和测区检校场检校，位置匹配精度优于0.3 m (庞勇等，2012)。

2 FOTO纹理计算方法和生物量建模 2.1 技术路线

如图 3所示，本研究基于FOTO纹理方法从高分辨率航片中提取FOTO纹理因子，使用机载LiDAR反演的森林地上生物量作为参考值进行建模，对比分析不同波段、不同计算窗口下FOTO纹理因子的生物量反演能力。采取以下2种方法估算生物量：

1) 对CCD影像的3个单波段影像，即红、绿、蓝波段分别采用FOTO算法，3个波段影像共提取得到9个纹理因子 (每个波段影像的PC1、PC2、PC3) 用于建立回归方程，并对森林生物量进行预测。

2) 对红、绿、蓝3个波段取平均生成一幅单波段灰度影像，通过FOTO算法从该波段平均影像中提取出3个纹理因子 (波段平均影像的PC1、PC2、PC3)，与生物量建立回归方程并对生物量进行预测。

2.2 FOTO纹理指数算法

FOTO纹理指数是一组将快速傅里叶变换算法与主成分分析结合提取的指数，包含森林冠层的周期信息及空间异质性。根据FOTO纹理指数算法原理，首先选择一个固定的窗口尺寸遍历整幅影像，将影像分割成若干方形小块。窗口大小W_S定义为 (单位：m)：

${W_s} = N \times S{\rm{。}}$

(1)

式中：N为X或Y方向上的像素点个数；S为图像的分辨率。

再对分割产生的每一小块影像进行二维快速傅里叶变换 (2D-FFT) 处理，傅里叶变换提供了一种将图像从空间域转换到频率域进行频谱分析的方法，将复杂信号看作不同频率、不同方向和不同振幅的正余弦信号的叠加。图像经过2D-FFT变换后得到的频谱是复数矩阵，复数的实部a_pq与虚部b_pq统称为傅里叶分解系数，是复杂的图像信号经过傅里叶变换后变成简单余弦信号和正弦信号的幅值：

${a_{pq}} = {N^2}\sum\limits_{x = 1}^N {\sum\limits_{y = 1}^N {{M_{xy}}\cos \left[ {2{\rm{ \mathsf{ π} }}\left( {px + qy} \right)/N} \right]} } ;$

(2)

${b_{pq}} = {N^2}\sum\limits_{x = 1}^N {\sum\limits_{y = 1}^N {{M_{xy}}\sin \left[ {2{\rm{ \mathsf{ π} }}\left( {px + qy} \right)/N} \right]} } {\rm{。}}$

(3)

对于一个固定的正余弦振动方向，频率p或q分别代表正余弦波形在影像内的重复次数。假设给定的影像为n×n维灰度影像，根据奈奎斯特频率可得到最高正余弦频率为p_max=n/2和q_max=n/2。FOTO中计算的傅里叶径向频谱 (Fourier radial spectra，简称r-spectra)(Couteron, 2002; Couteron et al., 2006) 用来表征影像频率环内各方向的平均能量，r表示给定某一范围的影像内正余弦波在该范围内拟合出的树冠周期次数。在本研究中，选择前29个傅里叶谐波能量 (Ploton et al., 2012)。随着谐波频率逐渐趋于奈奎斯特频率，幅值越来越小并趋近于0。

$I\left( r \right) = \left[ {\sum\limits_\theta {\left( {a_{pq}^2 + b_{pq}^2} \right)} } \right]/\left( {k{\sigma ^2}} \right);$

(4)

$\left\{ {\begin{array}{*{20}{l}} {\quad {r^2} = {p^2} + {q^2}}\\ {\theta = \arctan \;\left( {p/q} \right)\;} \end{array}} \right\}{\rm{。}}$

(5)

式中：r=1, 2, …, 29；θ=0°~360°；k为每个频率圆环内点的个数；σ²为图像方差。

图像r-spectra反映了图像频率域中每一谐波频率环内的能量大小，r-spectra曲线可作为描述森林中树冠周期的指标。

在一定区域内，如果该区域中树冠较大的成熟林木占主要成分，则在图像中表现为树冠纹理周期长，其对应的r-spectra曲线会在中低频区域有较高的数值，高频区域的数值较低，曲线呈陡峭锯齿状；若窗口内以树冠较小的幼年林木为主，则影像中树冠的纹理周期较短，r-spectra的低频和高频区域数值无明显起伏，曲线略显平缓 (Couteron et al., 2006；Singh，2012；Ploton et al., 2013)。采用120 m×120 m窗口从原影像中随机选出4种不同尺寸的树冠图像 (图 4)，从计算所得的r-spectra曲线可看出，当林木树冠较大时，波峰会在中低频区域出现，当林木树冠较小时，则会在中高频区域有较高值，符合上述描述规律。

根据FOTO算法得出所有窗口影像的r-spectra值后放入同一矩阵内，对该r-spectra矩阵进行标准主成分分析 (principle component analysis, PCA)。结果显示，前3个主成分分量占总体变量成分的84%以上，故最终采用前3个主成分作为从图像中提取的纹理因子。

2.3 基于FOTO纹理指数的生物量建模分析

将基于FOTO算法从影像中提取的纹理因子作为自变量，激光雷达反演的参考生物量作为因变量，采用多元逐步回归方法建立生物量模型。进行生物量建模时采取十折交叉验证，即每次取其中9份样本与生物量进行多元逐步回归分析，剩余1份样本进行反演精度检验，该过程进行10次直至遍历全部样本 (Picard et al., 1984)。在某些复杂情况下，交叉验证可获得预测误差的近似无偏估计 (Efron et al., 1983)，十折交叉验证可充分利用样本有效信息获得稳定可靠的估测模型，并可在一定程度上避免过拟合问题。

模型评价指标通常采用R²与RMSE，十次交叉验证结果的R²与RMSE平均值作为最终模型评价指标。R²为决定系数，值在0~1之间，当自变量对因变量的解释能力越强时，R²越接近于1；RMSE为均方根误差值，RMSE越小，表示模型预测能力越好。

3 结果与分析 3.1 FOTO纹理因子的提取

根据FOTO原理，可从每个单波段影像中提取得到3个纹理因子，图 5为CCD影像中红波段提取到的FOTO因子PC1、PC2和PC3(窗口尺寸选择100 m×100 m)，其余波段提取的纹理影像类似于红波段。将CCD影像的9个FOTO纹理因子和波段平均影像的3个FOTO纹理因子分别与森林生物量进行相关性分析，统计结果如表 1所示，所有统计值都在0.01水平上显著。从表 1可看出，FOTO方法提取的纹理因子与生物量的相关性高低排序为PC2、PC1、PC3，且相关系数均大于0.3，其中蓝波段的PC2与生物量的相关系数最高，为0.636，红波段的PC3与生物量的相关系数最低，为0.309。PC1与森林生物量在0.01水平上呈显著负相关，而PC2和PC3与森林生物量在0.01水平上呈显著正相关。

3.2 纹理窗口尺寸选择

FOTO纹理因子计算时首先需要选择计算窗口尺寸，且窗口尺寸的选择直接影响最终生成的生物量图的空间分辨率。本文在FOTO方法基础上分析采用5种不同窗口尺寸 (60 m×60 m，80 m×80 m，100 m×100 m，120 m×120 m和150 m×150 m) 对生物量建模与反演的影响。通过对激光雷达生物量图进行相似性区域分割，将分割斑块均值作为提取的生物量结果。CCD影像中提取出的9个纹理因子经过多元逐步回归与生物量建立估测模型，结果如表 2所示。

从表 2可看出，5种窗口尺寸下得到的纹理自变量与生物量的R²均达到0.67以上，选入模型的变量均在0.01水平显著。这表明FOTO纹理因子对生物量的解释能力较强，且R²随着窗口尺寸变大从0.67增加到0.81，RMSE从66.98 t·hm^-2减少至46.78 t·hm^-2。图 6显示了5种窗口尺寸下FOTO预测生物量值与实际生物量值的十折精度检验结果，窗口尺寸为60 m×60 m，80 m×80 m，100 m×100 m，120 m×120 m和150 m×150 m时对生物量的估测精度分别为67.3%，73.4%，74.4%，78.3%和80.9%。图 6显示对于生物量值小于100 t·hm^-2的区域出现高估现象，在生物量大于400 t·hm^-2时开始存在不同程度的低估，60 m×60 m和80 m×80 m窗口尺寸下的预测精度低于100 m×100 m以上窗口尺寸的结果，低值高估与高值低估的情况更多。

同时，利用上述5种不同窗口尺寸对波段平均影像进行FOTO处理，经过多元逐步回归将3个纹理因子与生物量建立估测模型并得出模型表达式。从表 3可知，模型R²在0.57~0.71之间，RMSE从77.12 t·hm^-2减少至58.51 t·hm^-2，选入模型的变量与生物量在0.01水平上显著。图 7显示了5种窗口尺寸下的反演精度检验结果，其预估精度分别为58.2%，62.1%，64.3%，67.4%和70.9%，同样显示出当窗口尺寸为60 m×60 m与80 m×80 m时，模型拟合精度较低，易出现更多的错误估算情况。

3.3 森林生物量反演结果

综合考虑生物量反演尺度与反演精度等因素，80 m×80 m~120 m×120 m之间为温带森林的最优FOTO窗口尺寸。利用提取的9个纹理因子参与多元逐步回归模型对试验区机载CCD影像进行反演，100 m×100 m窗口生成的反演公式为：

$\begin{array}{l} {\rm{AGB}} = 263.087 + 18247.724 \times {\rm{R\_PC1}} - 8162.608 \times \\ {\rm{R\_PC2}} - {\rm{11193}}{\rm{.077}} \times {\rm{R\_PC3}} - {\rm{18279}}{\rm{.432}} \times \\ {\rm{G\_PC1}} + 791.536 \times {\rm{G\_PC2 + 13864}}{\rm{.787}} \times \\ {\rm{G\_PC3}} - {\rm{1709}}{\rm{.723}} \times {\rm{B\_PC1}} + 1799.860 \times {\rm{B\_PC2。}} \end{array}$

(6)

式中：R_PC*，G_PC*，B_PC*分别表示红、绿、蓝波段的FOTO纹理因子。

生物量反演效果如图 8所示。凉水国家自然保护区平均生物量大于260 t·hm^-2，回归以及验证结果表明，FOTO方法提取的纹理因子与生物量相关性较好，预测精度为74.41%，RMSE为50.55 t·hm^-2。

3.4 结果分析

本研究对高分辨率航空CCD影像的3个波段以及波段平均影像分别采用FOTO算法提取纹理因子参与生物量建模并反演，分析了5种FOTO窗口尺寸对生物量反演结果的影响。结果显示，FOTO纹理参数均与森林生物量在0.01水平显著，且在5种不同窗口尺寸下，3个波段CCD影像提取的9个纹理因子对生物量的估测精度分别67.3%，73.4%，74.4%，78.3%和80.9%，波段平均影像提取的3个FOTO纹理因子对生物量的估测精度分别为58.2%，62.1%，64.3%，67.4%和70.9%。对本研究而言，CCD影像的9个FOTO纹理因子与生物量进行模型回归并参与反演的统计精度优于波段平均影像的3个纹理因子的结果。

CCD影像红、绿、蓝波段间存在的相关性使得提取的纹理参数在参与建模时会导致多重共线性问题，模型的预测功能减弱。经统计回归分析，模型自变量存在方差膨胀因子VIF>10的情况，即存在多重共线性关系。因参数估计量皆显著，且回归分析的目的仅用于预测，因此对此多重共线性可不做处理 (林清泉，2009)。波段平均影像则对红、绿、蓝波段进行光谱数值平均，生成的“拟全色”影像弱化了波段间的相关性，得到3个纹理因子参与回归，统计结果显示VIF接近1即所有解释变量几乎完全无关，生成的模型较稳定。在实际应用中，在满足精度评价指标RMSE与拟合程度R²的要求时，选择单波段进行FOTO运算能兼具精度与效率，同时预测模型比较稳定。对多波段进行FOTO运算能得到更好的统计精度，但伴随着更高的计算复杂度，且由于可能存在多重共线性关系会导致预测模型不稳定。

3.5 结果讨论

FOTO方法的应用始于解决热带地区高密度生物量的估算饱和问题，并在全球典型热带地区的生物量估算研究中取得了较好的估测效果，如准确估测了法属圭亚那地区的热带红树林和陆地热带雨林 (Proisy et al., 2007a; 2007b；2011)、南美亚马逊盆地的热带雨林 (Barbier et al., 2010)、中非地区的热带雨林、马来西亚婆罗洲地区的热带森林 (Singh，2012) 以及印度南部的西高止山脉热带植被 (Ploton et al., 2012) 的森林生物量以及其他森林结构参数，在某些样地生物量大于450 t·hm^-2的地区并未出现估值饱和的现象，且根据估算模型得出的预测精度较高，多数情况R²达到0.8以上 (Proisy et al., 2011)。

FOTO方法本质上是FFT与PCA的有效结合 (Singh，2012)。FFT提供了从频域角度分析信号的方法，尤其当影像中存在明显的树冠周期信息时，通过FFT方法分析影像的周期特征非常具有优势；PCA将r-spectra变量提取前N个重要特征量表征树冠大小信息。FOTO纹理因子通过将二者巧妙结合，把林木的树冠大小通过其在影像上表现出的纹理特征与森林生物量紧密结合起来，成为一种有效估算生物量的新方法。

FOTO运算过程中选择的窗口尺寸大小会影响生物量反演结果。Ploton等 (2012)指出，FOTO选择的最佳窗口尺寸应大于影像中最大树冠尺寸的5倍，因此该研究在估算热带森林时采用125 m (文中最大冠层估测为25 m) 的窗口。对凉水试验区样地调查结果分析显示 (Meng et al., 2016)，样地冠层平均尺寸为 (4.36±1.64) m，影像目视判读结果显示研究区内尚存在冠层直径为7.6~13.2 m的林木，且由于林木疏密度呈现的林木冠层与间隙的纹理表现不同，综合考虑选取80~120 m作为温带林试验区的FOTO窗口尺寸。与单木和林分追求的精准反演结果不同，FOTO是一种大尺度的生物量反演，当选择的窗口尺寸不同，FOTO输出结果的分辨率会随着窗口尺寸发生变化，并对反演精度产生影响：当窗口尺寸太小时，如果区域内存在特大冠幅林木，则会影响树冠结构周期性的分析结果；当窗口尺寸过大时，则会产生较大的估测尺度使得生成的生物量结果图的空间分辨率降低，并纳入冗余的地形纹理信息等 (Couteron et al., 2006)。在平衡建模精度与产品输出空间分辨率因素后，研究最终选择生成100 m分辨率的生物量预测图。

森林树冠结构的周期性越明显，则纹理参数与生物量的相关性越强 (Eckert，2012)。当森林中存在较大的林窗、林分间隙或者因自然干扰因素形成的倒木或斜木，则在分析影像时会出现树冠纹理的“伪周期”，给估算结果带来一定影响。且不同的太阳高度角、传感器拍摄角度以及地形起伏的影响，使得树冠间的阴影歪曲了真实的树冠周期，也会降低估算精度。

4 结论

本研究利用FOTO方法从高分辨率机载影像提取纹理参数，与森林生物量进行相关性分析，并采用多元逐步回归建立黑龙江省凉水国家自然保护区的森林生物量遥感估测模型。结果如下：

1) FOTO纹理对温带针阔混交林区的生物量反演精度较高，R²达0.81，RMSE为46.78 t·hm^-2，表明从影像中获取的FOTO纹理因子对北方温带混交林区的生物量反演有极大潜力。

2) 对于温带森林反演，可根据估算精度需求选择80~120 m作为FOTO窗口尺寸。

3) CCD影像3个波段的反演精度优于波段平均影像的反演精度，波段平均影像的自变量因子少且模型较前者更稳定。

FOTO方法具有入选变量少、不需利用其他多源数据 (如气象、地形数据和植被信息等) 协同反演且利用多元线性回归即能得到较高精度的反演结果等优点，随着我国高分辨率对地观测计划 (高分-2号影像的空间分辨率为0.8 m) 和国内外一些商业高分辨率遥感卫星数据 (高景一号分辨率为0.5 m、WorldView-3的为0.3 m) 的广泛获取，FOTO纹理指数对提高生物量的估算精度和效率具有十分重要的意义。