中国森林病虫害共有8 000多种，经常造成危害的有200多种，森林病虫害大面积发生，危害极为严重，每年发生面积800万hm²左右，损失木材生长量1 700万m³，造成经济损失50亿元以上(武红敢等，2004)。松材线虫病(Bursaphelenchus xylophilus)是松树的一种毁灭性病害，该病也称作松材线虫枯萎病和松材线虫萎蔫病，于20世纪80年代传进我国(王琼瑜，1992)，目前主要分布在我国的江苏、安徽、山东、浙江、广东、湖北(恩施)、湖南、香港等地以及日本、朝鲜、韩国、美国、加拿大、墨西哥、葡萄牙等国家和地区。松材线虫病的主要寄主包括日本赤松(Pinus densiflora)、云南松(P. yunnanensis)、思茅松(P. kesiya)、马尾松(P. massoniana)、欧洲黑松(P. nigra)、火炬松(P. taeda)、湿地松(P. elliottii)、琉球松(P. luchuensis)、白皮松(P. bungeana)等80种松属植物和雪松(Cedrus deodara)、香脂冷松(Abies balsamea)、日本落叶松(Larix leptolepis)、黑云杉(Picea mariana)等27种非松属植物，而松褐天牛(Monochamus alternatus)是该病的主要传媒昆虫(王明旭，2007)。

监测是松材线虫治理的关键环节，特别是松材线虫病爆发前的早期预测预警。王明旭(2007)从管理系统、监测系统、生态恢复系统等方面论述了湖南松材线虫病早期预测预警系统的构建和实践；石进等(2006)采用航空遥感结合地面GPS定位相结合的方法对有松材线虫病症状的松林进行准确、快速定位；王震等(2007)对马尾松4种不同受害类型的林木进行反射光谱测量，表明不同受害类型马尾松光谱曲线差异明显，在绿光区，随受害程度加深叶绿素和水分减少，其反射率逐渐减小；而在红光区，随受害加深，其反射率逐渐增大，为研究松材线虫入侵过程及相应的监测管理提供了基础。

高光谱遥感在研究植物胁迫方面的广泛应用，可见光区的“蓝边”(蓝过渡到绿)、绿峰、“黄边”(绿过渡到红)、红光低谷及红光过渡到近红外的“红边”是描述植被色素状态和健康状况的重要指示波段(张金恒等，2003)。“红边”(680~780 nm)是植物光谱曲线最明显的特征，试验和理论都证明，红边位置REP(red edge position)随叶绿素含量、叶面积指数(LAI)、生物量(biomass)、年龄、植被的健康程度、季节而变化，当植被健康并有高叶绿素含量和高LAI，红边位置将往长波方向移动，当遭受病虫害或萎黄病发病时，LAI低的时候，红边位置将向短波方向移动(Helmi et al., 2006)。因此，REP成为植被胁迫和衰老的指示器，常用它用来研究植物的生长和健康状况(Jago et al., 1995; Dawson et al., 1998; Boochs et al., 2003)。

松材线虫病发病早期，病害木自身无外观症状，其早期发病的探测、诊断方面的研究还比较少。已有的研究也仅对比已知不同受害程度马尾松反射光谱的变化情况(王震等，2007)。在未知情况下，反射光谱曲线及特征参数等能否对发病早期植株进行探测值得研究，而且对松材线虫的防治也更具实践意义。

分形维数测量可以理解空间数据的复杂性，其构形、格局，尤其是距离、曲率和形状的测量能够提供相关结构与过程的潜在信息(Lam et al., 1993)。周子勇等(2005)指出，高光谱遥感光谱曲线具有分形特征，这为高感光遥感分形测量提供了基础；贝淑坤(1997)计算了不同药品红外光谱盒维数，简述了分形维数与分子结构的关系，指出光谱波动越剧烈，分形维数越高；张平等(2007)探讨了分形理论在太赫兹光谱识别中的应用，计算了光谱曲线非常相似的红参(Radix ginseng)和高丽参光谱曲线盒维数，结果表明盒维数能很好的将2者区分开。本研究在计算发病早期马尾松高光谱数据REP等特征参数的基础上，采用分形理论来分析发病早期马尾松高光谱数据，探讨其对马尾松材线虫病发病早期探测预测的可行性，为马尾松林病虫害治理、生态系统恢复等提供基础。

