CN105023011A - 一种基于hmm的农作物物候动态估计方法 - Google Patents

Abstract

本发明涉及一种基于HMM的农作物物候动态估计方法，包括：建立HMM混合模型，获得包括：初始概率、局部状态转移矩阵、观测概率矩阵；向HMM混合模型输入农作物的多源特征，多源特征包括：NDVI均值、分维值、有效积温；通过HMM混合模型计算出农作物各发育阶段的百分比。本发明采用的分维值为遥感影像的纹理特征，可用于指示农作物生育过程的变化，可直接采用提取的多远特征作为数据输入，无须额外的过滤或平滑处理，农作物发育阶段百分比的精度在±12.91%至±16.14%之间，而遥感光谱逐像元方法最优的为±24.6%，所获得的结果较好，可实现玉米作物物候信息的动态估计。

Description

一种基于HMM的农作物物候动态估计方法

技术领域

本发明涉及一种基于HMM的农作物物候动态估计方法。

背景技术

目前，基于分形的农作物物候遥感特征提取受传统地面作物物候信息单点调查手段费时、费力，且无法大范围操作的限制，遥感手段成为农作物物候信息提取的热点。建立遥感影像纹理特征与农作物物候的联系，是遥感物候检测首当其冲的。分维作为影像纹理粗糙度的一种表述，可反映农作物发育过程中NDVI影像纹理的变化。建立分维与农作物物候的关联，可为物候检测提供理论基础。另外，考虑到遥感影像上农作物耕地呈块状零星分布的特点，需利用分形乘积的原理，设计并实现一种针对遥感影像不规则ROI的分维估计算法。

目前，大多数作物物候检测方法，仅依赖单源特征(比如：NDVI指数或有效积温等)，且只能检测出少量特定的物候期(比如：变绿期、成熟期、衰落期和休眠期)。集成遥感光谱值、影像纹理和地面气象站变量等多源多特征，有望提高农作物物候检测的种类、准确度和实时化程度。

发明内容

本发明的目的是提供一种基于HMM的农作物物候动态估计方法。

为达到上述目的，本发明采用的技术方案是：

一种基于HMM的农作物物候动态估计方法，包括：

1)：建立HMM混合模型，HMM混合模型中包括获得：

a、初始概率：农作物在各发育阶段发生的概率；

b、局部状态转移矩阵：农作物在整个生育周期过程中的各发育阶段转移概率；

c、观测概率矩阵：在给定发育阶段状态的情况下，各观测变量发生的概率；

2)：向HMM混合模型输入农作物的多源特征，多源特征包括：

a、NDVI均值；

b、分维值；

c、有效积温；

3)、通过HMM混合模型计算出农作物各发育阶段的百分比。

优选地，初始概率通过统计该地区同一时期的历史数据获得。

优选地，获取局部状态转移矩阵包括：

1)、选取同一时间点的多年归一化地面物候调查数据数据，计算其均值；

2)、计算出每个时间点的状态转移概率矩阵。

优选地，获取NDVI均值包括：

1)、影像合成：将每日的MODIS NDVI时间序列影像合成为每周的NDVI时间序列影像；

2)、影像掩膜：剔除MODIS NDVI时间序列影像上非农作物类型的像元。

进一步优选地，对剔除MODIS NDVI时间序列影像上非农作物类型的像元进行估算获得。

优选地，有效积温获取包括：获取离散点位的地表日最低和最高温度；将离散点位的地表日最低和最高温度转化为行政区划级别面状温度数据；对行政区划级别的日最低和最高温度进行日温校正，并计算出有效积温。

其中：

初始概率(或状态先验概率)：表示起始时刻，不考虑观测值的情况下，各发育阶段发生的概率。

局部状态转移矩阵：作物在整个生育周期过程中的各发育阶段转移概率随时间变化而改变。

观测概率：在给定发育阶段状态的情况下，各观测变量发生的概率。

NDVI指数：根据植被在遥感影像中光谱响应的差异及动态变化原理，作为衡量植被绿度及生长状况的指数。

