气候变化预测以及趋势分析方法

摘要

本发明实施例提供一种气候变化预测以及趋势分析方法。本发明方法，包括：处理器接收气候变量数据以及与所述气候变量数据对应的时间值；所述处理器采用模糊聚类算法计算所述气候变量数据的渐进变化与跳跃变化；所述处理器根据所述渐进变化与所述跳跃变化对气候变化进行预测以及趋势分析。本发明实施例准确体现气候变化中跳跃变化，有利于进一步准确预测资源的合理利用及区域可持续发展，有助于分析对应该跳跃变化的成因。

说明书

技术领域

本发明实施例涉及计算机计算领域，尤其涉及一种气候变化预测以及趋势分析方法。

背景技术

近年来，洪涝、干旱和冰雨等极端气候引起了全世界各地人们的关注。造成该现象的原因有很多。例如，冷暖气流的急剧变化等自然因素，同时，人类对森林的滥砍滥伐、汽车的尾气排放等人类活动的人为因素也占很大的比例。

现有技术对此采用混沌-人工神经网络方法对气候变化进行预测以及趋势分析，该方法是结合混沌理论和人工神经网络，互相抵消两者之间的不足之处。分析流域气温、降雨和径流量的变化过程。进而对于气候的变化预测以及趋势分析。

该种方法对于气候的变化体现了气候变化的整体趋势，但是无法准确地体现气候跳跃变化。

发明内容

本发明实施例提供一种气候变化预测以及趋势分析方法，以克服现有中无法准确地对于气候变化中的跳跃变化作出预测以及趋势分析的问题。

本实施例提供了一种气候变化预测以及趋势分析方法，包括：

处理器接收气候变量数据以及与所述气候变量数据对应的时间值；

所述处理器采用模糊聚类算法计算所述气候变量数据的渐进变化与跳跃变化；

所述处理器根据所述渐进变化与所述跳跃变化对气候变化进行预测以及趋势分析。

进一步地，所述处理器采用模糊聚类算法计算所述气候变量数据的渐进变化与跳跃变化，包括：

步骤一、处理器初始化所述模糊聚类算法参量，设置所述模糊聚类算法终止阈值、初始分割阶数和最小分割阶数；

步骤二、所述处理器求得隶属度、混合概率，并判断所述混合概率是否大于0，若是，则直接更新聚类参数以及所述隶属度，若否，则降低分割阶数后更新聚类参数以及所述隶属度；

步骤三、所述处理器判断相邻两次迭代过程的模糊划分矩阵的相对误差是否大于终止阈值，若否，则返回步骤二，若是，则根据最小信息长度准则选择最优分割阶数；

步骤四、所述处理器将所述最优分割阶数与最优聚类参数代入模糊聚类算法得出所述气候变量数据的渐进变化与跳跃变化。

进一步地，所述处理器初始化所述模糊聚类算法参量，包括：

处理器采用分段常函数得到在时间坐标集合上k_max个初始模糊聚类中心为

$>{\mu }_k=\frac{{\Sigma }_{i=1}^N{\left[s_k\left(t_i\right)\right]}^m{t}_i}{{\Sigma }_{i=1}^N{\left[s_k\left(t_i\right)\right]}^m},1\le k\le k_{\max }---\left(1\right)$>

其中，N为气候变量数据的个数，m为模糊加权指数且m＝2，t_i为气候变量数据对应的时间，p＝int(N/k_max)，当k＝k_max，令kp＝N；

处理器将时间坐标集合上的初始偏差定义为

$>{\sigma }_k^2=\frac{{\Sigma }_{i=1}^N{\left[s_k\left(t_i\right)\right]}^m(t_i-{\mu }_k){(t_i-{\mu }_k)}^T}{{\Sigma }_{i=1}^N{\left[s_k\left(t_i\right)\right]}^m},1\le k\le k_{\max }---\left(2\right)$>

其中，N为气候变量数据的个数，m为模糊加权指数且m＝2，t_i为气候变量数据对应的时间，μ_k为在时间坐标集合上的初始模糊聚类中心，s_k(t_i)为如上所述分段常函数；

所述处理器从气候变量数据集合中随机选取k_max个数据作为初始中心，得到

v_k∈{x_i,1≤i≤N},1≤k≤k_max   (3)

其中，N为气候变量数据的个数，x_i为气候变量数据；

所述处理器令所述气候变量数据集合中的初始协方差与单位矩阵成正比，得到

$>{\Sigma }_k={\sigma }^2{I}_{k_{\max }\times k_{\max }}---\left(4\right)$>

其中，为k_max阶单位矩阵，$>{\sigma }^2=\frac{1}{10}\mathrm{trace}\left(\frac{1}{N}\underset{i=1}{\overset{N}{\Sigma }}(x_i-m){(x_i-m)}^T\right),$>为气候变量数据的平均值；

所述处理器令初始混合概率α_k与所述k_max成反比，得到

α_k＝1/k_max,1≤k≤k_max   (5)

其中，k_max为初始聚类数。

进一步地，所述处理器将所述最优分割阶数代入模糊聚类算法得出所述气候变量数据的渐进变化与跳跃变化，包括：

处理器采用高斯函数得到在时间坐标集合上k_max个模糊分割β_k(t_i)，公式为