1 研究区域及数据获取 1.1 研究区域与测量仪器

研究地区为浙江省宁波市象山县四洲头镇(图 1)。共设3块样地进行调查，分别为1号样地：121°45′20.7″ E，29°22′58.4″ N；2号样地：121°44′49.6″ E，29°23′3.5″ N；3号样地：121°45′23.7″ E，29°22′51.8″ N。

野外光谱测量采用浙江林学院资源监测与信息管理实验室的美国制造的便携式ASD (analyzed spectral devices), FieldSpec Full Range野外光谱测量仪。其基本指标为：光谱范围：350~2 500 nm；光谱分辨率：3 nm (350~1 000 nm)，10 nm (1 000~2 500 nm)；扫描时间：100 ms；采样间隔：1.4 nm(350~1 050 nm)，2 nm(1 000~2 500 nm)；波长精度：±1 nm；标准参考板：聚四氟乙烯标准白板；前方探头为光纤探头。

1.2 试验设计与测量方法

所选3个样地的马尾松测量之前均处于健康状态。2007年7月17日—8月22日，每3天对样地的每棵马尾松进行反射光谱测量(遇阴雨，则停测)，并记录健康状况。如果发现某棵树部分针叶开始发黄或发现枯枝，由于测死枝的结果已失去探测预测的价值仍然测该树的活枝。

为减少太阳辐射影响，测量时间选在11: 00—13: 00，样木活枝针叶测量在1 min内完成(参考板和目标)。测量时将待测样本平铺于黑布之上(黑布反射率接近0，可以抑制地面背景的影响)，用三角架固定光谱仪探头垂直向下测量。每次测量之前先进行参考板测量和自动优化，然后设置仪器自动获取同一目标10组光谱数据，并取平均值做为这棵树的反射光谱值。另外，由于H₂O和CO₂ 对太阳能的吸收，野外测量光谱在1 400和1 900 nm附近的反射率都是随机值，需要做简单的剔除处理。

2 研究方法 2.1 高光谱参数

本研究主要采用红边位置、绿峰的反射高度、红谷的吸收深度等3个参数描述马尾松高光谱特征。红边通常指680~780 nm间反射率突然变化的区域，它是由于叶绿素的强吸收和近红外叶子内部散射导致的高反射所形成的(Horler et al., 1980；1983; Miller et al., 1990; Dawson et al., 1998; Helmi et al., 2006; Moses et al., 2006), 如图 2 C, D之间的曲线。红边位置指红边范围光谱反射率的拐点，或红边区光谱斜率最大的地方(Dawson et al., 1998; Moses et al., 2006)，目前关于红边位置的计算方法很多，包括一阶导数法(吴彤等，2007)、4点线性插值法(Helmi et al., 2006; Dawson et al., 1998; Moses et al., 2006)、基于导数光谱的拉格朗日3点插值法(Dawson et al., 1998)、倒高斯模型(吴继友等，1997)、多项式拟合法(Moses et al., 2006; Pu et al., 2003)以及线性外推法(Moses et al., 2006)等。鉴于导数光谱分析技术对光照强度、噪声影响敏感性较低，同时能增强光谱特征等(浦瑞良等，2000；Helmi et al., 2006)，采用导数光谱分析技术REP，即一阶导数最大值所对应的波长来计算红边位置，导数光谱的技术方法如下(Dawson et al., 1998; 范文义等，2004; 张风丽等，2005):

(1)

式中R_i, λ_i分别为第i波段的反射率和波长。

绿峰反射高度H_G、红谷吸收深度H_R分别按公式(2)，(3)计算(吴彤等，2007)，式中R_A, R_B, R_C, R_D分别代表图 2中A，B，C，D处的反射率；λ_A, λ_B, λ_C, λ_D为相应波长。

(2)

(3)