分维值：用于衡量NDVI影像中农作物类型像元的粗糙度，并作为影像纹理异质性的指标来反映农作物的生长状况。

有效积温：在作物生长发育所需要的其他条件均得到满足时，在一定温度范围内，气温和发育速度成正相关，并且要积累到一定的温度总和，才能完成其发育期，这个温度的累积数称为积温；每种作物都有其生长的下限温度，当温度高于下限温度时，它才能生长发育，这个对作物生长发育起有效作用的高出的温度值，称作有效温度；作物在整个生育期内有效温度的累积，即有效积温。

由于上述技术方案运用，本发明与现有技术相比具有下列优点和效果：

本发明采用的分维值为遥感影像的纹理特征，可用于指示农作物生育过程的变化，可直接采用提取的多远特征作为数据输入，无须额外的过滤或平滑处理，农作物发育阶段百分比的精度在±12.91％至±16.14％之间，而遥感光谱逐像元方法最优的为±24.6％，所获得的结果较好，可实现玉米作物物候信息的动态估计。

附图说明

附图1为本实施例中HMM混合模型；

附图2为本实施例中玉米在整个生育周期过程中的各发育阶段转移示意图。

具体实施方式

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

本实施例以玉米为例，其物候动态估计方法包括：

1)、建立HMM混合模型，

假设HMM混合模型参数为λ(N,M,∏,A,B)：

T：观测值长度，对应每个时间周期的周数，

N：状态数目，对应物候期的数目，

M：每个状态可能的观测值数，对应特征个数，

S：隐形状态S＝{S₁,…,S_N}，

q：状态序列q＝{q₁,…,q_T}，

O：观测值序列O＝{O₁,…,O_T}，

∏：初始时刻状态的概率分布，比如S_j状态的概率可记为π_j，

A：状态转移矩阵(transition probability matrix)，可以表示成a_i,j＝P(q_t＝S_j|q_t-1＝S_i)，其中(1≤i,j≤N)，对应某时刻从一物候期转换到另一物候期的概率。

B：观测值的概率分布(emission probability matrix)，可以表示成b_j(t)＝P(O_t|q_t＝S_j)，(2≤t≤T)，对应某一物候期、特征的分布概率。

如图1为本实施例中HMM混合模型，其中本实施例中隐形状态有8个，它们依次为预备期(pre-season)、播种期(planted)、出苗期(emerged)、抽丝期(silking)、乳熟期(dough)、腊熟期(dent)、完熟期(mature)和收割期(harvested)。其中，预备期被定义为玉米种子播种前的一段时期，它仅仅提供模型设计时的便利性，不产生最终的结果，根据模型变量的定义，全概率可以表示为：

P(q₁＝S_r,…,q_t＝S_j,O₁,…,O_t|λ)＝P(O₁,…,O_t|q₁＝S_r,…,q_t＝S_j|λ)·P(q₁＝S_r,…,q_t＝S_j|λ)  (1-1)；

假设:1)O₁,…,O_t观测值序列是一个典型的一阶(first-order)马尔可夫链，也就是说q_t仅由q_t-1决定；2)t时刻的观测值O_t只由t时刻的状态q_t决定(独立输出假设)，为此对于全概率：

$P(q_1=S_r,\mathrm{...},q_t=S_j,O_1,\mathrm{...},O_t|\mathrm{\λ})={\mathrm{\π}}_{q_1}\·\underset{k=2}{\overset{t}{\Π}}{a}_{q_{t-1,}{q}_t}\left(k\right)\·\underset{k=1}{\overset{t}{\Σ}}{b}_{q_t}\left(O_k\right)---(1-2);$

由于HMM可以看成是马尔可夫模型和混合模型的组合，故观测值可以看做是由N个状态共同作用而成，HMM混合模型中，两个嵌入的随机过程对应两条链：显性链(观测值)和隐性链(生育期)，该随机过程反映了作物生育期监测和数据观测本身存在的不确定性，每个时间节点上的观测值都可以看做是多个隐含物候期共同作用的结果，假设观测量形成N个簇，那么其可以构模成N个组分的混合，记每个组分的密度分布函数为b_i(O_t)，那么观测值概率可以表示为：