$>{\beta }_k\left(t_i\right)=\frac{A_k\left(t_i\right)}{{\Sigma }_{j=1}^{k_{\max }}{A}_j\left(t_i\right)},1\le k\le k_{\max }---\left(6\right)$>

其中t_i为气候变量数据对应的时间，μ_k为在时间坐标集合上的初始模糊聚类中心，μ_k为时间坐标集合上的初始偏差；

进一步地，所述处理器求得隶属度、混合概率,包括：

采用公式

$>u_{i,k}^{\left(l\right)}=\frac{{[{\alpha }_k^{\left(l\right)}p\left(t_i\right|{\eta }_k^{\left(l\right)})p\left(x_i\right|{\eta }_k^{\left(l\right)})]}^{m-1}}{{\Sigma }_{j=1}^{k_{\max }}{\left[{\alpha }_j^{\left(l\right)}p\left(t_i\right|{\eta }_j^{\left(l\right)})p\left(x_i\right|{\eta }_j^{\left(l\right)})\right]}^{m-1}}---\left(7\right)$>

计算隶属度，其中，所述为第k个聚类的混合概率，p(t|η_k)为在时间坐标集合上的第k个高斯函数，p(x|η_k)为在气候变量数据集合上第k个高斯函数，m为模糊加权指数且m＝2；

采用公式

$>{\alpha }_k^{(l+1)}=\frac{\max \{0,\left(\underset{i=1}{\overset{N}{\Sigma }}{u}_{i,k}^{\left(l\right)}\right)-\frac{c}{2}\}}{\underset{j=1}{\overset{k_{\max }}{\Sigma }}\max \{0,\left(\underset{i=1}{\overset{N}{\Sigma }}{u}_{i,k}^{\left(l\right)}\right)-\frac{c}{2}\}},\{{\alpha }_1^{(l+1)},...,{\alpha }_{k_{\max }}^{(l+1)}\}=\frac{\{{\alpha }_1^{(l+1)},...,{\alpha }_{k_{\max }}^{(l+1)}\}}{{\Sigma }_{k=1}^{k_{\max }}{\alpha }_k^{(l+1)}}---\left(8\right)$>

计算混合概率，其中，所述u_i,k为第i个数据属于第k个聚类的隶属度，c为每个聚类成员所含参数的维数；

所述则更新聚类参数以及所述隶属度，包括：

处理器采用期望值最大算法将所述聚类参数以及隶属度更新为

$>{\mu }_k^{(l+1)}=\frac{{\Sigma }_{i=1}^N{\left(u_{i,k}^{\left(l\right)}\right)}^m{t}_i}{{\Sigma }_{i=1}^N{\left(u_{i,k}^{\left(l\right)}\right)}^m}---\left(9\right)$>

$>{\left({\sigma }_k^2\right)}^{(l+1)}=\frac{{\Sigma }_{i=1}^N{\left(u_{i,k}^{\left(l\right)}\right)}^m(t_i-{\mu }_k^{\left(l\right)}){(t_i-{\mu }_k^{\left(l\right)})}^T}{{\Sigma }_{i=1}^N{\left(u_{i,k}^{\left(l\right)}\right)}^m}---\left(10\right)$>

$>v_k^{(l+1)}=\frac{{\Sigma }_{i=1}^N{\left(u_{i,k}^{\left(l\right)}\right)}^m{x}_i}{{\Sigma }_{i=1}^N{\left(u_{i,k}^{\left(l\right)}\right)}^m}---\left(11\right)$>

$>{\Sigma }_k^{(l+1)}=\frac{{\Sigma }_{i=1}^N{\left(u_{i,k}^{\left(l\right)}\right)}^m(x_i-v_k^{\left(l\right)}){(x_i-v_k^{\left(l\right)})}^T}{{\Sigma }_{i=1}^N{\left(u_{i,k}^{\left(l\right)}\right)}^m}---\left(12\right)$>

$>u_{i,k}^{(l+1)}=\frac{{\left[{\alpha }_k^{(l+1)}p\left(t_i\right|{\eta }_k^{(l+1)})p\left(x_i\right|{\eta }_k^{(l+1)})\right]}^{m-1}}{{\Sigma }_{j=1}^{k_{\max }}{\left[{\alpha }_j^{(l+1)}p\left(t_i\right|{\eta }_j^{(l+1)})p\left(x_i\right|{\eta }_j^{(l+1)})\right]}^{m-1}}---\left(13\right)$>

其中，所述u_i,k为第i个数据属于第k个聚类的隶属度，x_i为气候变量数据，t_i为气候变量数据x_i对应的时间，μ_k与分别为在时间坐标集合上的第k个高斯函数p(t|η_k)的中心和偏差，v_k与Σ_k分别为在气候变量数据集合上第k个高斯函数p(x|η_k)的中心和方差，m为模糊加权指数且m＝2。

进一步地，所述根据最小信息长度准则选择最优分割阶数，包括：

处理器根据最小信息长度准则，采用公式

