松材线虫病的高光谱综合分析早期探测方法

摘要

一种松材线虫病的高光谱综合分析早期探测方法。在松材线虫病发病季节，定期采取松树针叶样本用保鲜袋密封带回，在暗室内在人工光源下用野外光谱辐射仪测定样本的光谱数据，计算红边位置、绿峰反射高度、红光吸收深度、红绿指数、水绿指数、分形维数，继而综合分析这些数值，得出探测结论。本发明的有益效果是：用接近连续的包括可见光、近红外、部分中红外的高光谱数据，从细微特征到整体特征探测松树因患病而改变的光谱特性，微观与宏观兼顾；综合应用各种特征指数，充分考虑特征间的联系。其中红绿指数、水绿指数为本发明首次提出，它们分别反映了植物的光合作用和含水率情况。

说明书

【技术领域】

本发明属于遥感技术应用领域，具体是一种基于地面非成像高光谱数据综合分析的松材线虫病早期探测方法。

【背景技术】

松材线虫病(bursaphelenchus xylophilus)有松树“癌症”之喻，目前在我国的江苏、安徽、浙江、广东、湖北等南方大部分省都有分布。松材线虫病的主要寄主包括马尾松(P.massoniana)、云南松(P.yunnanensis)、思茅松(P.kesiya)、欧洲黑松(P.nigra)、火炬松(P.taeda)、湿地松(P.elliottii)等80种松属植物和雪松(Cedrus deodara)、香脂冷松(Abies balsamea)、黑云杉(Picea mariana)等27种非松属植物，而松褐天牛(Monochamus alternatus)是该病的主要传媒昆虫。该病自上个世纪八十年代传入我国以来，对松林及其森林生态环境造成了极大的危害，同时使一些重要风景区面临严重威胁。目前已有技术对发病初期的病害木进行抢救治疗，如挂瓶输液，因而对一些珍贵古树名木如黄山迎客松等，如能在发病初期发现，及时采取措施，仍有挽救的可能，因此，对珍贵树木来说，早期探测意义重大。

高光谱遥感在植物胁迫方面的研究具有很好的应用，可见光区的“蓝边”(蓝过渡到绿)、绿峰、“黄边”(绿过渡到红)、红光低谷及红光过渡到近红外的“红边”，是描述植物色素状态和健康状况的重要指示波段，实验证明，红边位置随叶绿素含量、叶面积指数、生物量、年龄、植物的健康程度、季节而变化，当植物健康叶绿素含量高时，红边位置往长波方向移动，当遭受病虫害时，红边位置向短波方向移动。因此目前多采用红边位置、绿峰高度、红谷吸收深度等参数分析植物的健康状况。如图2所示，反映植物健康状况的另一个重要因子是植物的含水率，中红外波段是研究植物含水率的重要波段。通常当植物得病后，蒸腾作用下降，含水率降低，反映在中红外波段上就是反射率增强。这些观点和反映出的现象已为业界学者所公认。

现有植物病虫害高光谱探测方面的技术和研究存在下列不足：一是仅使用少数几个孤立波段的数据，对植物患病前后光谱特性改变的研究不够全面；二是孤立应用一些参数或特征，对它们的相互联系考虑较少；三是对曲线形状的改变，缺乏一个整体衡量参数。

【发明内容】

针对现有技术的上述不足，本发明要解决的技术问题是提供一种松材线虫病的高光谱综合分析早期探测方法。

解决上述技术问题的技术方案是按如下步骤进行：

(1)松树针叶反射率连续测量：

a、待探测松树确定：侧重于具有较高价值、需要特别保护的珍贵松树作待探测对象，总数在10株以上；

b、测定期与测定时间的确定：地处北半球的我国选择每年4-10月松材线虫病发病季节，采样时间9:00-11:00，每3天采样1次；

c、样本采集：每棵树采一位于树木中上部且阳光能照射到的末端枝，放入保鲜袋密封，及时带回室内测定；

d、光谱测定：暗室内在人工光源下用便携式野外光谱辐射仪测定每个样本的高光谱数据，将仪器设定为一个样本10次重复；要求光谱辐射仪的物理光谱分辨率10nm以上，插值后输出光谱分辨率1nm，波长范围350-2000nm；

(2)光谱数据预处理：利用仪器提供的软件将测定的数据转换为普通可读写格式；求每个样本10次重复数据的算术平均数，作为该样本的光谱数据用于后面的计算；

(3)计算红边位置RE：红边位置是670nm-780nm波长范围内光谱反射率拐点处的波长，采用近似一阶导数法计算RE：