$P\left(O_t\right)=\underset{i=1}{\overset{N}{\Σ}}{\mathrm{\π}}_i\left(t\right)\·b_i\left(O_t\right)---(1-3);$

其中，π_i(t)可以看作是第i个组分的权重，

HMM混合模型至少包含三种概率：初始概率、局部状态转移矩阵、观测概率矩阵，其中：

初始概率：玉米的播种期主要集中在第14、15周，因此，第13周的玉米作物一般处于预备期，初始概率可以通过统计该地区在同一时期的历史记录，比如：起始时间若为第13周，则计算同期可用地面物候调查数据记录的均值，将其作为初始概率。

局部状态转移矩阵：随玉米作物的生长而变化，例如：在玉米作物生长的初期，其最有可能出现在生长初期对应的物候期，在转移矩阵中，从当前发育阶段转移到其自身的概率比较高，转移到下一发育阶段的概率比较低，但是随着时间的推移，转移到下一发育阶段的能力逐渐增强，而转移到自身的能力逐渐减弱，直至消亡。一般来说，作物的物候期转移变化依赖于生物物理机理和外在影响玉米作物生长的因素。前者主要体现在玉米作物自身的特性上，比如：育种对这种特性的改变；后者主要受外在条件，比如：土壤性质、地表高程、日照、温度、降雨、人类活动等。

如图2所示：具体来说及转移概率a_i,j(t)随时间t变化而改变，我们假设玉米发育阶段在整个生育过程中为单向的，发育阶段S_i仅能转移到其自身，或其下一发育阶段S_i+1。a_i,j(t)可通过公式(1-4)直接利用归一化的地面物候调查数据数据进行计算：

$\left(t\right)=\left\{\begin{array}{c}i,ifi=j=N,{\mathrm{\π}}_N(t-1)\≠0;\\ 1-\underset{k=i+1}{\overset{N}{\Σ}}({\mathrm{\π}}_k\left(t\right)-{\mathrm{\π}}_k(t-1))/{\mathrm{\π}}_i(t-1),\\ ifi=j\≠N,{\mathrm{\π}}_N(t-1)\≠0;\\ \underset{k=i+1}{\overset{N}{\Σ}}({\mathrm{\π}}_k\left(t\right)-{\mathrm{\π}}_k(t-1))/{\mathrm{\π}}_i(t-1),\\ ifi=j-1,{\mathrm{\π}}_N(t-1)\≠0;\\ 0,else\end{array}\right.---(1-4);$

其中，公式中有4种条件限制，前3种条件限制分别对应于最后一发育阶段(收割期)转移到其自身、除了最后一发育阶段的其他阶段转移到其自身、当前发育阶段转移到其下一发育阶段，如果q_t-1＝S₆，则除了a_6,6(t)和a_6,7(t)的其他转移概率均为0，a_6,6(t)和a_6,7(t)分别利用公式(1-4)的第2和3种条件限制进行计算，两者之和为1。

具体为：1)选取同一时间点的多年归一化地面物候调查数据数据，计算其均值；2)利用公式(1-4)计算出每个时间点的状态转移概率矩阵。

观测概率矩阵：观测值随着作物物候的交替连续变化，对于一块研究区域来说，某一时间节点，可能出现多个发育阶段共存的状况。为此，将某一时刻的观测值看成是由多状态观测值的混合。HMM混合模型由若干个服从某分布的分量相互组成，通过混合分布模型来逼近样本的真实分布。混合高斯函数的权重由行政区划单元内各发育阶段玉米作物面积比例决定，即对应于发育阶段先验概率π_i(t)。行政区划单元内观测变量的概率密度函数可由多变量高斯分布函数表示。也就是说，公式(1-3)中，P(O_t)是各子高斯分布函数b_i(O_t)的线性组合，而b_i(O_t)可由参数向量均值μ_i和协方差矩阵Σ_i表示。b_i(O_t)可表示为：