$>\left(\begin{array}{c}\mathrm{MessLen}({\eta }^{(l+1)},Z)=-\underset{i=1}{\overset{N}{\Sigma }}\log \underset{k:{\alpha }_k>0}{\Sigma }{\alpha }_{k}p\left(t_i\right|{\eta }_k^{(l+1)})p\left(x_i\right|{\eta }_k^{(l+1)})\\ +\frac{c}{4}\underset{k:{\alpha }_k>0}{\Sigma }\log \left(\frac{N{\alpha }_k}{12}\right)+\frac{k_{\mathrm{nz}}}{2}\log \frac{N}{12}+\frac{k_{\mathrm{nz}}(c+1)}{2}.\end{array}\right)---\left(14\right)$>

求得最小信息长度，其中，所述c为每个聚类成员所含参数的个数，α_k为第k个聚类的混合概率，N为气候变量数据的个数，k_nz为当前聚类个数，p(t|η_k)为在时间坐标集合上的第k个高斯函数，p(x|η_k)为在气候变量数据集合上第k个高斯函数。

进一步地，所述步骤三之后，还包括：

处理器判断若相邻两次迭代过程的模糊划分矩阵的相对误差低于终止阈值时，删除具有最低概率的聚类成员并返回所述步骤二，查找更小的信息长度值。

本发明实施例处理器接收气候变量数据以及与所述气候变量数据对应的时间值，采用模糊聚类算法计算所述气候变量数据的渐进变化与跳跃变化，并根据所述渐进变化与所述跳跃变化对齐气候变化进行预测以及趋势分析，准确地体现了气候变化中的跳跃变化，有利于进一步准确预测资源的合理利用及区域可持续发展，有助于分析对应该跳跃变化的成因。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作一简单地介绍，显而易见地，下面描述中的附图是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图1为本发明气候变化预测以及趋势分析方法流程图；

图2为本发明模糊聚类算法的流程图；

图3为本发明仿真数据示意图；

图4为本发明仿真数据的初始模糊分割示意图；

图5为本发明仿真数据的最小信息长度准则的演化示意图；

图6为本发明仿真数据的最优模糊分割结果示意图；

图7为本发明仿真数据的最优隶属度结果示意图；

图8为本发明仿真数据的最优明确分割结果示意图。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

图1为本发明气候变化预测以及趋势分析方法流程图，如图1所示，本实施例方法，包括：

步骤101、处理器接收气候变量数据以及与所述气候变量数据对应的时间值；

步骤102、所述处理器采用模糊聚类算法计算所述气候变量数据的渐进变化与跳跃变化；

步骤103、所述处理器根据所述渐进变化与所述跳跃变化对气候变化进行预测以及趋势分析。

具体来说，处理器接收气候变量数据，同时还接收与该气候变量数据相对应的时间值。例如，气候变化数据为花园口水温站提供的径流数据，并且该径流数据还对应标明年份的时间值。处理器根据该气候变量数据以及时间值，采用模糊聚类算法计算气候变量数据的渐进变化与跳跃变化。从而预测气候变化与分析气候变化的趋势。

假设数据可以被具有参数η＝{η₁,…,η_K}的多元(包括时间坐标)高斯混合分布有效建模，其中x_i＝[x_1,i,x_2,i,…,x_d,i]^T(1≤i≤N),t_i是数据x_i对应的时间点，K是聚类的数目，η_k是第k(1≤k≤K)个聚类参数构成的向量。

因为时间变量t_i独立于变量x_i，所以数据z_i到第k个聚类原型的距离测度D²(z_i,η_k)由下式表示：

$>\frac{1}{D^2(z_i,{\eta }_k)}=p\left({\eta }_k\right)p\left(z_i\right|{\eta }_k)=p\left({\eta }_k\right)p\left(t_i\right|{\eta }_k)p\left(x_i\right|{\eta }_k)---\left(1\right)$>

其中η_k表示第k个聚类原型的参数，表示第k个聚类的混合概率。p(t_i|η_k)表示在时间坐标集合上时间t_i属于第k个聚类的概率，且

$>p\left(t_i\right|{\eta }_k)=\frac{1}{\sqrt{2\pi {\sigma }_k^2}}\exp (-\frac{1}{2}\frac{{(t_i-{\mu }_k)}^2}{{\sigma }_k^2})---\left(2\right)$>

其中μ_k与分别表示在时间坐标集合上的第k个高斯函数p(t|η_k)的中心和方差,p(x_i|η_k)表示在数据集合上数据x_i属于第k个聚类的概率,且

$>p\left(x_i\right|{\eta }_k)=\frac{1}{{\left(2\pi \right)}^{\frac{d}{2}}\sqrt{\det \left({\Sigma }_k\right)}}\exp (-\frac{1}{2}{(x_i-v_k)}^T{\Sigma }_k^{-1}(x_i-v_k))---\left(3\right)$>

其中v_k与Σ_k分别表示在数据集合上第k个高斯函数p(x|η_k)的中心和方差。第k个聚类的原型参数K个聚类的所有参数的集合定义为η≡{η₁,…,η_K}。此时,概率密度函数定义为：

$>p\left(Z\right|\eta )=\underset{i=1}{\overset{N}{\Pi }}p\left(z_i\right|\eta )=\underset{i=1}{\overset{N}{\Pi }}\underset{k=1}{\overset{K}{\Sigma }}{\alpha }_{k}p\left(t_i\right|{\eta }_k)p\left(x_i\right|{\eta }_k)---\left(4\right)$>