$R_{\max }^{\prime }\left(\lambda \right)\approx \underset{\lambda \in \left[\mathrm{670,780}\right]}{\max }\left(\frac{R\left({\lambda }_{i+1}\right)-R\left({\lambda }_{i-1}\right)}{{\lambda }_{i+1}-{\lambda }_{i-1}}\right)$

RE＝R′_max^-1(R′_max(λ))

式中λ、R(λ)、R′_max(λ)和i分别为波长、波长为λ的光谱反射率、光谱反射率一阶导数最大值和波长序号，RE＝R′_max^-1(R′_max(λ))表示R′_max(λ)的反函数，即光谱反射率一阶导数最大值处的波长；

(4)计算绿峰反射高度H_G：

$H_G=1-\frac{|R_A+\frac{R_C-R_A}{{\lambda }_C-{\lambda }_A}({\lambda }_B-{\lambda }_A)|}{R_B}$

式中R_A、R_B和R_C分别为吸收特征起点、中心点和结束点处的光谱反射率，λ_A、λ_B和λ_C分别为反射特征起点、中心点和结束点处的波长，其值分别为560nm、500nm和670nm；

(5)计算红光吸收深度H_R：

$H_R=1-\frac{R_C}{|R_B+\frac{R_D-R_B}{{\lambda }_D-{\lambda }_B}({\lambda }_C-{\lambda }_B)|}$

式中R_B、R_C和R_D分别为吸收特征起点、中心点和结束点处的光谱反射率，λ_B、λ_C和λ_D分别为反射特征起点、中心点和结束点处的波长，其值分别为670nm、560nm和760nm；

(6)计算红绿指数RG：先计算500nm-600nm波长范围内即绿光区域反射率最大值GR_max和对应的波长Gλ_max，600nm-700nm波长范围内即红光区域反射率最小值RR_min和对应的波长Rλ_min：

${\mathrm{GR}}_{\max }=\underset{\lambda \in \left[\mathrm{500,600}\right]}{\max }R\left(\lambda \right)$

Gλ_max＝R^-1(GR_max)

${\mathrm{RR}}_{\min }=\underset{\lambda \in \left[\mathrm{600,700}\right]}{\min }R\left(\lambda \right)$

Rλ_min＝R^-1(RR_min)

按下式计算红绿指数RG：

$\mathrm{RG}=\frac{{\mathrm{RR}}_{\min }-{\mathrm{GR}}_{\max }}{{\mathrm{R\lambda }}_{\min }-{\mathrm{G\lambda }}_{\max }}$

(7)计算水绿指数WG：先计算1360nm-1570nm和1830nm-1960nm两个波长范围内的反射率最小值WR_min1、WR_min2和对应的波长Wλ_min1、Wλ_min2：

${\mathrm{WR}}_{\min 1}=\underset{\lambda \in \left[\mathrm{1360,1570}\right]}{\min }R\left(\lambda \right)$

Wλ_min1＝R^-1(WR_min1)

${\mathrm{WR}}_{\min 2}=\underset{\lambda \in \left[\mathrm{1830,1960}\right]}{\min }R\left(\lambda \right)$

Wλ_min2＝R^-1(WR_min2)

按下式计算水绿指数WG：

$\mathrm{WG}=\frac{({\mathrm{WR}}_{\min 1}+{\mathrm{WR}}_{\min 2})/2-{\mathrm{GR}}_{\max }}{({\mathrm{W\lambda }}_{\min 1}+W{\lambda }_{\min 2})/2-G{\lambda }_{\max }}$

(8)计算分形维数fd：采用方差法计算350nm-2000nm波长范围内反射光谱曲线的分形维数：

a、尺度r设定：尺度范围{2，4，6，…，64}，取起始尺度r＝2；

b、用以r为宽的矩形框首尾相接将光谱曲线覆盖起来，如图4所示。设第i个框内光谱曲线的最大值与最小值之差为H_i，若尺度r很小，H_i值就逼近第i个框内曲线的长度，则等价的测度数即方差V(r)的表达式为

$V\left(r\right)=\frac{\mathrm{\Sigma r}{H}_i}{r^2}=\frac{\Sigma H_i}{r}$

c、令r＝r+2，若r＞64，按步骤(8)的d方法处理，否则按步骤(8)的b方法处理；

d、以ln(V(r))为因变量、ln(V(r))为自变量建立直线回归方程，由直线的斜率k得到光谱曲线的分形维数.fd＝2-k；

(9)病害早期诊断：据上述计算各数值，用纵横比较法进行诊断，纵向比较是将某树当前的数据和它以前的相应数据相比较，横向比较是将某树的数据和它树的同一时间相应数据相比较；

诊断依据如下：

a、若红边位置蓝移，即RE值变小，则树木有患病可能；

b、若绿峰反射高度降低，即H_G值变小，则树木有患病可能；

