CN106708782B - 基于小波分析的区域病虫害检测诊断判别方法 - Google Patents

Abstract

本发明涉及数据遥感分析技术领域，是一种基于小波分析的区域病虫害检测诊断判别方法，包括以下步骤：第一步对获取胡杨虫害的试验数据进行抽取采样分析；第二步采用小波变换对采样数据进行初步分析；第三步定义并构建NDVI时间波谱曲线；第四步对NDVI时间波谱曲线进行滤噪处理；第五步对胡杨虫害进行增强和分离分析；第六步建立胡杨食叶害虫的信号检测模型；第七步对检测出的虫害信号进行实时动态检测。本发明基于森林病虫害遥感监测的原理与物理基础，借助NDVI时间序列数据，定义并构建含有虫害信息的胡杨NDVI时间波谱曲线，得到NDVI随时间变化的演变特性与规律，有效的对植物长势以及病虫害进行监测，避免胡杨长势逐渐衰弱。

Description

基于小波分析的区域病虫害检测诊断判别方法

技术领域

本发明涉及林业病虫害数据遥感分析技术领域，是一种基于小波分析的区域病虫害检测诊断判别方法。

背景技术

胡杨作为一种在河岸、湖泊等湿润条件下生长的荒漠河岸林，其独特的生理结构使其能够在盐碱、风沙和干旱的恶劣环境中生存，是干旱荒漠区近顶级的天然乔木群落，是荒漠生态系统的重要组成部分，在维护荒漠区生态平衡、防风固沙、调节气候、改善生态环境等方面发挥着重要的生态功能。但是近年来，由于人为过度引水和气候变化，导致胡杨生存环境恶化，长势逐渐衰弱，其正常的生理活动减弱，抵御病虫害的能力下降，使得病虫害蔓延迅速，开始大面积肆掠，其中尤以尺蠖虫害最为严重，严重时森林失叶如同遭遇火灾，而且有进一步向绿洲扩散的趋势。因此，研究胡杨尺蠖虫害特征及其分布规律对胡杨林的保护具有十分重要的意义。

在尺蠖研究方面，许多学者从尺蠖的生存环境、生理机制、生活习性、生活史规律、分布范围、种群动态变化、造成的危害状况、预测防治等方面已经做了大量研究。但是胡杨林区地处荒漠、交通不便，运用传统的实地调查方法难以对其进行大范围、深入的研究。遥感技术具有受地面条件限制少、获取资料速度快、周期短、用途广泛的优势，目前已成为国际上监测森林病虫危害最先进的手段之一，其在森林病虫害探测领域的应用已被证明是可行并且有效的，但目前这一技术尚处于探索发展阶段，仍然存在着遥感监测参数少，数据源及辅助信息的利用不足，技术与方法缺乏，针对性差等问题；致使病虫害有效监测时间段较短，监测精度较低，且难以监测早期损害。因此，尚需进一步加强基础理论研究，充分利用森林不同受害阶段的生理生化和形态指标，构造适宜的遥感监测指数，将更多数学方法融入病虫害遥感监测研究。

发明内容

本发明提供了一种基于小波分析的区域病虫害检测诊断判别方法，克服了上述现有技术之不足，其能有效解决现有的胡杨尺蠖虫害因没有有效的检测诊断判别方法造成胡杨长势逐渐衰弱的问题。

本发明的技术方案是通过以下措施来实现的：基于小波分析的区域病虫害检测诊断判别方法，包括以下步骤：

第一步，获取胡杨虫害的试验数据，对试验数据进行抽取采样分析，之后进入第二步；

第二步，采用小波变换对采样数据进行初步分析，之后进入第三步；

第三步，定义并构建NDVI时间波谱曲线，之后进入第四步；

第四步，对NDVI时间波谱曲线进行滤噪处理，之后进入第五步；

第五步，对胡杨虫害进行增强和分离分析，之后进入第六步；

第六步，建立胡杨食叶害虫的信号检测模型，之后进入第七步；

第七步，对检测出的虫害信号进行实时动态检测，提取胡杨尺蠖的危害信息。

下面是对上述发明技术方案的进一步优化或/和改进：

上述在第二步中，小波变换的分析步骤如下：

(1)将信号采用一组基函数分解成不同尺度的细节信号，所述基函数是通过将小波母函数ψ(t)进行伸缩或平移得到的，设其伸缩因子为a，平移因子为b，则平移伸缩后的函数ψ_a,b(t)为：

其中：a为尺度参数，b为平移参数。尺度参数a确定了小波函数的时域宽度，平移参数b确定了小波函数的中心位置；

(2)对于任意函数f(t)∈L²(R)连续小波变换为：

其中(W_ψf)(a,b)为小波变换系数，小波变换的实质是对信号的时-频联合分析，在分析高频信号时，时间窗变小；而在分析低频信号时，时间窗变大；

(3)对于尺度及位移均离散变化的小波序列选取a＝a₀ ^j，b＝n₃b₀a₀ ^j，其中a₀>1，b₀>0，j、n₃是整数，定义离散小波变换为：

上述在第三步中，定义并构建NDVI时间波谱曲线包括如下步骤：