在时间坐标集合上时间点t_i属于第k个聚类的隶属度β_k(t_i)∈[0,1]表示时间点t_i属于第k个时间序列片段区间的程度,即

$>A_k\left(t_i\right)=\exp (-\frac{1}{2}\frac{{(t_i-{\mu }_k)}^2}{{\sigma }_k^2}),{\beta }_k\left(t_i\right)=\frac{A_k\left(t_i\right)}{{\Sigma }_{j=1}^K{A}_j\left(t_i\right)}---\left(5\right)$>

本发明实施例处理器接收气候变量数据以及与所述气候变量数据对应的时间值，采用模糊聚类算法计算所述气候变量数据的渐进变化与跳跃变化，并根据所述渐进变化与所述跳跃变化对齐气候变化进行预测以及趋势分析，准确地体现了气候变化中的跳跃变化，有利于进一步准确预测资源的合理利用及区域可持续发展，有助于分析对应该跳跃变化的成因。

图2为本发明模糊聚类算法的流程图，如图2所示，本实施例方法，包括：

步骤201、处理器初始化所述模糊聚类算法参量，设置所述模糊聚类算法终止阈值、初始分割阶数和最小分割阶数；

具体来说，处理器采用分段常函数得到在时间坐标集合上k_max个初始模糊聚类中心为

$>{\mu }_k=\frac{{\Sigma }_{i=1}^N{\left[s_k\left(t_i\right)\right]}^m{t}_i}{{\Sigma }_{i=1}^N{\left[s_k\left(t_i\right)\right]}^m},1\le k\le k_{\max }---\left(6\right)$>

其中，N为气候变量数据的个数，m为模糊加权指数且m＝2，t_i为气候变量数据对应的时间，p＝int(N/k_max)，当k＝k_max，令kp＝N；

处理器将时间坐标集合上的初始偏差定义为

$>{\sigma }_k^2=\frac{{\Sigma }_{i=1}^N{\left[s_k\left(t_i\right)\right]}^m(t_i-{\mu }_k){(t_i-{\mu }_k)}^T}{{\Sigma }_{i=1}^N{\left[s_k\left(t_i\right)\right]}^m},1\le k\le k_{\max }---\left(7\right)$>

其中，N为气候变量数据的个数，m为模糊加权指数且m＝2，t_i为气候变量数据对应的时间，μ_k为在时间坐标集合上的初始模糊聚类中心，s_k(t_i)为如上所述分段常函数；

所述处理器从气候变量数据集合中随机选取k_max个数据作为初始中心，得到

v_k∈{x_i,1≤i≤N},1≤k≤k_max   (8)

其中，N为气候变量数据的个数，x_i为气候变量数据；

所述处理器令所述气候变量数据集合中的初始协方差与单位矩阵成正比，得到

$>{\Sigma }_k={\sigma }^2{I}_{k_{\max }\times k_{\max }}---\left(9\right)$>

其中，为k_max阶单位矩阵，$>{\sigma }^2=\frac{1}{10}\mathrm{trace}\left(\frac{1}{N}\underset{i=1}{\overset{N}{\Sigma }}(x_i-m){(x_i-m)}^T\right),$>为气候变量数据的平均值；

所述处理器令初始混合概率α_k与所述k_max成反比，得到

α_k＝1/k_max,1≤k≤k_max   (10)

其中，k_max为初始聚类数。

步骤202、所述处理器求得隶属度、混合概率，并判断所述混合概率是否大于0，若是，则直接更新聚类参数以及所述隶属度，若否，则降低分割阶数后更新聚类参数以及所述隶属度；

具体来说，采用公式

$>u_{i,k}^{\left(l\right)}=\frac{{[{\alpha }_k^{\left(l\right)}p\left(t_i\right|{\eta }_k^{\left(l\right)})p\left(x_i\right|{\eta }_k^{\left(l\right)})]}^{m-1}}{{\Sigma }_{j=1}^{k_{\max }}{\left[{\alpha }_j^{\left(l\right)}p\left(t_i\right|{\eta }_j^{\left(l\right)})p\left(x_i\right|{\eta }_j^{\left(l\right)})\right]}^{m-1}}---\left(11\right)$>

计算隶属度，其中，所述为第k个聚类的混合概率，p(t|η_k)为在时间坐标集合上的第k个高斯函数，p(x|η_k)为在气候变量数据集合上第k个高斯函数，m为模糊加权指数且m＝2；

采用公式

$>{\alpha }_k^{(l+1)}=\frac{\max \{0,\left(\underset{i=1}{\overset{N}{\Sigma }}{u}_{i,k}^{\left(l\right)}\right)-\frac{c}{2}\}}{\underset{j=1}{\overset{k_{\max }}{\Sigma }}\max \{0,\left(\underset{i=1}{\overset{N}{\Sigma }}{u}_{i,k}^{\left(l\right)}\right)-\frac{c}{2}\}},\{{\alpha }_1^{(l+1)},...,{\alpha }_{k_{\max }}^{(l+1)}\}=\frac{\{{\alpha }_1^{(l+1)},...,{\alpha }_{k_{\max }}^{(l+1)}\}}{{\Sigma }_{k=1}^{k_{\max }}{\alpha }_k^{(l+1)}}---\left(12\right)$>

计算混合概率，其中，所述u_i,k为第i个数据属于第k个聚类的隶属度，c为每个聚类成员所含参数的维数；