2.2 反射光谱曲线分形测量

曲线分形维数的计算方法很多，如盒子法、方差法、结构函数法、协方差加权、量规法等(Lam et al., 1993; 葛世荣等，1997; Zhu et al., 2000)。反射光谱曲线局部变异会对整个曲线带来影响，故选择方差法计算反射光谱曲线的分形维数(葛世荣等，1997)：

首先以宽为r的矩形框首尾相接将表面轮廓曲线覆盖起来，令第i个框内表面轮廓的最大值与最小值之差为H_i，若尺度r很小，则H_i值就逼近第i个框内曲线的长度。则等价的测度数V(r)的表达式为:

(4)

然后改变尺度r，并按(4)计算该尺度下的方差；如此继续，最后V(r)与r在双对数坐标进行最小二乘回归计算，则由斜率α便可以得到曲线的分形维数D如下：

(5)

3 结果与分析 3.1 光谱特征参数分析

3个样地不同样木红边位置、绿峰反射高度和红谷吸收深度的变化曲线如图 3。

从8月6日开始1号样地的10和12号样木红边位置降低，5和19号样木从8月12开始，红边也有下降趋势；2号样地的4棵样木的红边位置从8月7号开始较大的波动，2号和12号样木红边位置大幅降低，但12日后又升高，11和21号样木在8月12日到8月17日之间呈下降趋势，8月17日4棵样木的红边位置相对初始时间(7月18日)均有一定程度的降低；3号样地的4棵样木从8月6日开始，红边位置虽稍有改动，但均比初始时间低。

根据测量和调查记录，1和2号样地内的样木在测量时间内，肉眼未发现有病症状；3号样地的1和29号样木在8月初出现病症。后期补充调查发现：1)3号样地的1和29号样木分别于8月31日和9月8日死亡；2)10月13日，1号样地的12和15号样木，2号样地的4，15，18号样木，以及3号样地的14，28，30号样木均死亡。

红边位置降低(蓝移)说明植被健康状况开始恶化，根据以上分析，3号样地的1，29号样木与此相符，且从发现病症到死亡仅经历1个月左右的时间。1号样地的12号样木，3号样地的14，28号样木，在测量期间没有出现病症的情况下于测量后2个月死亡，其红边位置同样有下降趋势，但存在波动。树木出现病虫害时，其活枝还能在一段时间内进行正常的光合作用，因此，波动不可避免；另外，就3号样地的4棵样木和1号样地的12号样木来看，死亡经历的时间越短，红边位置降低幅度越大，且波动越小。

植被受损后，其绿峰高度和红谷吸收深度均会减小(吴彤等，2007)。1，2号样地绿峰反射高度从8月12日开始，基本呈下降趋势，3号样地的波动性较大，但8月12日后，其绿峰反射高度也基本低于初始值。

红谷吸收深度方面，1号样地的10号样木8月12日后呈下降趋势，2号样木8月13日后开始降低，3号样地8月4日后开始下降，但也存在起伏，如29号样木。

在未知是否会发生松材线虫病的情况下，通过红边位置、绿峰反射高度以及红谷吸收深度变化分析，以及后期补充调查，这3个参数特别是红边位置能对病虫害的发生起到一定的预测作用。但是也出现不同程度的波动，这种波动可能与病害程度以及出现症状到死亡所经历的时间有关。

3.2 光谱曲线分形分析

红边位置、绿峰反射高度、红谷吸收深度主要在450~780 nm光谱范围间提取，为保持一致性，光谱分形维数也在这个范围内按公式(4)、(5)计算，结果如图 4。

1，2号样地分形维数整体上具有逐渐降低的趋势，特别是在8月12日以后。3号样地分形维数在8月6日以后基本呈下降趋势。

3.2.1 分形维数意义及曲线构形变化机理