(1)设周期为T，f(t)表示NDVI时间波谱曲线，则f(t)表示为：

f(t)＝f(t+n₁T) n₁＝0,1,2,…N-1 (4)

若NDVI时间波谱曲线满足条件则说明该时间波谱绝对可积，可以进行傅里叶分析和小波分析；

(2)对NDVI时间波谱曲线进行特征分析得到NDVI生长曲线波形为接近正弦或者余弦的周期性曲线；

(3)对NDVI时间波谱曲线解构采用mallat算法，运用db6小波函数对时间波谱进行多层次的小波分析，具体包括以下过程：

(a)采用mallat算法，运用db6小波函数对时间波谱进行小波分析；

(b)采用DB小波基再次对时间波谱进行多尺度分析；

(c)根据小波变换信号的高频信息和低频信息确定分解尺度；

(d)利用以下公式选择消失距：

小波的消失距定义为：

若则小波具有N阶消失距；

(4)对NDVI时间波谱曲线成分进行分析，包括以下过程：

(a)借助小波分析工具对NDVI时间波谱曲线进行成分分析，提取其趋势成分；

(b)将周期成分确认是NDVI时间波谱曲线的基波，代表胡杨正常生长的规律；

(c)确认NDVI时间波谱曲线中出现明显的锯齿现象是否由随机成分对NDVI时间波谱曲线造成干扰。

上述在第四步中，对NDVI时间波谱曲线的滤噪包括以下步骤：

(1)对NDVI时间波谱曲线的分解及低频成分的提取是通过小波变换进行分解，将高频系数过滤，提取分解后的低频系数；

(2)采用加法模型对NDVI时间波谱曲线进行嫁接与重构，加法模型公式如下：

Y＝T+C+R (6)

式中，Y为NDVI时间波谱曲线；T为长期趋势成分；C为周期成分；R为灾变成分。

上述在第五步中，对胡杨虫害信息的增强与分离包括以下步骤：

(1)通过matlab小波多尺度分解，提取高频系数，经过多尺度重构获得不同层次的小波分量，区分NDVI时间波谱曲线中的突变信息；

(2)对突变信息表现明显的小波分量作进一步分析，确定尺蠖虫害出现的时间和特征。

上述在第六步中，建立胡杨食叶害虫的信号检测模型包括以下步骤：

(1)建立判别模型，其判别函数的公式如下：

y＝b₀+b₁x₁+b₂x₂+…+b_ix_i (7)

式中：y为判别函数值；x_i为判别变量；b_i为相应的判别系数，将所要判别对象的变量代入上述公式(7)的判别函数中，求出该对象的y值并判别所属类别；

(2)通过使用贝叶斯判别方法进行判别分析，根据自变量x建立判别函数

y＝Int(-6.378+231.149x) (8)

(3)采用Press’s Q标准对判别效果进行检验，表达式为：

其中，N为样本总量，n₂为被正确分类的样本量，K为组数；

(4)采用平均误差法对检测结果进行误差分析，平均相对误差计算公式如下：

式中，i表示第i个共振点，F_i为第i个共振点发生的日期，Q_i为第i个共振点实际发生日期，N为共振点的总个数。

本发明以叶尔羌河中下游沿岸为典型区域，以胡杨尺蠖为研究对象，根据其生活史规律，结合野外采样和实地调查，基于森林病虫害遥感监测的原理与物理基础，根据胡杨食叶害虫的生活史规律，结合野外采样和实地调查，借助NDVI时间序列数据，定义并构建含有虫害信息的胡杨NDVI时间波谱曲线，通过分析波谱曲线中各生长周期内NDVI的表现特征，阐明曲线的构造特点及其组成成分，得到NDVI随时间变化的演变特性与规律。借助小波分析技术对曲线进行平滑和重整，滤除干扰信息，实现波谱曲线的滤噪，根据时间波谱曲线中存在的胡杨尺蠖虫害突变信息，利用其共振现象，突出诊断点和突变点，实现虫害信息的放大与增强，分离出虫害突变信号，并对其进行放大和增益，突出虫害信息，构建检测诊断判别模型检出虫害信号，为病虫害预测提供更精细的信息输出。

附图说明

附图1为本发明的试验数据采集区位置示意图。

附图2为本发明的流程图。

附图3为本发明胡杨NDVI时间波谱曲线。

附图4为本发明NDVI时间波谱曲线长期趋势成分示意图。

附图5为本发明NDVI时间波谱曲线周期成分示意图。

附图6为本发明NDVI时间波谱曲线含有虫害信息的随机成分示意图。

附图7为本发明NDVI时间波谱曲线未处理与滤噪后对比示意图。

附图8为本发明重构后的NDVI时间波谱曲线示意图。

附图9为本发明的小波解构示意图。

附图10为本发明2001-2014年3月份平均气温折线图。

附图11为本发明共振点幅值直方图。

具体实施方式

本发明不受下述实施例的限制，可根据本发明的技术方案与实际情况来确定具体的实施方式。

下面结合实施例及附图对本发明作进一步描述：

如附图1、2、3、4、5、6、7、8、9所示，基于小波分析的区域病虫害检测诊断判别方法，包括以下步骤：