所述则更新聚类参数以及所述隶属度，包括：

处理器采用期望值最大算法将所述聚类参数以及隶属度更新为

$>{\mu }_k^{(l+1)}=\frac{{\Sigma }_{i=1}^N{\left(u_{i,k}^{\left(l\right)}\right)}^m{t}_i}{{\Sigma }_{i=1}^N{\left(u_{i,k}^{\left(l\right)}\right)}^m}---\left(13\right)$>

$>{\left({\sigma }_k^2\right)}^{(l+1)}=\frac{{\Sigma }_{i=1}^N{\left(u_{i,k}^{\left(l\right)}\right)}^m(t_i-{\mu }_k^{\left(l\right)}){(t_i-{\mu }_k^{\left(l\right)})}^T}{{\Sigma }_{i=1}^N{\left(u_{i,k}^{\left(l\right)}\right)}^m}---\left(14\right)$>

$>v_k^{(l+1)}=\frac{{\Sigma }_{i=1}^N{\left(u_{i,k}^{\left(l\right)}\right)}^m{x}_i}{{\Sigma }_{i=1}^N{\left(u_{i,k}^{\left(l\right)}\right)}^m}---\left(15\right)$>

$>{\Sigma }_k^{(l+1)}=\frac{{\Sigma }_{i=1}^N{\left(u_{i,k}^{\left(l\right)}\right)}^m(x_i-v_k^{\left(l\right)}){(x_i-v_k^{\left(l\right)})}^T}{{\Sigma }_{i=1}^N{\left(u_{i,k}^{\left(l\right)}\right)}^m}---\left(16\right)$>

$>u_{i,k}^{(l+1)}=\frac{{\left[{\alpha }_k^{(l+1)}p\left(t_i\right|{\eta }_k^{(l+1)})p\left(x_i\right|{\eta }_k^{(l+1)})\right]}^{m-1}}{{\Sigma }_{j=1}^{k_{\max }}{\left[{\alpha }_j^{(l+1)}p\left(t_i\right|{\eta }_j^{(l+1)})p\left(x_i\right|{\eta }_j^{(l+1)})\right]}^{m-1}}---\left(17\right)$>

其中，所述u_i,k为第i个数据属于第k个聚类的隶属度，x_i为气候变量数据，t_i为气候变量数据x_i对应的时间，μ_k与分别为在时间坐标集合上的第k个高斯函数p(t|η_k)的中心和偏差，v_k与Σ_k分别为在气候变量数据集合上第k个高斯函数p(x|η_k)的中心和方差，m为模糊加权指数且m＝2。

步骤203、所述处理器判断相邻两次迭代过程的模糊划分矩阵的相对误差是否大于终止阈值，若否，则返回步骤二，若是，则根据最小信息长度准则选择最优分割阶数；

步骤204、所述处理器将所述最优分割阶数与最优聚类参数代入模糊聚类算法得出所述气候变量数据的渐进变化与跳跃变化。

具体来说，基于模糊聚类的最小信息长度准则：

$>\hat{\eta }=\arg \underset{\eta }{\min }\mathrm{MessLen}(\eta ,Z)---\left(18\right)$>

且

其中c是每个聚类成员所含参数的个数(即η_k的维数)，N是数据集合Z所含数据的总数。当α_k等于零,MML的惩罚项将趋于-∞。通过只编码概率不为零的聚类成员参数来解决这个问题。令表示具有非零概率的聚类成员的数目,则

$>\mathrm{MessLen}(\eta ,Z)=-\underset{i=1}{\overset{N}{\Sigma }}\log \underset{k:{\alpha }_k>0}{\Sigma }{\alpha }_{k}p\left(t_i\right|{\eta }_k)p\left(x_i\right|{\eta }_k)+\frac{c}{4}\underset{k:{\alpha }_k>0}{\Sigma }\log \left(\frac{N{\alpha }_k}{12}\right)+\frac{K_{\mathrm{nz}}}{2}\log \frac{N}{12}+\frac{K_{\mathrm{nz}}(c+1)}{2}---\left(20\right)$>

所以,基于模糊聚类的MML准则为

$>\hat{\eta }=\arg \underset{\eta }{\min }\{-\underset{i=1}{\overset{N}{\Sigma }}\log \underset{k:{\alpha }_k>0}{\Sigma }{\alpha }_{k}p\left(t_i\right|{\eta }_k)p\left(x_i\right|{\eta }_k)+\frac{c}{4}\underset{k:{\alpha }_k>0}{\Sigma }\log \left(\frac{N{\alpha }_k}{12}\right)+\frac{K_{\mathrm{nz}}}{2}\log \frac{N}{12}+\frac{K_{\mathrm{nz}}(c+1)}{2}\}---\left(21\right)$>

模糊聚类删减策略具体为：

EM算法最小化代价函数有下面的M步骤：

$>{\alpha }_k=\frac{\max \{0,\left(\underset{i=1}{\overset{N}{\Sigma }}{u}_{i,k}\right)-\frac{c}{2}\}}{\underset{j=1}{\overset{K}{\Sigma }}\max \{0,\left(\underset{i=1}{\overset{N}{\Sigma }}{u}_{i,k}\right)-\frac{c}{2}\}},k=\mathrm{1,2},...,K,---\left(22\right)$>