分形测量可以解释空间数据的复杂性、形状的变化等(Lam et al., 1993)，这种变化可以通过分形维数定量分析。分形维数是所有分形对象的主要特征之一，它量度分形物体的不规则和断裂程度(文志英等，1999)，同时分形维数也能定量的反映某一构形的复杂性及填充空间的能力，分维数越大，构形越复杂，其填充空间的能力越强；反之，分维数越小，构形越简单，其填充空间的能力越小。当植被受病虫害影响时，红边蓝移、红谷抬高、绿峰降低，随着受害程度的加剧，可见光到近红外反射光谱曲线的波峰、波谷特征逐渐消失，逐渐变平、拉直(梅安新等，2001；季荣等，2003)，因此，光谱曲线本身的几何构形、形状、曲率等发生变化，曲线变的越来越简单，而近似一条直线，分形维数也会随着曲线的简化而降低。

为了说明曲线复杂程度变化与分形维数变化之间的关系，进一步从理论上模拟了具有不同复杂程度的分形布朗曲线如图 5，其中h为hurst指数，它是决定曲线粗糙程度的指标。曲线的分形维数D=2-h，当h=1/2时；曲线最不稳定，当h>1/2时曲线趋于简单；当h＜1/2时，曲线趋于复杂，有较大的分形维数值(戴前伟等，1998；杨山等，2002)。从图 5可以看出，当分形维数为1.9时，曲线很粗糙其且复杂，而当分形维数为1.25时，曲线光滑且简单。

由图 4可以看出，3号样地4棵样木分形后期分形维数呈下降趋势，其变化与红边位置后期的变化基本吻合，但分形维数变化的波动明显比红边位置小，稳定性好。如果3号样地分析是正确的，那么2和3号样地分形维数的下降就能对松材线虫并的发生起到预测作用：如1号样地第12号样木后期死亡；另外，尽管2号样地4棵样木在后期的观测没有发现病情，但相邻的4，15，18号样木却死亡，是否也有警示作用，毛者将在下面进行分析验证。

3.2.2 分形维数与光谱特征参数之间的关系

结合图 2及上面的分析，在450~780nm之间，红谷吸收深度和绿峰反射高度的变化是影响曲线构形的主要因素；而红边变化主要是倾斜程度和红边位置左右移动，其直线构形基本不发生变化。因此，笔者分析了绿峰高度和红谷吸收深度与分形维数之间的相关关系(表 1)。

除1号样地的第19棵样木外，分形维数与绿峰反射高度均呈正相关关系。其中3号样地的4棵树中，1，29 2棵样木的相关系数明显高于14和28号样本；结合先前分析，这种相关系数的变化可能与病害程度有关。另外，1号样地的第12棵树的相关系数也较高。

分形维数与红谷吸收深度相关关系方面，2号样地的4棵树基本呈负相关或不相关；而与1，2号样地基本成正相关关系(除1号样地的第19和2号样地的第28号样木)。

另外，分形维数与绿峰反射高度的相关系数高于与红谷吸收深度之间的关系，就本研究而言，为说明绿峰对马尾松病害预测的效果更好，且绿峰变化对曲线分形维数的变化更敏感，笔者进一步将同一样地与分形维数呈正相关样木的所有观测值作为整体进行分析。图 6是3个样地绿峰高度、红吸收深度与分形维数之间的关系模型(因2号样地红谷吸收深度与分形基本呈负相关或不相关，故没有对其进行分析)。从图 6可以更直观的看出分形维数与绿峰反射高度的相关系数同样高于与红谷吸收深度之间的关系。这些关系也说明：如果用绿峰反射高度和红谷吸收深度指示植被的健康状况，结合分形维数意义及曲线构形变化分形机理，分形维数的变化同样适合于对植被健康状况的预测。

3.2.3 补充调查死亡树木的验证分析

鉴于1号样地的12，15号样木，2号样地的4，15，18号样木，以及3号样地的14，28，30号样木在后期补充调查时均死亡。我们按上述方法计算这5棵树反射光谱曲线的分形维数，分析发现它们的分析维数在后期也呈下降趋势；在此基础上进一步建立了它们与绿峰反射高度之间的关系如图 7。可以看出，分形维数与绿峰反射高度具有线性关系，这进一步验证了上面的分析，也再次说明植被反射光谱曲线分形分析对其健康状况具有一定的预测能力。