第一步，获取胡杨尺蠖虫害的试验数据，对试验数据进行抽取采样分析，之后进入第二步；

第二步，采用小波变换对采样数据进行初步分析，之后进入第三步；

第三步，定义并构建NDVI时间波谱曲线，之后进入第四步；

第四步，对NDVI时间波谱曲线进行滤噪处理，之后进入第五步；

第五步，对胡杨尺蠖虫害进行增强和分离分析，之后进入第六步；

第六步，建立胡杨食叶害虫的信号检测模型，之后进入第七步；

第七步，对检测出的虫害信号进行实时动态检测，提取胡杨尺蠖的危害信息。

这里第一步中的试验数据为野外调研数据，研究区位于夏河和夏马勒林场，地处天山南麓，塔里木河主要支流叶尔羌河中下游冲积平原，地理坐标东经77° 22′ 30″ -79°56′ 15″，北纬38° 47′ 30″ -40° 17′ 30″，东西长218km，南北宽134km，总面积2.1741×10⁴km²。根据研究区的地形地貌和胡杨分布状况，按照典型抽样原则，选取7个样点，用GPS记录各样点的经纬度坐标及高程，并于2009-2015年的每年2月24日—5月10日到各样点进行详细调查，主要记录每日蛹(总蛹数、活蛹、死蛹)、羽化成虫、产卵总数、孵化、幼虫(1龄、2龄)、地表温度等指标。在塔里木盆地严酷的气候背景下，胡杨尺蠖一年只发生一代，且大多数幼虫只能生长到1龄虫，少量可到2龄虫，一般3月下旬幼虫开始危害，到4月中下旬结束。此后至洪水期，受害胡杨二次萌芽，开始正常生长。因此，根据这一生活史规律，可以构建基于NDVI的胡杨生长时间波谱曲线，并从中提取尺蠖虫害危害信息。

这里使用归一化植被指数即NDVI获取试验数据记录，NDVI依据多光谱数据，利用植物光谱中的近红外与可见光红波段两个典型的波段值，经线性和非线性组合构成的对植被有一定指示意义的各种数值。NDVI作为一个重要的遥感参数，是检测植被覆盖度、植被长势的最佳因子。植被指数与生长量、叶面积指数、叶绿素含量等有较强的相关性，植物的长势、覆盖度及季相动态变化直接对应植被指数的数量变化，因此，植被指数不仅可以应用于植被分类，还可以应用于检测植被覆盖度、生物量估算、监测植物的长势以及监测病虫害等。

如附图10所示，这里选用2001～2014年的MOD13Q1第1波段即“250m16-day NDVI”用来构建NDVI时间波谱曲线，检测研究区胡杨林尺蠖虫害信息，由于该数据为16天最大值合成数据，因此最大程度的排除了云、沙尘暴等干扰信息的影响。覆盖研究区的MOD13Q1影像有4幅，轨道号为h23v04，h23v05，h24v04，h24v05，每年23期，共计322期1288幅。对MOD13Q1原始数据进行进行镶嵌、投影变换、重采样后，利用研究区胡杨林资源矢量对影像进行掩膜处理，获取研究区内胡杨林的NDVI数据，同时利用提取工具得到各采样点每期的NDVI值，以构建典型区域的NDVI时间波谱曲线。本发明基于森林病虫害遥感监测的原理与物理基础，根据胡杨食叶害虫的生活史规律，结合野外采样和实地调查，借助NDVI时间序列数据，定义并构建含有虫害信息的胡杨NDVI时间波谱曲线，通过分析波谱曲线中各生长周期内NDVI的表现特征，阐明曲线的构造特点及其组成成分，得到NDVI随时间变化的演变特性与规律。借助小波分析技术对曲线进行平滑和重整，滤除干扰信息，实现波谱曲线的滤噪，根据时间波谱曲线中存在的胡杨尺蠖虫害突变信息，利用其共振现象，突出诊断点和突变点，实现虫害信息的放大与增强，分离出虫害突变信号，并对其进行放大和增益，突出虫害信息，构建检测诊断判别模型，检出虫害信号。

可根据实际需要，对上述基于小波分析的区域病虫害检测诊断判别方法作进一步优化或/和改进：

如附图1、2、3、4、5、6、7、8、9所示，在第二步中，小波变换的分析步骤如下：

(1)将信号采用一组基函数分解成不同尺度的细节信号，所述基函数是通过将小波母函数ψ(t)进行伸缩或平移得到的，设其伸缩因子为a，平移因子为b，则平移伸缩后的函数ψ_a,b(t)为：

其中：a为尺度参数，b为平移参数；尺度参数a确定了小波函数的时域宽度，平移参数b确定了小波函数的中心位置；

(2)对于任意函数f(t)∈L²(R)连续小波变换为：

其中(W_ψf)(a,b)为小波变换系数，小波变换的实质是对信号的时-频联合分析，在分析高频信号时，时间窗变小；而在分析低频信号时，时间窗变大；

这里连续小波系数具有较大的冗余量，这种冗余量的存在加大了计算量，因此为了能在计算机中实现并减少计算量，对连续小波必须进行离散化。