$b_i\left(O_i\right)=N_i\left(O_i\right|{\mathrm{\μ}}_i,{\Σ}_i)=\frac{1}{\sqrt{{\left(2\mathrm{\π}\right)}^d}\left|{\Σ}_i\right|}\·\exp (-\frac{{(O_t-{\mathrm{\μ}}_i)}^{\′}\·{\Σ}_i^{-1}\·(O_t-{\mathrm{\μ}}_i)}{2})---(1-5);$

其中，d对应于特征向量空间的维度，在本案例中因选取了三种特征作为观测值，故d＝3。在本模型中，μ_i和Σ_i为全局HMM模型参数，也就是说该参数与时间向量无关。第i个组份的权重(混合系数)π_j(t)可由地面物候调查数据获得，即仅仅μ_i和Σ_i为未知变量。给定一个观测值序列O₁,…,O_T，我们可以利用最大似然法估计μ_i和Σ_i。假设参数空间Θ＝{μ,Σ}，那么似然函数的log函数形式，可表示为：

$\ln L\left(\mathrm{\Θ}\right|O_t)=\underset{t=1}{\overset{T}{\Σ}}\ln \underset{j=1}{\overset{N}{\Σ}}{\mathrm{\π}}_j\left(t\right)\·N_j\left(O_t\right|{\mathrm{\μ}}_j,{\Σ}_j)---(1-6);$

高斯函数的参数可利用EM算法(期望值最大，Expectation Maximum)进行估计。EM算法主要分为两个步骤：估计步骤E-step和最大化步骤M-step：首先利用样本对参数进行估计，然后在M-step中将需要估计的参数最大化(通常是求其最大似然估计)，不断地迭代此两个步骤，直到收敛，第q+1次迭代的和可表示为：

${\mathrm{\μ}}_j^{q+1}=\frac{{\Σ}_{t=1}^T{O}_t\·{\mathrm{\β}}_j^q\left(t\right)}{{\Σ}_{t=1}^T{\mathrm{\β}}_j^q\left(t\right)}---(1-7);$

${\Σ}_j^{q+1}=\frac{{\Σ}_{t=1}^T{\mathrm{\β}}_j^q\left(t\right)\·(O_t-{\mathrm{\μ}}_j^{q+1})\·{(O_t-{\mathrm{\μ}}_j^{q+1})}^{\′}}{{\Σ}_{t=1}^T{\mathrm{\β}}_j^q\left(t\right)}---(1-8);$

其中：

${\mathrm{\β}}_j^q\left(t\right)=E\left({\mathrm{\π}}_j\left(t\right)\right|O_t;{\mathrm{\Θ}}_j)=\frac{{\mathrm{\π}}_j\left(t\right)\·N\left(O_t\right|{\mathrm{\μ}}_j^q,{\Σ}_j^q)}{{\Σ}_{i=1}^N{\mathrm{\π}}_i\left(t\right)\·N\left(O_t\right|{\mathrm{\μ}}_i^q,{\Σ}_i^q)}---(1-9);$

其中，j＝1,…,N；t＝1,…,T；

对于迭代收敛条件通常可采用如下两种方式：1)|L(X|Θ)-L(X|Θ′)|<ε,其中， $L\left(X\right|\mathrm{\&Theta;})=$ $\log {\&Pi;}_{t=1}^T{\&Sigma;}_{j=1}^N{\mathrm{\&omega;}}_j\left(t\right)\&CenterDot;P(x_t;{\mathrm{\&Theta;}}_j)={\&Sigma;}_{t=1}^T\log {\&Sigma;}_{j=1}^N{\mathrm{\&omega;}}_j\left(t\right)\&CenterDot;P(x_t;{\mathrm{\&Theta;}}_j);$ 2)|Θ-Θ′|<ε。通常，ε＝10^-5。

本实施例中：观测值概率由8个高斯分布函数组成，EM参数估计中我们首先假定 利用公式(1-7)初始化公式(1-8)初始化利用EM算法迭代直至收敛，取最后一次有效迭代作为混合高斯函数的参数。