分形分析是一种多尺度的分析方法，能够用来分析数据的各种空间现象，描述现实世界详细的等级层次(Lam et al., 1993; Migue et al., 2005)。随着尺度r的变化，曲线的细节可以在不同的层次上得以表现。因此，相对与绿峰反射高度、红谷吸收深度等常规光谱参数而言，分形维数是对曲线整体变化的一个综合分析和定量描述，因健康状况引起的光谱曲线变化是曲线本身的几何构形、形状、曲率等的综合变化，而不是局限于某个或某几个特征的参数的变化。

4 结论与讨论

1) 对于病害木，如3号样地，其活枝反射光谱红边位置在后期有明显降低趋势，当然也存在一些起伏，但起伏的大小可能与发生病症到死亡所经历的时间长短有一定关系。

2) 相对于初始观测结果，在观测后期，绿峰反射高度和红谷吸收深度有不同程度的降低。结合红边位置，这3个参数对松材线虫的发生具有一定预测作用，这不同于在已知病害程度情况下，马尾松反射光谱特征参数的显著变化(王震等，2007)，因为，出现症状后，其活枝仍能进行光合作用。

3) 分形分析表明：3个样地随机选择的样木到观测后期，分形维数也呈降低趋势。结合分形维数的意义、反射光谱曲线变化的机理以及不同复杂程度的分形布朗运动曲线的模拟可知，在450~780 nm波长范围内反射光谱的曲线分形维数的变化是植被健康状况发生变化的综合描述，对马尾松松材线虫病害的发生具有一定的指示作用。

4) 分形维数与绿峰反射高度和红谷吸收深度存在正相关关系，且与绿峰的相关关系高于红谷吸收深度，是否说明绿峰的变化对马尾松健康状况的预测更有效，需要进一步研究。另外，对比光谱特征参数分析和后期调查，相关系数的大小与发现病害到死亡所经历的时间有一定的关系，经历时间越短，相关关系越大，反之亦然。这种相关不是简单的线性关系，因为分形维数是对非线性现象的综合定量描述。

5) 补充调查死亡样木反射光谱曲线的分形分析及与绿峰反射高度的关系模型，进一步验证分形维数对马尾松健康状况的分析是可行的。

6) 对于2号样地随机抽取的4棵树，其分形维数有下降趋势，但分形维数与红谷吸收深度呈负相关或不相关，这种相反的关系是否说明对于健康马尾松，红谷吸收深度对健康状况的变化不敏感；而分形维数与绿峰高度呈正相关，且在整体上与绿峰高度系数大于红谷，是否说明绿峰高度对马尾松松材线虫病的预测更有效，需要深入研究。

7) 马尾松发生松材线虫症状后，其活枝依然会进行光合作用，因此反射光谱会出现各种波动，即使健康植被，波动也不可避免，这就增加了通过高光谱对马尾松松材线虫症进行早期探测的难度。增加生物物理参数如叶绿素含量、含水量等可能会改进预测结果；另外，在双对数坐标通过V(r)与r绘图计算分形维数时，无标度区的选择会影响计算结果，有时V(r)与r之间可能出现分段线性的关系，说明存在多重分形特征(周子勇等，2005)，多重分形分析是否能够改善单分形分析结果，需要进一步分析。

Carter等(1994)通过胁迫方法研究不同植被光谱信息变化，指出在肉眼发现症状的16天前其光谱特征就有变化。本次研究表明：对没有施加任何胁迫或未知病害及程度的情况下，高光谱反射光谱及其特征参数对病害的预测还存在不确定性因素，如活枝反射光谱红边位置起伏波动等，这方面将结合生化参数如叶绿素、含水量等，进行深入研究。分形分析是非线性分析的新方法和思路，高光谱曲线具有分形特征(周子勇等，2005)。所以，植被光谱曲线分形分析是可行的。当植被受病虫害并随着受害程度加剧时，曲线变异机理与分形维数的意义一致。因此，分形维数在病虫害早期发生的预测中具有一定潜力，相关的研究还有待进一步深入。