(3)对于尺度及位移均离散变化的小波序列选取b＝n₃b₀a₀ ^j，其中a₀>1，b₀>0，j、n₃是整数，定义离散小波变换为：

如附图1、2、3、4、5、6、7、8、9所示，在第三步中，定义并构建NDVI时间波谱曲线包括如下步骤：

(1)设周期为T，f(t)表示NDVI时间波谱曲线，则f(t)表示为：

f(t)＝f(t+n₁T) n₁＝1,2,3,… (4)

若NDVI时间波谱曲线满足条件则说明该时间波谱绝对可积，可以进行傅里叶分析和小波分析；

这里的傅里叶分析和小波分析的实质是把f(t)这个波形分解成许多不同频率的正弦波或者余弦波的叠加。每个正弦和或余弦和乘以不同的系数。无论函数多么复杂，只要是周期的，并且满足了某些数学条件，都可以用这样的和来表示。因此NDVI时间波谱满足以下性质：线性性质、位移性质、微分性质、积分性质、乘积定理、能量积分。

(2)如附图3所示，对NDVI时间波谱曲线进行特征分析得到NDVI生长曲线波形为接近正弦或者余弦的周期性曲线；

这里根据分析，胡杨随着气候的季节性变化，从春季开始萌芽展叶，对应的NDVI值开始逐渐上升，到胡杨生长最旺盛的时期，对应的NDVI值达到最大值，然后随着胡杨叶片的变黄、枯萎，NDVI值开始逐渐下降，到了冬季胡杨进入休眠期，NDVI值出现最小值，整体上波形接近正弦或者余弦曲线。

在附图3中的胡杨NDVI时间波谱曲线存在明显的锯齿，许多点都存在突升或突降的现象，究其原因，胡杨在生长过程中，由于受到水汽、云、雪、干旱、火灾、病虫害、收割、采伐、地表的反射以及大气的不确定性影响，导致数据中存在许多噪声，这些噪声影响着NDVI能量值的大小，反映到时间波谱曲线上，就形成了突升或突降的锯齿，与胡杨正常的生长曲线不符。

(3)对NDVI时间波谱曲线解构采用mallat算法，运用db6小波函数对时间波谱进行多层次的小波分析，具体包括以下过程：

(a)采用mallat算法，运用db6小波函数对时间波谱进行小波分析；

在1989年，mallat在小波变换多分辨率分析理论与图像处理的应用研究中受到塔式算法的启发，提出了信号的塔式多分辨率分析与重构的快速算法称为马拉特(mallat)算法。其解构过程中，每经过一次解构，近似序列被解构为低一级的近似序列和细节序列，但两者的长度均为输入序列的一半，即数据总量保持不变。因此本发明选择马拉特算法进行解构。在小波分解过程中，本发明利用matlab7.0软件，采用mallat算法，运用db6小波函数对时间波谱进行8层小波分析。基于小波分析原理，当序列样本容量为N时，解构的层数最多不能超过log₂N层，本发明中时间波谱样本容量为N＝322，故最多可以解构至第8([log₂322])层，通过对比试验的解构标准误差，确定经过8次解构，误差标准差最小，且能满足对趋势变化以及突变成分进行提取的需要。

(b)如表2所示，采用DB小波基进行多尺度分析；

在比较小波基的特性时，正交性决定了操作的简单和便于理解的程度；紧支撑性决定了小波基的时频局部化能力和衰减的快慢；对称性能避免因相位的偏差而引起的失真；消失距的大小决定了频域局部化能力的大小等。本发明根据小波基的特性，进行多次试验，通过试验结果的满意程度最终DB小波是小波分析学者Inrid Daubechies构造的小波函数，简写为dbN，N是小波的阶数。小波Ψ(t)和尺度函数φ(t)中的支撑区为2N-1，Ψ(t)的消失矩为N。除N＝1外，dbN不具有对称性(即非线性相位)。dbN没有明确的表达式(除N＝1外)，但转换函数h的平方模是明确的。

(c)根据小波变换信号的高频信息和低频信息确定分解尺度；

小波变换的实质就是通过多尺度分析将信号分为高频信息和低频信息。理论上，一个数据长度为N的信号可进行log₂N层的分解，而每次分解后的数据长度会降低一半。在小波分析过程中，分解的层数越多，信号的高频信息和低频信息就会分解的越彻底，但是计算量随着分解层数的增多变得越大，就要占用一定的存储空间。另一方面，分解层数的多少严重影响去噪效果。本研究经过多次测试，发现当分解尺度为8时，既能保存信号中的有用成分，也能有效的去除无用成分，因此最终确定分解尺度为8。

(d)选择消失距利用以下公式：

小波的消失距定义为：

若则小波具有N阶消失距；