2)、获得玉米的多源特征，作为HMM混合模型的输入部分。

多源特征包括NDVI均值、分维值、和有效积温这些全局特征组成。其中：

NDVI均值我们计算其均值。针对每日的MODIS NDVI(采用MODIS获得的植被覆盖指数)时间序列影像，包括：1)影像合成，即利用MVC(Maximum Value Composting)将每日的MODIS NDVI时间序列影像合成为每周的MODIS NDVI时间序列影像；2)影像掩膜，即利用NASS的CDL数据作为掩膜板，剔除NDVI影像上非玉米作物类型的像元。

分维值并作为影像纹理异质性的指标来反映玉米作物生长状况，不同玉米田块处于不同的物候期，故其粗糙度的变化也伴随整个生长周期，即分维时间序列反映了玉米作物的发育阶段在整个生长周期中的时空变化，分维值利用DR-DBC算法进行估计，该计算同样作用于掩膜处理后的MODIS NDVI影像。

对作物生长发育起重要作用的两个常规环境因子包括光照和温度，而现今杂交玉米更加偏重于温度的影响，而较弱的受制于光照因素，作物生长发育通常需要一定的温度(热量)条件，在作物生长发育所需要的其他条件均得到满足时，在一定温度范围内，气温和发育速度成正相关，并且要积累到一定的温度总和，才能完成其发育期，这个温度的累积数称为积温。每种作物都有其生长的下限温度。当温度高于下限温度时，它才能生长发育。这个对作物生长发育起有效作用的高出的温度值，称作有效温度。作物在整个生育期内有效温度的累积，即有效积温。为此，我们给出有效积温的计算公式：

$AGDDs\left(t\right)=\underset{t}{\&Sigma;}(\frac{T_{\max }\left(t\right)-T_{\min }\left(t\right)}{2}-T_{base})---(1-10);$

其中，T_base表示基准温度，或下限温度，T_max(t)和T_min(t)表示校正后的日最高和最低温度。可以从地面气象站获取地表日最低和最高温度。

其中：

t：积日(Day of Year，DOY)，对应一年的第t日。

和单站日最低和最高温度，对应第i个气象站、第t个积日的观测。

和行政区划级别日最低和最高温度，对于某行政区划内、第t个积日的温度。

T_max(t)和T_min(t)：校正后的行政区划级别日最高和最低温度。

然而,从气象站获取的是离散点位信息。要将离散点位观测转换成行政区划级别的有效积温，需要经过以下两步操作：1)将离散点位观测和转换为行政区划级别面状温度数据和2)根据行政区划级别的日最低和最高温度，和日温度校正规则，根据公式(1-10)计算出有效积温。

具体为先利用行政区划内的气象站观测，通过泰森多边形加权法，计算出和每个气象站的面积权重系数w_i＝A_i/A_total，其中，A_total为某行政区划的面积，A_i为该行政区划受第i个气象站影响的面积，以的计算为例的计算方法相似)，给出如下公式：

${\stackrel{\&OverBar;}{T}}_{\min }\left(t\right)={\&Sigma;}_{i=1}^n{w}_i\&CenterDot;T_{\min }^i\left(t\right)---(1-11);$

其中，n表示行政区划内有效的气象站数量，根据日温度校正规则和公式(1-10)计算出有效积温。

3)、通过HMM混合模型计算出农作物各发育阶段的百分比。

在给定HMM模型参数λ、观测值序列O₁,…,O_t，求隐含状态j的后验概率P(q_t＝S_j|O₁,…,O_t)。在HMM中，已知模型参数和某一特定观测值序列，求解当前状态的概率分布，即P(q_t＝S_j|O₁,…,O_t)，被称为过滤(filtering)问题，使用前向法(forward algorithm)解决。现将P(q_t＝S_j,O₁,…,O_t)简写为κ_j(t)，用于表示给定观测值序列O₁,…,O_t，t时刻隐含S_j发生的概率，则：

