马尾松反射光谱曲线分形特征研究

马尾松反射光谱曲线分形特征研究

中国有8000多种森林病虫害，其中200种容易造成危险。森林病虫害大而复杂。每年发生约800万件hm，损失木材18770万件hm，经济损失超过50亿元（吴红喷泉等，2004）。松材线虫病(Bursaphelenchusxylophilus)是松树的一种毁灭性病害,该病也称作松材线虫枯萎病和松材线虫萎蔫病,于20世纪80年代传进我国(王琼瑜,1992),目前主要分布在我国的江苏、安徽、山东、浙江、广东、湖北(恩施)、湖南、香港等地以及日本、朝鲜、韩国、美国、加拿大、墨西哥、葡萄牙等国家和地区。松材线虫病的主要寄主包括日本赤松(Pinusdensiflora)、云南松(P.yunnanensis)、思茅松(P.kesiya)、马尾松(P.massoniana)、欧洲黑松(P.nigra)、火炬松(P.taeda)、湿地松(P.elliottii)、琉球松(P.luchuensis)、白皮松(P.bungeana)等80种松属植物和雪松(Cedrusdeodara)、香脂冷松(Abiesbalsamea)、日本落叶松(Larixleptolepis)、黑云杉(Piceamariana)等27种非松属植物,而松褐天牛(Monochamusalternatus)是该病的主要传媒昆虫(王明旭,2007)。监测是松材线虫治理的关键环节,特别是松材线虫病爆发前的早期预测预警。王明旭(2007)从管理系统、监测系统、生态恢复系统等方面论述了湖南松材线虫病早期预测预警系统的构建和实践;石进等(2006)采用航空遥感结合地面GPS定位相结合的方法对有松材线虫病症状的松林进行准确、快速定位;王震等(2007)对马尾松4种不同受害类型的林木进行反射光谱测量,表明不同受害类型马尾松光谱曲线差异明显,在绿光区,随受害程度加深叶绿素和水分减少,其反射率逐渐减小;而在红光区,随受害加深,其反射率逐渐增大,为研究松材线虫入侵过程及相应的监测管理提供了基础。高光谱遥感在研究植物胁迫方面的广泛应用,可见光区的“蓝边”(蓝过渡到绿)、绿峰、“黄边”(绿过渡到红)、红光低谷及红光过渡到近红外的“红边”是描述植被色素状态和健康状况的重要指示波段(张金恒等,2003)。“红边”(680～780nm)是植物光谱曲线最明显的特征,试验和理论都证明,红边位置REP(rededgeposition)随叶绿素含量、叶面积指数(LAI)、生物量(biomass)、年龄、植被的健康程度、季节而变化,当植被健康并有高叶绿素含量和高LAI,红边位置将往长波方向移动,当遭受病虫害或萎黄病发病时,LAI低的时候,红边位置将向短波方向移动(Helmietal.,2006)。因此,REP成为植被胁迫和衰老的指示器,常用它用来研究植物的生长和健康状况(Jagoetal.,1995;Dawsonetal.,1998;Boochsetal.,2003)。松材线虫病发病早期,病害木自身无外观症状,其早期发病的探测、诊断方面的研究还比较少。已有的研究也仅对比已知不同受害程度马尾松反射光谱的变化情况(王震等,2007)。在未知情况下,反射光谱曲线及特征参数等能否对发病早期植株进行探测值得研究,而且对松材线虫的防治也更具实践意义。分形维数测量可以理解空间数据的复杂性,其构形、格局,尤其是距离、曲率和形状的测量能够提供相关结构与过程的潜在信息(Lametal.,1993)。周子勇等(2005)指出,高光谱遥感光谱曲线具有分形特征,这为高感光遥感分形测量提供了基础;贝淑坤(1997)计算了不同药品红外光谱盒维数,简述了分形维数与分子结构的关系,指出光谱波动越剧烈,分形维数越高;张平等(2007)探讨了分形理论在太赫兹光谱识别中的应用,计算了光谱曲线非常相似的红参(Radixginseng)和高丽参光谱曲线盒维数,结果表明盒维数能很好的将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(analyzedspectraldevices),FieldSpecFullRange野外光谱测量仪。其基本指标为:光谱范围:350～2500nm;光谱分辨率:3nm(350～1000nm),10nm(1000～2500nm);扫描时间:100ms;采样间隔:1.4nm(350～1050nm),2nm(1000～2500nm);波长精度:±1nm;标准参考板:聚四氟乙烯标准白板;前方探头为光纤探头。1.2活枝激发对比测量所选3个样地的马尾松测量之前均处于健康状态。2007年7月17日—8月22日,每3天对样地的每棵马尾松进行反射光谱测量(遇阴雨,则停测),并记录健康状况。如果发现某棵树部分针叶开始发黄或发现枯枝,由于测死枝的结果已失去探测预测的价值仍然测该树的活枝。为减少太阳辐射影响,测量时间选在11:00—13:00,样木活枝针叶测量在1min内完成(参考板和目标)。测量时将待测样本平铺于黑布之上(黑布反射率接近0,可以抑制地面背景的影响),用三角架固定光谱仪探头垂直向下测量。每次测量之前先进行参考板测量和自动优化,然后设置仪器自动获取同一目标10组光谱数据,并取平均值做为这棵树的反射光谱值。另外,由于H2O和CO2对太阳能的吸收,野外测量光谱在1400和1900nm附近的反射率都是随机值,需要做简单的剔除处理。2学习方法2.1导数光谱分析技术本研究主要采用红边位置、绿峰的反射高度、红谷的吸收深度等3个参数描述马尾松高光谱特征。红边通常指680～780nm间反射率突然变化的区域,它是由于叶绿素的强吸收和近红外叶子内部散射导致的高反射所形成的(Horleretal.,1980;1983;Milleretal.,1990;Dawsonetal.,1998;Helmietal.,2006;Mosesetal.,2006),如图2C,D之间的曲线。红边位置指红边范围光谱反射率的拐点,或红边区光谱斜率最大的地方(Dawsonetal.,1998;Mosesetal.,2006),目前关于红边位置的计算方法很多,包括一阶导数法(吴彤等,2007)、4点线性插值法(Helmietal.,2006;Dawsonetal.,1998;Mosesetal.,2006)、基于导数光谱的拉格朗日3点插值法(Dawsonetal.,1998)、倒高斯模型(吴继友等,1997)、多项式拟合法(Mosesetal.,2006;Puetal.,2003)以及线性外推法(Mosesetal.,2006)等。鉴于导数光谱分析技术对光照强度、噪声影响敏感性较低,同时能增强光谱特征等(浦瑞良等,2000;Helmietal.,2006),采用导数光谱分析技术REP,即一阶导数最大值所对应的波长来计算红边位置,导数光谱的技术方法如下(Dawsonetal.,1998;范文义等,2004;张风丽等,2005):dRdλ=R(λi+1)-R(λi-1)λi+1-λi-1‚(1)dRdλ=R(λi+1)−R(λi−1)λi+1−λi−1‚(1)式中Ri,λi分别为第i波段的反射率和波长。绿峰反射高度HG、红谷吸收深度HR分别按公式(2),(3)计算(吴彤等,2007),式中RA,RB,RC,RD分别代表图2中A,B,C,D处的反射率;λA,λB,λC,λD为相应波长。ΗG=1-|RA+RC-RAλC-λA(λB-λA)|RB‚(2)ΗR=1-RC|RB+RD-RBλD-λB(λC-λB)|。(3)2.2反射光谱分形维数法曲线分形维数的计算方法很多,如盒子法、方差法、结构函数法、协方差加权、量规法等(Lametal.,1993;葛世荣等,1997;Zhuetal.,2000)。反射光谱曲线局部变异会对整个曲线带来影响,故选择方差法计算反射光谱曲线的分形维数(葛世荣等,1997):首先以宽为r的矩形框首尾相接将表面轮廓曲线覆盖起来,令第i个框内表面轮廓的最大值与最小值之差为Hi,若尺度r很小,则Hi值就逼近第i个框内曲线的长度。则等价的测度数V(r)的表达式为:V(r)=∑rΗir2=∑Ηir‚(4)然后改变尺度r,并按(4)计算该尺度下的方差;如此继续,最后V(r)与r在双对数坐标进行最小二乘回归计算,则由斜率α便可以得到曲线的分形维数D如下:D=2-α。(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～780nm光谱范围间提取,为保持一致性,光谱分形维数也在这个范围内按公式(4)、(5)计算,结果如图4。1,2号样地分形维数整体上具有逐渐降低的趋势,特别是在8月12日以后。3号样地分形维数在8月6日以后基本呈下降趋势。3.2.1光谱分形维数变化分形测量可以解释空间数据的复杂性、形状的变化等(Lametal.,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,292棵样木的相关系数明显高于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。可以看出,分形维数与绿峰反射高度具有线性关系,这进一步验证了上面的分析,也再次说明植被反射光谱曲线分形分析对其健康状况具有一定的预测能力。分形分析是一种多尺度的分析方法,能够用来分析数据的各种空间现象,描述现实世界详细的等级层次(Lametal.,1993;Migueetal.,2005)。随着尺度r的变化,曲线的细节可以在不同的层次上得以表现。因此,相对与绿峰反射高度、红谷吸收深度等常规光谱参数而言,分形维数是对曲线整体变化的一个综合分析和定量描述,因健康状况引起的光谱曲线变化是曲线本身的几何构形、形状、曲率等的综合变化,而不是局限于某个或某几个特征的参数的变化。4研究问题及方法1)对于病害木,如3号样地,其活枝反射光谱红边位置在后期有明显降低趋势,当然也存在一些起伏,但起伏的大小可能与发生病症到死亡所经历的时间长短有一定关系。2)相对于初始观测结果,在观测后期,绿峰反射高度和红谷吸收深度有不同程度的降低。结合红边位置,这3个参数对松材线虫的发生具有一定预测作用,这不同于在已知病害程度情况下,马尾松反射光谱特征参数的显著变化(王震等,2007),因为,出现症状后,其活枝仍能进行光合作用。3)分形分析表明:3个样地随机选择的样木到观测后期,分形维数也呈降低趋势。结合分形维数的意义、反射光谱曲线变化的机理以及不同复杂程度的分形布朗运动曲线的模拟可知,在450～780nm波长范围内反射光谱的曲线分形维数的变化是植被健康状况发生变化的综合描述,对马尾松松材线虫病害的发生具有一定的指示作用。4)分形维数与绿峰反射高度和红谷吸收深度存在正相关关系,且与绿峰的相关关系高于红谷吸收深度,是否说明绿峰的变化对马尾松健康状况的预测更有效,需要进一步研究。另外,对比光谱特征参数分析和后期调查,相关系数的大小与发现