当小波具有N阶消失距时，所有多项式的小波系数均为零，也就是小波分解的高频信息部分均为0，由此可见，在进行突变信号的分析时，为了能有效地快速地检测出突变点，小波基应具有尽可能高的消失距；经过试验，本发明最终选取的消失距为6。本发明借助MODIS13Q1影像，定义并构建胡杨NDVI时间波谱曲线，根据NDVI时间波谱曲线分析各生长周期和每个生长周期内NDVI的表现特征，阐明曲线的构造特点。分析并阐明其组成成分，得到NDVI随时间变化的演变特性与规律。这里的NDVI能量信号随着时间的变化形成NDVI时间波，把NDVI时间波按照波长或频率，递增或递减排列，则构成了NDVI时间波谱。对于MODISNDVI影像数据中任意的一个像元来说，其NDVI值的大小就表示NDVI能量的大小。同一像元的NDVI能量在不同的时间刻度是不同的，用x轴表示时间t，y轴表示NDVI能量，则某一像元在某一时刻的NDVI能量用笛卡尔坐标表示。

本发明利用胡杨在NDVI能量特征上的响应，借助2001-2014年共322期的胡杨NDVI时间序列数据，时间采样间隔为16天，输入NDVI能量，输入时间T，构建研究区长时间序列含有虫害信息的胡杨NDVI时间波谱曲线(附图3)，确定其波长、频率、振幅、周期等相关指标，并进行形态分析。

(4)对NDVI时间波谱曲线成分进行分析，包括以下过程：

(a)如附图4所示，借助小波分析工具对NDVI时间波谱曲线进行成分分析，提取其趋势成分；

这里的长期趋势成分指在一定的时间尺度范围内，在自然地理演进过程中起主要和决定性作用的因素，这类因素使得过程的变化趋势长期沿着一定的方向发展，呈现出某种长期的变化趋势。NDVI在2001至2012年期间呈平稳上升趋势，2012年以后则略有下降。说明胡杨生长环境有所改善，引起NDVI的上升。经研究发现，气候暖湿化趋势比较明显有助于植物生长，本文中NDVI上升趋势与气候暖湿化趋势相吻合。因此NDVI时间波谱曲线中的趋势成分在一定程度上反映了该区域的气候变化趋势。

(b)如附图5所示，将周期成分确认是NDVI时间波谱曲线的基波，代表胡杨正常生长的规律；

NDVI时间波谱曲线的最高峰和最低峰出现的时间符合胡杨的生长规律，周期性成分指从低到高或者从高到低的周而复始的一种有规律的变动。周期性波动是升降相间、涨落交替的变动，与趋势变动不同，它不是沿着单一的方向持续运动。

(c)如附图6所示，确认NDVI时间波谱曲线中出现明显的锯齿现象是否由随机成分对NDVI时间波谱曲线造成干扰。

随机成分指在正常的地理过程发生的一些不规则、偶然的、变动的因素，由于它们的影响使得过程的发展变化呈现出无规律的、不规则的状态。表现在NDVI时间波谱曲线中是由旱灾、火灾、沙尘暴、病虫害、收获、采伐等引起NDVI的突降，在曲线中出现明显的锯齿现象。

其中旱灾一般出现在夏季(每年的7月-8月)，病虫害一般出现在每年的春季(3月-4月)，这两种灾害属于可以确定的灾害因素。而沙尘暴、火灾、采伐等属于不确定的灾变因素，虽然在每年频繁发生，但是没有具体的规律可循。另外，影响NDVI时间波谱曲线出现锯齿现象的还有云、雪、水汽、地表的反射以及大气的不确定性等因素，这些因素都属于随机因素，都会对NDVI时间波谱曲线造成干扰。

如附图1、2、3、4、5、6、7、8、9所示，在第四步中，对NDVI时间波谱曲线的滤噪包括以下步骤：

(1)对NDVI时间波谱曲线的分解及低频成分的提取是通过小波分解，将高频系数过滤，提取分解后的低频系数；

这里对胡杨NDVI时间波谱曲线进行小波分解，提取滤除了干旱、火灾、病虫害、收获采伐、沙尘暴等随机因素影响的低频成分曲线；小波分析法对NDVI时间波谱曲线去噪实际上是对信号进行低通滤波的过程。

本发明选择db6小波基函数并确定消失距与分解的层次，对NDVI时间波谱曲线进行分解；其次将分解后的高频系数自动过滤，只提取各个层的低频系数；最后对提取的各个层的低频系数采用mallat算法进行多尺度重构，获取去噪后时间波谱曲线。

(2)采用加法模型对NDVI时间波谱曲线进行嫁接与重构，加法模型公式如下：

Y＝T+C+R (6)

式中，Y为NDVI时间波谱曲线；T为长期趋势成分；C为周期成分；R为灾变成分。

本发明中将此低频成分曲线与3月下旬至5月上旬幼虫危害时段的原始胡杨时间波谱曲线成分进行嫁接移植，重构得到既含有春季尺蠖危害信息又滤除了其他时段随机成分的胡杨NDVI时间波谱曲线。

如附图1、2、3、4、5、6、7、8、9所示，在第五步中，对胡杨虫害信息的增强与分离包括以下步骤：

(1)通过matlab小波多尺度分解，提取高频系数，经过多尺度重构获得不同层次的小波分量，区分NDVI时间波谱曲线中的突变信息；