对于k:α_k>0,更新参数可得

$>{\mu }_k=\frac{{\Sigma }_{i=1}^N{u_{i,k}^{m}t}_i}{{\Sigma }_{i=1}^N{u}_{i,k}^m}---\left(23\right)$>

$>{\sigma }_k^2=\frac{{\Sigma }_{i=1}^N{u}_{i,k}^m(t_i-{\mu }_k){(t_i-{\mu }_k)}^T}{{\Sigma }_{i=1}^N{u}_{i,k}^m}---\left(24\right)$>

$>v^k=\frac{{\Sigma }_{i=1}^N{u_{i,k}^{m}x}_i}{{\Sigma }_{i=1}^N{u}_{i,k}^m}---\left(25\right)$>

$>{\Sigma }_k=\frac{{\Sigma }_{i=1}^N{u}_{i,k}^m(x_i-v_k){(x_i-v_k)}^T}{{\Sigma }_{i=1}^N{u}_{i,k}^m}---\left(26\right)$>

其中u_i,k在E步骤的更新公式如下

$>u_{i,k}=\frac{1}{{\Sigma }_{j=1}^{k_{\max }}{[D(z_i,{\eta }_k)/D(z_i,{\eta }_j)]}^{2/(m-1)}}=\frac{{\left[{\alpha }_{k}p\left(t_i\right|{\eta }_k)p\left(x_i\right|{\eta }_k)\right]}^{m-1}}{{\Sigma }_{j=1}^{k_{\max }}{\left[{\alpha }_{j}p\left(t_i\right|{\eta }_j)p\left(x_i\right|{\eta }_j)\right]}^{m-1}}---\left(27\right)$>

在删除运算中采用针对混合模型的分量形式的EM算法(CEM²)。这个算法并不同时更新所有的聚类成员参数,而是序贯更新聚类成员的参数：更新α₁和η₁,计算

更新α₂和η₂,重新计算等等。当删除一个聚类成员时(α_k(l+1)＝0),CEM²算法重新分配概率密度来增加其余聚类成员的生存机会。这样初始聚类的数目可以任意大,而不产生任何问题。CEM²算法的计算负担仅稍重于EM算法。

图3为本发明仿真数据示意图，图3所示为黄河流域花园口水文站1919-2004年年径流数据。

图4为本发明仿真数据的初始模糊分割示意图。

具体来说，处理器采用分段常函数得到在时间坐标集合上k_max个初始模糊聚类中心为

$>{\mu }_k=\frac{{\Sigma }_{i=1}^N{\left[s_k\left(t_i\right)\right]}^m{t}_i}{{\Sigma }_{i=1}^N{\left[s_k\left(t_i\right)\right]}^m},1\le k\le k_{\max }---\left(28\right)$>

其中，N为气候变量数据的个数，m为模糊加权指数且m＝2，t_i为气候变量数据对应的时间，p＝int(N/k_max)，当k＝k_max，令kp＝N；

处理器将时间坐标集合上的初始偏差定义为

$>{\sigma }_k^2=\frac{{\Sigma }_{i=1}^N{\left[s_k\left(t_i\right)\right]}^m(t_i-{\mu }_k){(t_i-{\mu }_k)}^T}{{\Sigma }_{i=1}^N{\left[s_k\left(t_i\right)\right]}^m},1\le k\le k_{\max }---\left(29\right)$>

其中，N为气候变量数据的个数，m为模糊加权指数且m＝2，t_i为气候变量数据对应的时间，μ_k为在时间坐标集合上的初始模糊聚类中心，s_k(t_i)为如上所述分段常函数；

所述处理器从气候变量数据集合中随机选取k_max个数据作为初始中心，得到

v_k∈{x_i,1≤i≤N},1≤k≤k_max   (30)

其中，N为气候变量数据的个数，x_i为气候变量数据；

所述处理器令所述气候变量数据集合中的初始协方差与单位矩阵成正比，得到

$>{\Sigma }_k={\sigma }^2{I}_{k_{\max }\times k_{\max }}---\left(31\right)$>

其中，为k_max阶单位矩阵，$>{\sigma }^2=\frac{1}{10}\mathrm{trace}\left(\frac{1}{N}\underset{i=1}{\overset{N}{\Sigma }}(x_i-m){(x_i-m)}^T\right),$>为气候变量数据的平均值；

所述处理器令初始混合概率α_k与所述k_max成反比，得到

α_k＝1/k_max,1≤k≤k_max   (32)

其中，k_max为初始聚类数。

进一步地，所述处理器将所述初始模糊聚类参数代入模糊聚类算法得出所述气候变量数据的初始模糊分割，包括：

处理器采用高斯函数得到在时间坐标集合上k_max个模糊分割β_k(t_i)，公式为

$>{\beta }_k\left(t_i\right)=\frac{A_k\left(t_i\right)}{{\Sigma }_{j=1}^{k_{\max }}{A}_j\left(t_i\right)},1\le k\le k_{\max }---\left(33\right)$>

其中t_i为气候变量数据对应的时间，μ_k为在时间坐标集合上的初始模糊聚类中心，为时间坐标集合上的初始偏差。

图5为本发明仿真数据的最小信息长度准则的演化示意图。

具体来说，基于模糊聚类的最小信息长度准则：