c、若红光吸收深度降低，即H_R值变大，则树木有患病可能；

d、若红绿指数RG增大，则树木有患病可能；

e、若水绿指数WG增大，则树木有患病可能；

f、若分形维数fd变小，则树木有患病可能；

综合上述6种情况，做出最后结论。

本发明具有下列有益效果：一是用接近连续的包括可见光、近红外、部分中红外的高光谱数据，从细微特征到整体特征探测松树因患病而改变的光谱特性，微观与宏观兼顾；二是综合应用红边位置、绿峰反射高度、红光吸收深度、红绿指数、水绿指数、分形维数等特征数据进行探测，充分考虑各种特征数据间的联系。上述一和二使得松材线虫病早期探测准确率较高；三是用分形维数刻画整体曲线特征，提供了一个整体衡量参数，简洁直观。

【附图说明】

图1为本发明的流程图。

图2为三种不同健康状况马尾松针叶反射光谱曲线图。

图3为典型植物可见光到近红外区域反射光谱曲线图。

图4为分形维数计算示意图。

图5为马尾松松材线虫病早期探测分形维数变化趋势模型图。

【具体实施方式】

下面结合实施例并参照附图对本发明作进一步详述：

需要说明的是，本发明以马尾松为例，将患松材线虫病简称为患病。本发明的流程如图1所示。具体步骤如下：

(1)松树针叶反射率连续测量：

a、待探测松树确定：侧重于具有较高价值、需要特别保护的珍贵松树作待探测对象，总数在10株以上；

b、测定期与测定时间的确定：地处北半球的我国选择每年4-10月松材线虫病发病季节，采样时间9:00-11:00，每3天采样1次；

c、样本采集：每棵树采一位于树木中上部且阳光能照射到的末端枝，放入保鲜袋密封，及时带回室内测定；

d、光谱测定：暗室内在人工光源下用便携式野外光谱辐射仪测定每个样本的高光谱数据，将仪器设定为一个样本10次重复；要求光谱辐射仪的物理光谱分辨率10nm以上，插值后输出光谱分辨率1nm，波长范围350-2000nm，本例采用美国制造的便携式ASD(Analyzed Spectral Devices，FieldSpec Full Range)野外光谱辐射仪，其基本指标：光谱范围：350～2500nm；光谱分辨率：3nm(350～1000nm)，10nm(1000～2500nm)；扫描时间：100ms；采样间隔：1.4nm(350～1050nm)，2nm(1000～2500nm)；波长精度：±1nm；标准参考板：聚四氟乙烯标准白板；前方探头为光纤探头；

(2)光谱数据预处理：利用仪器提供的软件将测定的数据转换为普通可读写格式；求每个样本10次重复数据的算术平均数，作为该样本的光谱数据用于后面的计算；

(3)计算红边位置RE：红边位置是670nm-780nm波长范围内光谱反射率拐点处的波长，采用近似一阶导数法计算RE：

$R_{\max }^{\prime }\left(\lambda \right)\approx \underset{\lambda \in \left[\mathrm{670,780}\right]}{\max }\left(\frac{R\left({\lambda }_{i+1}\right)-R\left({\lambda }_{i-1}\right)}{{\lambda }_{i+1}-{\lambda }_{i-1}}\right)$

RE＝R′_max^-1(R′_max(λ))

式中λ、R(λ)、R′_max(λ)和i分别为波长、波长为λ的光谱反射率、光谱反射率一阶导数最大值和波长序号，RE＝R′_max^-1(R′_max(λ))表示R′_max(λ)的反函数，即光谱反射率一阶导数最大值处的波长；

当松树发生病害后，红边位置蓝移，即往短波方向移动；

(4)计算绿峰反射高度H_G：

$H_G=1-\frac{|R_A+\frac{R_C-R_A}{{\lambda }_C-{\lambda }_A}({\lambda }_B-{\lambda }_A)|}{R_B}$

式中R_A、R_B和R_C分别为吸收特征起点、中心点和结束点处的光谱反射率，λ_A、λ_B和λ_C分别为反射特征起点、中心点和结束点处的波长，其值分别为560nm、500nm和670nm，见图3；

绿峰反射高度反映松树叶子的叶绿素含量状况，健康松树叶绿素含量高，绿峰反射高度较高，患病后叶绿素含量下降，绿峰反射高度降低；

(5)计算红光吸收深度H_R：

$H_R=1-\frac{R_C}{|R_B+\frac{R_D-R_B}{{\lambda }_D-{\lambda }_B}({\lambda }_C-{\lambda }_B)|}$

式中R_B、R_C和R_D分别为吸收特征起点、中心点和结束点处的光谱反射率，λ_B、λ_C和λ_D分别为反射特征起点、中心点和结束点处的波长，其值分别为670nm、560nm和760nm，见图3；