如附图10所示，胡杨生长季一般开始于3月上旬(即每年的t5时期)，结束于11月下旬(即每年的t21时期)，而尺蠖幼蛹于3月中旬孵化，3月下旬到4月是其危害盛期。分析2001年NDVI时间波谱曲线与通过小波解构后重整的细节分量D1可知，胡杨生长期开始后NDVI值迅速上升，到t6时期(3月下旬)却突然下降，所对应的细节分量D1波动幅值为0.0240，到t7(4月上旬)时期又继续上升，经实地调查发现，造成胡杨NDVI值突然下降的原因是尺蠖啃食胡杨叶片，使胡杨叶片残缺不全，进而危害胡杨林，导致该时期NDVI时间波谱下降。

(2)对突变信息表现明显的小波分量作进一步分析，确定尺蠖虫害出现的时间和特征。

这里，理论上NDVI时间波谱曲线应与其植被生长曲线相吻合，如果出现NDVI值突然下降而导致与其植被生长曲线相反，则下降点可以称为突变点。通常来说突变部分反应了检测对象的某些重要特征，含有较丰富的高频信号，而小波具有自适应的时-频局部化功能，通过多尺度分解后，使NDVI时间波谱曲线的突变部分在某些高频分量上的表现幅度增大，其与正常信息在高频部分的表现正好形成鲜明的对比，从而能有效的区分波谱曲线中的突变信息。

通过matlab小波多尺度分析，可得到第8层近似系数(低频系数)和第8层近似分量(低频分量)，以及8层细节系数(高频系数)与这8层细节系数重整的细节分量(高频分量)。

如附图9所示：D8、D7、D6、D5、D4、D3、D2、D1按8层不同的细节系数重构后的8层细节分量，A8是近似系数重构后的第8层近似分量，若x是被分解的NDVI时间波谱曲线，那么x可表示为x＝A8+D8+D7+D6+D5+D4+D3+D2+D1，db6小波函数的解构层数越高，接收的时间波谱信号的近似信息(低频信息)和细节信息(高频信息)就越清晰，且近似信息越接近原始信号。经过对比发现第1层的细节信息(高频信息)是最清楚的，也是最能反应突变信息的一层细节分量，所以将第一层细节信息(高频信息)单独提取出来，并进行详细的分析，进而判断时间波谱曲线中突变信息的位置和时刻，即胡杨尺蠖虫害信息发生的具体时间和特征。

由附图9所示，在第一层细节分量信息(细节分量D1)中，NDVI时间波谱幅值分布不均匀，但其大部分都分布在一定的区域之内，每个频率上的值都是对称分布的。时间波谱中的细节信息里存在突变信息，这些突变信息出现在不同的时间位置上，在波动幅值较大的地方，出现“共振点”，“共振点”的峰值有大有小，本发明可以通过“共振点”判定胡杨尺蠖出现的时间。其中，第一个“共振点”出现在第6期，第二个“共振点”出现在第121期，第三个“共振点”出现在第146期，第四个“共振点”出现在167期，第五个“共振点”出现在第190期，第六个“共振点”出现在第213期，第七个“共振点”出现在第238期，第八个“共振点”出现在第260期，第九“共振点”出现在第305期。在这几期NDVI时间波谱中，造成“共振点”出现的原因是NDVI值突然下降引起的，而“共振点”幅值上下波动的大小则是NDVI值下降的大小引起的。经过对比和调查，发现突变点出现的时间与胡杨尺蠖孵化啃食胡杨树叶发生的时间基本吻合，所以对以上几个“共振点”进行具体分析。

胡杨NDVI时间波谱中“共振点”出现的时间位置，遥感影像上表现出NDVI值突然下降，通过matlab小波解构，在细节分量D1上能明显的识别出“共振点”信息。共振点的幅值有正有负，主要原因是本期的NDVI能量值与其前一期、后一期的NDVI能量值的差异大小决定的，如果本期的NDVI值与前一期的NDVI能量值差异较大，那就共振点的幅值为正，如果与后一期的NDVI能量值差异较大，则共振点的幅值为负。根据共振点幅值的大小对尺蠖虫害的危害程度进行分级，共分为4级，即0级：无虫害；1级：轻度发生；2级：中度发生；3级：重度发生。

如附图11所示，对细节分量D1进行直方图统计可发现，其波动幅值在0.00423左右的期数为232期，波动幅值在0.01014-0.0150范围的期数为24期，波动幅值在-0.0017--0.0150范围的期数为48期，波动幅值大于0.0150的期数为9期，波动幅值小于-0.0150的期数为9期。根据此统计范围，将出现“共振点”时其幅值的波动大小与直方图统计范围进行对比，“共振点”幅值处在大于0.0150、小于-0.0150的范围之中，而出现“共振点”的期数为9期，与直方图统计信息相符，根据此分析结果可得，波动幅值在大于-0.0150、小于0.0150范围内的没有病虫害，波动幅值在大于0.0150、小于-0.0150范围内的有病虫害信息。

根据共振点幅值的大小对尺蠖虫害的危害程度进行分级，共分为4级，即0级：无虫害；1级：轻度发生；2级：中度发生；3级：重度发生。共振点的幅值也相应的分为4级，其标准见表4所示。

如附图1、2、3、4、5、6、7、8、9所示，在第六步中，建立胡杨食叶害虫的信号检测模型包括以下步骤：