$>\hat{\eta }=\arg \underset{\eta }{\min }\mathrm{MessLen}(\eta ,Z)---\left(34\right)$>

且

其中c是每个聚类成员所含参数的个数(即η_k的维数),N是数据集合Z所含数据的总数。当α_k等于零,MML的惩罚项将趋于-∞。通过只编码概率不为零的聚类成员参数来解决这个问题。令表示具有非零概率的聚类成员的数目,则

$>\mathrm{MessLen}(\eta ,Z)=-\underset{i=1}{\overset{N}{\Sigma }}\log \underset{k:{\alpha }_k>0}{\Sigma }{\alpha }_{k}p\left(t_i\right|{\eta }_k)p\left(x_i\right|{\eta }_k)+\frac{c}{4}\underset{k:{\alpha }_k>0}{\Sigma }\log \left(\frac{N{\alpha }_k}{12}\right)+\frac{K_{\mathrm{nz}}}{2}\log \frac{N}{12}+\frac{K_{\mathrm{nz}}(c+1)}{2},---\left(36\right)$>

所以,基于模糊聚类的MML准则为

$>\hat{\eta }=\arg \underset{\eta }{\min }\{-\underset{i=1}{\overset{N}{\Sigma }}\log \underset{k:{\alpha }_k>0}{\Sigma }{\alpha }_{k}p\left(t_i\right|{\eta }_k)p\left(x_i\right|{\eta }_k)+\frac{c}{4}\underset{k:{\alpha }_k>0}{\Sigma }\log \left(\frac{N{\alpha }_k}{12}\right)+\frac{K_{\mathrm{nz}}}{2}\log \frac{N}{12}+\frac{K_{\mathrm{nz}}(c+1)}{2}\}.---\left(37\right)$>

图6为本发明仿真数据的最优模糊分割结果示意图，如图6所示，本实施例方法将所述最优分割阶数与最优聚类参数代入模糊聚类算法得出最优模糊分割结果。

具体来说，处理器采用高斯函数得到在时间坐标集合上k_opt个模糊分割β_k(t_i)，公式为

$>{\beta }_k\left(t_i\right)=\frac{A_k\left(t_i\right)}{{\Sigma }_{j=1}^{k_{\max }}{A}_j\left(t_i\right)},1\le k\le k_{\mathrm{opt}}---\left(38\right)$>

其中t_i为气候变量数据对应的时间，μ_k为在时间坐标集合上的初始模糊聚类中心，为时间坐标集合上的初始偏差。

图7为本发明仿真数据的最优隶属度结果示意图，如图7所示，本实施例方法将所述最优分割阶数与最优聚类参数代入模糊聚类算法得出最优隶属度结果。

具体来说，采用公式

$>u_{i,k}=\frac{{[{\alpha }_{k}p\left(t_i\right|{\eta }_k)p\left(x_i\right|{\eta }_k)]}^{m-1}}{{\Sigma }_{j=1}^{k_{\max }}{\left[{\alpha }_{j}p\left(t_i\right|{\eta }_j)p\left(x_i\right|{\eta }_j)\right]}^{m-1}}---\left(39\right)$>

计算隶属度，其中，所述为第k个聚类的混合概率，p(t|η_k)为在时间坐标集合上的第k个高斯函数，p(x|η_k)为在气候变量数据集合上第k个高斯函数，m为模糊加权指数且m＝2；

图8为本发明仿真数据的最优明确分割结果示意图，如图8所示，本实施例方法将所述最优分割阶数与最优聚类参数代入模糊聚类算法得出最优明确分割结果。

具体来说，采用公式

$>k^{*}=\arg \underset{1\le k\le k_{\mathrm{opt}}}{\max }{u}_{i,k},1\le i\le N---\left(40\right)$>

计算最优明确分割结果。

本实施例中，时间序列模糊分割算法的具体为：

(1)输入数据集合设置终止阈值ε，初始分割阶数k_max，最小分割阶数k_min。

(2)初始聚类参数其中其中μ_k与分别表示在时间坐标集合上的第k个高斯函数p(t|η_k)的中心和方差，v_k与Σ_k分别表示在数据集合上第k个高斯函数p(x|η_k)的中心和方差，

$>p\left(t_i\right|{\eta }_k)=\frac{1}{\sqrt{2\pi {\sigma }_k^2}}\exp (-\frac{1}{2}\frac{{(t_i-{\mu }_k)}^2}{{\sigma }_k^2})---\left(41\right)$>

$>p\left(x_i\right|{\eta }_k)=\frac{1}{{\left(2\pi \right)}^{\frac{d}{2}}\sqrt{\det \left({\Sigma }_k\right)}}\exp (-\frac{1}{2}{(x_i-v_k)}^T{\Sigma }_k^{-1}(x_i-v_k)).---\left(42\right)$>

(3)计算隶属度

$>u_{i,k}^{\left(l\right)}=\frac{{[{\alpha }_k^{\left(l\right)}p\left(t_i\right|{\eta }_k^{\left(l\right)})p\left(x_i\right|{\eta }_k^{\left(l\right)})]}^{m-1}}{{\Sigma }_{j=1}^{k_{\max }}{\left[{\alpha }_j^{\left(l\right)}p\left(t_i\right|{\eta }_j^{\left(l\right)})p\left(x_i\right|{\eta }_j^{\left(l\right)})\right]}^{m-1}}---\left(43\right)$>