$P(q_t=S_j|O_1,\mathrm{...},O_t)=\frac{{\mathrm{\&kappa;}}_j\left(t\right)}{{\&Sigma;}_{i=1}^N{\mathrm{\&kappa;}}_i\left(t\right)}---(1-12);$

下式中κ_j(1)由初始时刻先验π_j(已知)和P(O₁|q₁＝S_j)(可由计算出来的HMM参数求得)。同理，κ_j(2)的计算，只需要多考虑一个转换概率，而转换概率是可由多年统计数据计算得到。

π_j·P(O₁|q₁＝S_j)＝κ_j(1)             (1-13)；

$(\underset{i=1}{\overset{N}{\&Sigma;}}{\mathrm{\&kappa;}}_i\left(1\right)\&CenterDot;P(q_2=S_j|q_1=S_i))\&CenterDot;P\left(O_2\right|q_2=S_j)={\mathrm{\&kappa;}}_j\left(2\right)---(1-14);$

……

$(\underset{i=1}{\overset{N}{\&Sigma;}}{\mathrm{\&kappa;}}_i(t-1)\&CenterDot;P(q_t=S_j|q_{t-1}=S_i))\&CenterDot;P\left(O_t\right|q_t=S_j)={\mathrm{\&kappa;}}_j\left(t\right)---(1-15);$

由于P(q_t＝S_j|O₁,…,O_t)得到的是各发育阶段的百分比，对应于归一化后的地面物候调查数据数据。为了反求原始的地面物候调查数据数值，我们可采用下式计算求得：

${\mathrm{\&delta;}}_i\left(t\right)=\underset{k=1}{\overset{N}{\&Sigma;}}P(q_t=S_k|O_1,\mathrm{...},O_t)---(1-16);$

其中，i＝2,…,N。

上述实施例只为说明本发明的技术构思及特点，其目的在于让熟悉此项技术的人士能够了解本发明的内容并据以实施，并不能以此限制本发明的保护范围。凡根据本发明精神实质所作的等效变化或修饰，都应涵盖在本发明的保护范围之内。

Claims (6)

1. 一种基于HMM的农作物物候动态估计方法，其特征在于：包括：

1）：建立HMM混合模型，HMM混合模型中包括获得：

a、初始概率：农作物在各发育阶段发生的概率；

b、局部状态转移矩阵：农作物在整个生育周期过程中的各发育阶段转移概率；

c、观测概率矩阵：在给定发育阶段状态的情况下，各观测变量发生的概率；

2）：向HMM混合模型输入农作物的多源特征，多源特征包括：

a、NDVI均值；

b、分维值；

c、有效积温；

3）、通过HMM混合模型计算出农作物各发育阶段的百分比。

2. 根据权利要求1所述的一种基于HMM的农作物物候动态估计方法，其特征在于：初始概率通过统计该地区同一时期的历史数据获得。

3. 根据权利要求1所述的一种基于HMM的农作物物候动态估计方法，其特征在于：获取局部状态转移矩阵包括：

1）、选取同一时间点的地面物候调查数据，计算其均值；

2）、计算出每个时间点的状态转移概率矩阵。

4. 根据权利要求1所述的一种基于HMM的农作物物候动态估计方法，其特征在于：获取NDVI均值包括：

1）、影像合成：将每日的MODIS NDVI时间序列影像合成为每周的NDVI时间序列影像；

2）、影像掩膜：剔除MODIS NDVI时间序列影像上非农作物类型的像元。

5. 根据权利要求4所述的一种基于HMM的农作物物候动态估计方法，其特征在于：对剔除MODIS NDVI时间序列影像上非农作物类型的像元进行估算获得。

6. 根据权利要求1所述的一种基于HMM的农作物物候动态估计方法，其特征在于：有效积温获取包括：获取离散点位的地表日最低和最高温度；将离散点位的地表日最低和最高温度转化为行政区划级别面状温度数据；对行政区划级别的日最低和最高温度进行日温校正，并计算出有效积温。