(1)建立判别模型，其判别函数的公式如下：

y＝b₀+b₁x₁+b₂x₂+…+b_ix_i (7)

式中：y为判别函数值；x_i为判别变量；b_i为相应的判别系数，将所要判别对象的变量代入上述公式(7)的判别函数中，求出该对象的y值，然后再判别它应该属于哪个类别；

这里的判别分析是根据一批已知明确类别的样品观察数据的多项指标，制定出一个分类标准，用以指导未知类别的个体和归类。判别分析的基本模型就是判别函数，表示分组变量与判别变量的线性函数关系，常用的判别分析法主要包括距离判别法、贝叶斯判别法、费希尔判别法。本发明使用贝叶斯判别方法进行判别分析。根据尺蠖虫害分级标准，利用spss软件进行模型计算，得到如下表：

表5组统计量是各组和总体的每个变量的描述统计分析。

表6反映了判别函数的特征值、解释方差的比例和典型的相关系数。

表7是对判别函数的显著性检验，结果显示两个判别函数均具有统计学上的意义，即Sig小于0.05。

表8是非标准化判别系数，即费歇尔判别系数，是用来计算判别函数值的，将个体在各判别变量上的取值代入判别函数式，就可以计算出个体的判别函数值，根据此函数值决定个体所属的类别。

(2)通过使用贝叶斯判别方法进行判别分析，根据自变量x建立判别函数

y＝Int(-6.378+231.149x) (8)

此模型有效地将不同的监测点进行正确分类；这里的92.9％的变量被正确分类。

(3)采用Press’s Q标准对判别效果进行检验，表达式为：

其中，N为样本总量，n₂为被正确分类的样本量，K为组数；

这里的Press’s Q是一种将分类准确性与随机过程相比较的标准，这个指标在0.01显著水平上的临界值为6.63，在0.05显著水平上的临界值为3.84(服从自由度为1的卡方分布)。

根据模型分类结果可知，样本总量为14，正确分类的样本量为13，错误分类的样本量为1，组数为4，将数据代入式(9)得：

判别值为34.4，大于在0.01显著水平上的临界值6.63，判别显著，符合判别标准，确认判别模型可信；

根据判别模型分类结果可知，在2001-2004年期间，未发生尺蠖虫害的年份为2002年、2003年、2004年、2005年、2013年；轻度发生尺蠖虫害的年份为2001年、2014年；中度发生尺蠖虫害的年份为2008年、2009年、2010年、2011年、2012年；重度发生尺蠖虫害的年份为2006年、2007年。通过现地调查验证，判别结果与现地尺蠖发生的情况基本一致，判别结果可信。

(4)采用平均误差法对检测结果进行误差分析，平均相对误差计算公式如下：

式中，i表示第i个共振点，F_i为第i个共振点发生的日期，Q_i为第i个共振点实际发生日期，N为共振点的总个数。

通过小波分析对重构后的含有虫害信息的NDVI时间波谱曲线进行增强和放大，获取共振点，进而确定胡杨尺蠖发生的具体时间。一般来说，在获取共振点确定尺蠖发生的时间后，还必须对小波分析获取的发生时间的结果做检验评估，只有通过检验才能鉴定我们所使用的方法的优劣，并加以改进。

一般来说，获取的共振点的具体发生日期的准确度与误差成反比。误差越小，表示检测结果越准确。当误差接近0时，检测结果与实况差距最小，检测近于完全正确。一般认为在20％以下时，检测结果尚可；当在10％以下，检测结果较准确。将表10中的数据代入式(10)中，得等于8.76％，说明运用小波分析法对尺蠖虫害信息发生的时间和特征的提取结果较为准确。

以上技术特征构成了本发明的实施例，其具有较强的适应性和实施效果，可根据实际需要增减非必要的技术特征，来满足不同情况的需求。

表1 2009-2015年尺蠖生活史表

表2常见小波基函数的主要特性

表3共振点幅值统计表

表4 尺蠖虫害分级标准

级别	共振点幅值
		0	<0.015
1	0.015～0.030
		2	0.031～0.0450
3	>0.045

表5 组统计量

表6 特征值

表7 Wilks的Lambda

函数检验	Wilks的Lambda	卡方	df	Sig.
					1	.052	30.977	3	.000

表8 典型判别式函数系数

表9分类结果

表10尺蠖发生时间统计表

年份	发生范围	检测发生日期	实际发生日期
				2001	3月下旬	22	21
2006	3月下旬	22	19
				2007	4月下旬	23	21
2008	3月下旬	21	22
				2009	3月下旬	22	21
2010	3月下旬	22	23
				2011	4月下旬	23	21
2012	4月上旬	6	4
				2014	3月下旬	22	19
平均值		20.3	19

Claims (5)

1.一种基于小波分析的区域病虫害检测诊断判别方法，其特征在于包括以下步骤：

第一步，获取胡杨尺蠖虫害的试验数据，对试验数据进行抽取采样分析，之后进入第二步；

第二步，采用小波变换对采样数据进行初步分析，之后进入第三步；

第三步，定义并构建NDVI时间波谱曲线，包括如下步骤：