(4)更新混合概率

$>{\alpha }_k^{(l+1)}=\frac{\max \{0,\left(\underset{i=1}{\overset{N}{\Sigma }}{u}_{i,k}^{\left(l\right)}\right)-\frac{c}{2}\}}{\underset{j=1}{\overset{k_{\max }}{\Sigma }}\max \{0,\left(\underset{i=1}{\overset{N}{\Sigma }}{u}_{i,k}^{\left(l\right)}\right)-\frac{c}{2}\}}---\left(44\right)$>

$>\{{\alpha }_1^{(l+1)},...,{\alpha }_{k_{\max }}^{(l+1)}\}=\frac{\{{\alpha }_1^{(l+1)},...,{\alpha }_{k_{\max }}^{(l+1)}\}}{{\Sigma }_{k=1}^{k_{\max }}{\alpha }_k^{(l+1)}}---\left(45\right)$>

(5)如果$>{\alpha }_k^{(l+1)}>0,$>更新聚类参数$>{\eta }_k^{(l+1)}=\{{\mu }_k^{(l+1)},{\left({\sigma }_k^2\right)}^{(l+1)},v_k^{(l+1)},{\Sigma }_k^{(l+1)}\}$>和隶属度

$>{\mu }_k^{(l+1)}=\frac{{\Sigma }_{i=1}^N{\left(u_{i,k}^{\left(l\right)}\right)}^m{t}_i}{{\Sigma }_{i=1}^N{\left(u_{i,k}^{\left(l\right)}\right)}^m}---\left(46\right)$>

$>{\left({\sigma }_k^2\right)}^{(l+1)}=\frac{{\Sigma }_{i=1}^N{\left(u_{i,k}^{\left(l\right)}\right)}^m(t_i-{\mu }_k^{\left(l\right)}){(t_i-{\mu }_k^{\left(l\right)})}^T}{{\Sigma }_{i=1}^N{\left(u_{i,k}^{\left(l\right)}\right)}^m}---\left(47\right)$>

$>v_k^{(l+1)}=\frac{{\Sigma }_{i=1}^N{\left(u_{i,k}^{\left(l\right)}\right)}^m{x}_i}{{\Sigma }_{i=1}^N{\left(u_{i,k}^{\left(l\right)}\right)}^m}---\left(48\right)$>

$>{\Sigma }_k^{(l+1)}=\frac{{\Sigma }_{i=1}^N{\left(u_{i,k}^{\left(l\right)}\right)}^m(x_i-v_k^{\left(l\right)}){(x_i-v_k^{\left(l\right)})}^T}{{\Sigma }_{i=1}^N{\left(u_{i,k}^{\left(l\right)}\right)}^m}---\left(49\right)$>

$>u_{i,k}^{(l+1)}=\frac{{\left[{\alpha }_k^{(l+1)}p\left(t_i\right|{\eta }_k^{(l+1)})p\left(x_i\right|{\eta }_k^{(l+1)})\right]}^{m-1}}{{\Sigma }_{j=1}^{k_{\max }}{\left[{\alpha }_j^{(l+1)}p\left(t_i\right|{\eta }_j^{(l+1)})p\left(x_i\right|{\eta }_j^{(l+1)})\right]}^{m-1}}---\left(50\right)$>

(6)如果降低当前分割阶数k_nz，令k_nz＝k_nz-1。如果k_nz<k_min，算法结束；令k＝k+1，如果k<k_max，回到步骤(3)。

(7)计算最小信息长度准则

$>\left(\begin{array}{c}\mathrm{MessLen}({\eta }^{(l+1)},Z)=-\underset{i=1}{\overset{N}{\Sigma }}\log \underset{k:{\alpha }_k>0}{\Sigma }{\alpha }_{k}p\left(t_i\right|{\eta }_k^{(l+1)})p\left(x_i\right|{\eta }_k^{(l+1)})\\ +\frac{c}{4}\underset{k:{\alpha }_k>0}{\Sigma }\log \left(\frac{N{\alpha }_k}{12}\right)+\frac{k_{\mathrm{nz}}}{2}\log \frac{N}{12}+\frac{k_{\mathrm{nz}}(c+1)}{2}.\end{array}\right)---\left(51\right)$>

(8)如果||U^(l+1)-U^(l)||≥ε，令l＝l+1，回到第(3)步。

(9)如果MessLen(η^(l+1),Z)≤MessLen_min，令

MessLen_min＝MessLen(η_(l+1),Z)   (52)

η_opt＝η_(l+1)   (53)

(10)令k^*＝argmin_k{α_k>0},k_nz＝k_nz-1。如果k_nz≥k_min，回到步骤(3)。

本发明实施例处理器接收气候变量数据以及与所述气候变量数据对应的时间值，采用模糊聚类算法计算所述气候变量数据的渐进变化与跳跃变化，并根据所述渐进变化与所述跳跃变化对齐气候变化进行预测以及趋势分析，准确地体现了气候变化中的跳跃变化，有利于进一步准确预测资源的合理利用及区域可持续发展，有助于分析对应该跳跃变化的成因。

最后应说明的是：以上各实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述各实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的范围。