植物光合作用需要吸收红光波段电磁波，健康松树光合作用强，吸收红光波段电磁波较多，红光吸收深度较深，患病后光合作用变弱，红光吸收深度变浅；

(6)计算红绿指数RG：先计算500nm-600nm波长范围内即绿光区域反射率最大值GR_max和对应的波长Gλ_max，600nm-700nm波长范围内即红光区域反射率最小值RR_min和对应的波长Rλ_min：

${\mathrm{GR}}_{\max }=\underset{\lambda \in \left[\mathrm{500,600}\right]}{\max }R\left(\lambda \right)$

Gλ_max＝R^-1(GR_max)

${\mathrm{RR}}_{\min }=\underset{\lambda \in \left[\mathrm{600,700}\right]}{\min }R\left(\lambda \right)$

Rλ_min＝R^-1(RR_min)

按下式计算红绿指数RG：

$\mathrm{RG}=\frac{{\mathrm{RR}}_{\min }-{\mathrm{GR}}_{\max }}{{\mathrm{R\lambda }}_{\min }-{\mathrm{G\lambda }}_{\max }}$

患病后绿光区域反射率最大值GR_max降低，红光区域反射率最小值RR_min升高，导致RG增大；

(7)计算水绿指数WG：先计算1360nm-1570nm和1830nm-1960nm两个波长范围内的反射率最小值WR_min1、WR_min2和对应的波长Wλ_min1、Wλ_min2：

${\mathrm{WR}}_{\min 1}=\underset{\lambda \in \left[\mathrm{1360,1570}\right]}{\min }R\left(\lambda \right)$

Wλ_min1＝R^-1(WR_min1)

${\mathrm{WR}}_{\min 2}=\underset{\lambda \in \left[\mathrm{1830,1960}\right]}{\min }R\left(\lambda \right)$

Wλ_min2＝R^-1(WR_min2)

按下式计算水绿指数WG：

$\mathrm{WG}=\frac{({\mathrm{WR}}_{\min 1}+{\mathrm{WR}}_{\min 2})/2-{\mathrm{GR}}_{\max }}{({\mathrm{W\lambda }}_{\min 1}+W{\lambda }_{\min 2})/2-G{\lambda }_{\max }}$

1440nm附近和1880nm附近是两个重要的水吸收带，患病后松树蒸腾作用下降，含水率降低，这两个吸收带反射率增强，WG增大；

(8)计算分形维数fd：采用方差法计算350nm-2000nm波长范围内反射光谱曲线的分形维数：

a、尺度r设定：尺度范围{2，4，6，…，64}，取起始尺度r＝2；

b、用以r为宽的矩形框首尾相接将光谱曲线覆盖起来，如图4所示。设第i个框内光谱曲线的最大值与最小值之差为H_i，若尺度r很小，H_i值就逼近第i个框内曲线的长度，则等价的测度数即方差V(r)的表达式为

$V\left(r\right)=\frac{\mathrm{\Sigma r}{H}_i}{r^2}=\frac{\Sigma H_i}{r}$

c、令r＝r+2，若r＞64，按步骤(8)的d方法处理，否则按步骤(8)的b方法处理；

d、以ln(V(r))为因变量、ln(V(r))为自变量建立直线回归方程，由直线的斜率k得到光谱曲线的分形维数fd＝2-k；

健康松树其绿峰、红谷、近红外高反射区、中红外水吸收区特征明显，光谱曲线形状复杂，分形维数高，患病后绿峰降低、红谷抬高、近红外高反射区降低、中红外水吸收区抬高，光谱曲线形状复杂度降低，分形维数下降。以时间为横坐标，分形维数为纵坐标绘制曲线，若分形维数持续降低，如图5，表示松树开始病变；若与同时间测定的其它树木比较，分形维数明显较低，也表示该树有患病可能；

(9)病害早期诊断：据上述计算各数值，用纵横比较法进行诊断。纵向比较是将某树当前的数据和它以前的相应数据相比较，横向比较是将某树的数据和它树的同一时间相应数据相比较。

诊断依据如下：

a、若红边位置蓝移，即RE值变小，则树木有患病可能；

b、若绿峰反射高度降低，即H_G值变小，则树木有患病可能；

c、若红光吸收深度降低，即H_R值变大，则树木有患病可能；

d、若红绿指数RG增大，则树木有患病可能；

e、若水绿指数WG增大，则树木有患病可能；

f、若分形维数fd变小，则树木有患病可能；

综合上述6种情况，做出最后结论。

步骤(6)和步骤(7)所述有关红绿指数、水绿指数的概念、计算方法为本发明首次提出。