(1)设周期为T，f(t)表示NDVI时间波谱曲线，则f(t)表示为：

f(t)＝f(t+n₁T) n₁＝0,1,2,…N-1 (4)

若NDVI时间波谱曲线满足条件则说明该时间波谱绝对可积，能够进行傅里叶分析和小波分析；

(2)对NDVI时间波谱曲线进行特征分析得到NDVI生长曲线波形为接近正弦或者余弦的周期性曲线；

(3)对NDVI时间波谱曲线解构采用Mallat算法，运用db6小波函数对时间波谱进行多层次的小波分析，具体包括以下过程：

(a)采用Mallat算法，运用db6小波函数对时间波谱进行小波分析；

(b)采用DB小波基再次对时间波普进行多尺度分析；

(c)根据小波变换信号的高频信息和低频信息确定分解尺度；

(d)利用以下公式选择消失距：小波的消失距定义为：

若则小波具有N阶消失距；

(4)对NDVI时间波谱曲线成分进行分析，包括以下过程：

(a)借助小波分析工具对NDVI时间波谱曲线进行成分分析，提取其趋势成分；

(b)将周期成分确认是NDVI时间波谱曲线的基波，代表胡杨正常生长的规律；

(c)确认NDVI时间波谱曲线中出现的明显锯齿现象是否由随机成分对NDVI时间波谱曲线造成干扰；

之后进入第四步；

第四步，对NDVI时间波谱曲线进行滤噪处理，之后进入第五步；

第五步，对胡杨尺蠖虫害进行增强和分离分析，之后进入第六步；

第六步，建立胡杨食叶害虫的信号检测模型，包括以下步骤：

(1)建立判别模型，其判别函数的公式如下：

y＝b₀+b₁x₁+b₂x₂+…+b_ix_i (7)

式中：y为判别函数值；x_i为判别变量；b_i为相应的判别系数，将所要判别对象的变量代入上述公式(7)的判别函数中，即可求出该对象的y值并判别所属类别；

(2)通过使用贝叶斯判别方法进行判别分析，根据自变量x建立判别函数

y＝Int(-6.378+231.149x) (8)

(3)采用Press’s Q标准对判别效果进行检验，表达式为：

其中，N₁为样本总量，n₂为被正确分类的样本量，K为组数；

(4)采用平均相对误差法对检测结果进行误差分析，平均相对误差计算公式如下：

式中，i表示第i个共振点，F_i为第i个共振点发生的日期，O_i为第i个共振点实际发生日期，N₂为共振点的总个数；之后进入第七步；

第七步，对检测出的虫害信号进行实时动态检测，提取胡杨尺蠖的危害信息。

2.根据权利要求1所述的基于小波分析的区域病虫害检测诊断判别方法，其特征在于第二步中，小波变换的分析步骤如下：

(1)将信号采用一组基函数分解成不同尺度的细节信号，所述基函数是通过将小波母函数ψ(t)进行伸缩或平移得到的，设其伸缩因子为a，平移因子为b，则平移伸缩后的函数ψ_a,b(t)为：

其中：a为尺度参数，b为平移参数；尺度参数a确定了小波函数的时域宽度，平移参数b确定了小波函数的中心位置；

(2)对于任意函数f(t)∈L²(R)连续小波变换为：

其中，(W_ψf)(a,b)为小波变换系数，小波变换的实质是对信号的时-频联合分析，在分析高频信号时，时间窗变小；而在分析低频信号时，时间窗变大；

(3)对于尺度及位移均离散变化的小波序列选取a＝a₀ ^j，b＝n₃b₀a₀ ^j，其中a₀>1，b₀>0，j、n₃是整数，定义离散小波变换为：

3.根据权利要求1或2所述的基于小波分析的区域病虫害检测诊断判别方法，其特征在于第四步中，对NDVI时间波谱曲线的滤噪包括以下步骤：

(1)对NDVI时间波谱曲线的分解及低频成分的提取是通过小波变换进行分解，将高频系数过滤，提取分解后的低频系数；

(2)采用加法模型对NDVI时间波谱曲线进行嫁接与重构，加法模型公式如下：

Y＝T+C+R (6)

式中，Y为NDVI时间波谱曲线；T为长期趋势成分；C为周期成分；R为灾变成分。

4.根据权利要求1或2所述的基于小波分析的区域病虫害检测诊断判别方法，其特征在于第五步中，对胡杨尺蠖虫害信息的增强与分离包括以下步骤：

(1)通过matlab小波多尺度分解，提取高频系数，经过多尺度重构获得不同层次的小波分量，区分NDVI时间波谱曲线中的突变信息；

(2)对突变信息表现明显的小波分量作进一步分析，确定尺蠖虫害出现的时间和特征。

5.根据权利要求3所述的基于小波分析的区域病虫害检测诊断判别方法，其特征在于第五步中，对胡杨尺蠖虫害信息的增强与分离包括以下步骤：

(1)通过matlab小波多尺度分解，提取高频系数，经过多尺度重构获得不同层次的小波分量，区分NDVI时间波谱曲线中的突变信息；

(2)对突变信息表现明显的小波分量作进一步分析，确定尺蠖虫害出现的时间和特征。