利用云气候学数据预判目标区域未来月平均云量的方法

摘要

本发明提供了一种利用云气候学数据预判目标区域未来月平均云量的方法。该方法基于目标区域历史云量的平稳性和季节周期性特征，选用不同的预判模型来预判目标区域未来月平均云量，具有可靠性高、应用广泛的优势，可广泛应用于经济、气象、农业、生物科学等领域。

说明书

技术领域

本发明涉及卫星遥感技术领域，尤其涉及一种利用云气候学数据预判 目标区域未来月平均云量的方法。

背景技术

卫星遥感技术在生态环境保护、防灾减灾、资源调查、城市规划、地 图测绘等领域发挥了重要作用。随着遥感技术的发展，遥感图像空间分辨 率的不断提高，导致遥感数据特别是高空间分辨率数据量急剧增加，这对 于卫星资源的高效利用提出了更高的要求。由于天空中云层覆盖的影响， 光学传感器在获取遥感数据时经常无法获取满足用户需求的无云或少云 遥感数据。而多云遥感数据直接影响着遥感数据在各领域的应用效果。因 而如何预判未来一年月平均云量，为用户获取优质卫星遥感数据提供参考 就成为了一个重要问题。

国际卫星云气候研究计划(ISCCP)成立于1982年，作为世界气候研 究计划(WCRP)之一，它自1983年7月起通过气象卫星收集全球云空间分 布和时间分布信息。为了研究方便，ISCCP采用等面积投影方式，对全球 范围按照标准网格划分。例如D2数据集单个网格大小为2.5°×2.5°(约280 公里*280公里)，全球共划分为6596个网格。根据ISCCP多年来收集的 数据，国内外学者对云的空间和时间分布及其变化规律进行了大量的研究， 但基本不涉及未来云量的预判问题。

现有技术一般采用月平均云量年度平均的方法计算未来一年月平均 云量的方法。如图1所示，以ISCCP的D2数据集为例，该方法包括：

步骤S101，从D2数据集中获取某一网格的n年某月的月平均云量序 列C_i，j。

其中，i为网格编号，取值范围[1，6596]；j为月份，取值范围[1，12]， C_i，j为12×n个数组成的序列；

步骤S102，将n年的月平均云量平均值作为作为未来一年j月份 的平均云量预判值。

然而，在实现本发明的过程中，申请人发现上述的预判未来月云量的 方法存在如下缺陷：上述方法没有利用考虑到不同区域月平均云量序列的 特征，简单的取平均预判未来月云量导致预判结果准确度较低。

发明内容

(一)要解决的技术问题

鉴于上述技术问题，本发明提供了一种利用云气候学数据预判目标区 域未来月云量的方法，以提高预判结果的准确性。

(二)技术方案

本发明提供了一种利用云气候学数据预判目标区域未来月平均云量 的方法，该方法包括：步骤A，读取ISCCP的D2数据集中目标区域所在 网格的预设年份段的月平均云量序列C；步骤B，采用迪基-福勒检验法检 验所述月平均云量序列C的平稳性，如果所述月平均云量序列C的数据平 稳，执行子步骤D1；以及子步骤D1，利用自回归滑动平均模型对月平均 云量进行处理，得到未来月平均云量。

(三)有益效果

从上述技术方案可以看出，本发明利用云气候学数据预判目标区域未 来月云量的方法具有以下有益效果：

(1)根据目标区域的历史云量特征选择的不同的预判模型，能够适 应不同区域的云量变化特点，提高预判精度；

(2)自回归滑动平均(ARMA)模型、求和自回归滑动平均(ARIMA) 模型和季节性求和自回归滑动平均(SARIMA)模型能够充分利用时间序 列自身的规律性，得出较为准确的预判值，被广泛应用于经济、气象、农 业、生物科学等领域，本发明结合历史月平均云量序列特征，选择上述三 种不同的模型之一对目标区域的未来月平均云量进行预判，可以供遥感卫 星订购用户做参考，预判目标区域有效数据获取的可行性。

附图说明

图1为本发明实施例利用云气候学数据预判目标区域未来月云量的方 法的流程图；

图2为对于北京地区利用1984年到1996年12年数据预判1997年12 个月份的预判图；

图3为对于乌鲁木齐地区利用1984年到2008年24年数据预判2009 年12个月份的预判图。

具体实施方式

在本发明的一个示例性实施例中，提供了一种利用云气候学数据预判 目标区域未来月云量的方法。请参照图1，该方法包括四个步骤：

步骤A：读取ISCCP的D2数据集中目标区域所在网格的预设年份段 【1984-2012】的月平均云量序列C：

$C=\left(\begin{array}{c}{C}_{1984}\\ C_{1995}\\ ···\\ C_i\\ ···\\ C_{2012}\end{array}\right)---\left(1\right)$

其中，C为18×12行1列的向量，$C_i=\left(\begin{array}{c}{C}_{i1}\\ C_{i2}\\ ···\\ C_{\mathrm{ij}}\\ ···\\ C_{i12}\end{array}\right),$i为年份，j代表月份。

读取目标区域所在网格的云气候学数据中的月平均云量数据可以按 照ISCCP数据集的说明[W.B.Rossow，A.W.Walker，D.E.Beuschel and  M.D.Roiter，“International satellite cloud climatology project(ISCCP) documentation of new cloud datasets.”WMO/TD737，World Meteorological  Organization，1996，pp.115.]进行。

此外，为了满足预判的准确性，预设年份段应当大于等于10年的月 平均云量数据，即120个历史月平均云量观测值。

步骤B：采用迪基-福勒检验法(Dickey-Fuller test，DF)检验月平均 云量序列C的平稳性，如果该月平均云量序列C的数据平稳，执行步骤 D1；如果月平均云量序列C的数据不平稳，执行步骤C；

该检验月平均云量序列C的平稳性的步骤B进一步包括：

子步骤B1，建立检验的原假设H₀：序列为非平稳序列；相应的备择 假设H₁：序列为平稳序列；

子步骤B2，在原假设H₀的基础上，以月平均云量序列C为样本数据， 计算其检验统计量t_τ；

子步骤B3，查找DF统计量百分位表中与给定的显著性水平α相应的 分位点；

子步骤B4，将检验统计量t_τ与该分位点进行比较，如果小于该分位点， 则拒绝原假设，序列为平稳序列，如果大于该分位点，则接受该原假设， 该序列为非平稳序列。

步骤C，由自相关函数ACF(Autocorrelation function)判断月平均云 量序列C的季节周期性特征，如果不满足季节周期性特征，执行步骤D2； 如果满足季节周期性特征，执行步骤D3；

该判断月平均云量序列C的季节周期性特征的步骤C进一步包括：

子步骤C1，计算月平均云量序列C的第12阶自相关系数r₁₂；

$r_{12}=\frac{e_{12}}{e_0}=\frac{\frac{1}{N}{\Sigma }_{t=1}^{N-12}(c_t-\bar{c})(c_{t+12}-\bar{c})}{{\Sigma }_{t=1}^N{(c_t-\bar{c})}^2}---\left(2\right)$

其中，为月平均云量序列C的平均值，N为月平均云量序列C的元 素个数，c_t是月平均云量序列C第t个元素值，c_t+k是月平均云量序列C 第t+k个元素值，

子步骤C2，判断该第12阶自相关系数r₁₂与阈值的关系：如果该 第12阶自相关系数r₁₂大于阈值则说明该平均云量时间序列满足季节 周期性特征，执行步骤D2，否则说明该平均云量时间序列不满足季节周 期性特征，执行步骤D3；

步骤D，针对不同类型的月平均云量序列C建立相应的模型，进行月 平均云量的预判，该步骤D包括：

子步骤D1，对于满足平稳性的月平均云量序列C，利用自回归滑动 平均模型对月平均云量进行处理，得到未来月平均云量预判值；

自回归滑动平均(auto regression moving average，ARMA)模型，简 记为ARMA(p，q)，

c_t＝φ₁c_t-1+φ₂c_t-2+…+φ_pc_t-p+ε_t-θ₁ε_t-1-θ₂ε_t-2-…-θ_qε_t-q   (3)

其中，c_t为月平均云量序列，φ₁，φ₂，...，φ_p为自回归系数，ε_t，ε_t-1，...ε_t-q为未观测到的白噪声，θ₁，θ₂，...，θ_q为滑动平均系数。p为自回归系数的个 数；q为滑动平均系数的个数。

为了方便表示，引进自回归算子Φ(B)＝1-φ₁B-φ₂B²-…-φ_pB^p和滑动 平均算子θ(B)＝e_t-θ₁B-θ₂B²-…-θ_qB^q自回归滑动平均模型ARMA(p，q) 简写为

Φ(B)c_t＝Θ(B)ε_t         (4)

基于上述原理分析，该子步骤D1包括：

子分步骤D1a，使用以下参数初始值：......，

θ₁＝-0.90，θ₂＝-0.85，...，θ_q＝-0.95+0.05q；ε_t＝0

按照p＝1，2，...，12，q＝1，2，...12，共144组可能组合，分别建立自回归 滑动平均模型，每一自回归滑动平均模型的参数包含p阶自回归系数和q 阶滑动平均系数，即......、θ₁、θ₂、......、θ_q；

把月平均云量序列C的序列值及参数估计初始值代入下述方程：

e_t＝c_t-φ₁c_t-1-φ₂c_t-2-…-φ_pc_t-p-ε_t+θ₁ε_t-1+θ₂ε_t-2+…+θ_qε_t-q    (5)

子分步骤D1b，利用方程(5)针对144组组合分别建立相应的最小 化目标函数：

$S_c({\phi }_1,{\phi }_2,...,{\phi }_p,{\theta }_1,{\theta }_2,...,{\theta }_q)=\underset{t=2}{\overset{N}{\Sigma }}{e}_t^2---\left(6\right)$

利用高斯牛顿法对144个目标函数的每个函数 S_c(φ₁，φ₂，...，φ_p，θ₁，θ₂，...，θ_q)中的参数，即 ......、θ₂、......、θ_q进行最小二乘估计，得到144组 p阶自回归系数和q阶滑动平均系数参数值，即 ......、θ₁、θ₂、......、θ_q；

子分步骤D1c，利用残差分析方法，从144组参数值中选出最优的一 组参数组合作为最终的模型参数，即p阶自回归系数和q阶滑动平均系数， 确定最终的自回归滑动平均模型；

小步骤D1c1，对于每组组合的参数值，利用以下公式计算其残差；

e_t＝c_t-φ₁c_t-1-φ₂c_t-2-…-φ_pc_t-p-ε_t+θ₁ε_t-1+θ₂ε_t-2+…+θ_qε_t-q   (7)

其中e_t为残差，c_t为月平均云量序列C中第t个元素的值，其余项为 拟合值。

小步骤D1c2，计算每一组参数值对应残差的标准差；

小步骤D1c3选用赤池信息准则(Akaike Information criterion，AIC) 检验标准，即

其中，σ_e是残差序列e_t的标准差，n为月平均云量序列C的数量，p， q为上述144种组合中的p，q的取值，得到144个AIC数值；

小步骤D1c4，选取AIC最小的一组模型参数组合为最终模型参数值， 即p阶自回归系数和q阶滑动平均系数，确定最终的自回归滑动平均模型；

子步骤D1d中，将月平均云量序列C及确定的最终自回归滑动平均 模型参数代入以下预测方程：

$\hat{c}\left(r\right)={\Phi }_1\hat{c}(r-1)+{\Phi }_2\hat{c}(r-2)+···+{\Phi }_p\hat{c}(r-p)+{\Theta }_0-$

${\Theta }_{1}E\left(e_{N+r-1}\right|c_1,c_2,...,c_N)-{\Theta }_{2}E\left(e_{N+r-2}\right|c_1,c_2,...,c_N)-···-$

${\Theta }_{q}E\left(e_{N+r-q}\right|c_1,c_2,...,c_N)---\left(9\right)$

其中，$E\left(e_{N+r-1}\right|c_1,c_2,...,c_N)=\left(\begin{array}{cc}0& r>0\\ e_{N+r}& r\le 0\end{array}\right),$c₁，c₂，...，c_N为月平均云量 序列C的元素，，φ₁，φ₂，...，φ_p和θ₁，θ₂，...，θ_q为p阶自回归系数和q阶滑动 平均系数。e_N+r-1为残差。

计算得到未来第r个月的月平均云量预判值；

令r＝1，2，...，12，就可以分别得到未来12个月的月平均云量预判值。

子分步骤D1e，由以下公式计算相应月份给定置信区间的月平均云量 区间。

$(\hat{c}\left(r\right)-\sigma u_{\frac{\alpha }{2}},\hat{c}\left(r\right)+\sigma u_{\frac{\alpha }{2}})---\left(10\right)$

其中，为月平均云量序列预判值，σ为预测误差序列的标准差， 即r＝1，2，..，12.，α＝1-b，b为给定置信区间， 可根据α通过查正态分布表得到。

本实施例中，分别令α＝0.2，0.05，可得到置信区间为80％和95％ 的月平均云量区间。

子步骤D2，对于不满足平稳性，同时不满足季节性周期性的月平均 云量序列数据，建立求和自回归滑动平均(autoregressive integrated moving  average，ARIMA)模型对月平均云量序列数据进行处理，得到未来月平均 云量预判值；

求和自回归滑动平均(autoregressive integrated moving average)模型， 简记为ARIMA(p，d，q)模型：

$\Phi \left(B\right){▿}^d{c}_t=\Theta \left(B\right){ϵ}_t---\left(11\right)$

式中，表示对月平均云量序列C进行d次差分。差分后建立与步 骤C1相同的模型，求解自回归算子和滑 动平均算子θ(B)＝e_t-θ₁B-θ₂B²-…-θ_qB^q中的自回归系数和滑动平均系 数。方法与步骤D1a～D1d相同。

子分步骤D2a，对月平均云量序列C进行一次差分处理，生成一次差 分月平均云量序列：

$C_d=\begin{array}{c}\left(\begin{array}{c}{c}_2-c_1\\ c_3-c_2\\ c_4-c_3\\ ·······\\ c_{N-1}-c_{N-2}\\ c_N-c_{N-1}\end{array}\right)---\left(12\right)\end{array}$

子分步骤D2b，利用一次差分月平均云量序列C_d，进行步骤D1a～D1c 的操作，生成一次差分月平均云量序列C_d的月平均云量预判值，即 c_d(N+1)，c_d(N+2)，…，c_d(N+r)；

子分步骤D2c，进行差分还原操作得到最终月平均云量预判值，利用 以下公式计算未来12个月的月平均云量预判值：

$\left(\begin{array}{c}{c}_{N+1}=c_{d\left(N\right)},+c_N;\\ c_{N+2}=c_{d(N+1)}+c_{N+1;}\\ ...;\\ c_{N+12}=c_{d(N+11)}+c_{N+11};\end{array}\right)---\left(13\right)$

子分步骤D2d，得到最终月平均云量预判值后，利用分步骤D1e的操 作计算生成80％置信区间和95％置信区间。

子步骤D3，对于满足季节周期性特征的月平均云量序列C，建立季 节性求和自回归滑动平均(Seasonal autoregressive integrated moving  average，SARIMA)模型对月平均云量进行处理，得到未来月平均云量预 判值；

季节性求和自回归滑动平均(Seasonal autoregressive integrated moving  average，SARIMA)模型可以表示为：

${▿}^d{{▿}_S}^D{c}_t=\frac{\Theta \left(B\right){\Theta }_S\left(B\right)}{\Phi \left(B\right){\Phi }_S\left(B\right)}{ϵ}_t---\left(14\right)$

式中，表示对数据进行季节性差分，S表示季节周期，本方法中， 季节差分阶数D选择1。由于为每月一个数据点，因此季节周期S为12。 表示对月平均云量序列C进行d次差分。

季节性差分处理后建立与步骤D2相同的模型，求解差自回归算子 和滑动平均算子 θ(B)＝e_t-θ₁B-θ₂B²-…-θ_qB^q中的自回归系数和滑动平均系数的方法与 步骤D1相同。

子分步骤D3a，对月平均云量序列C进行季节性差分操作，生成季节 差分月平均云量序列C_S：

$c_s=\begin{array}{c}\left(\begin{array}{c}{c}_{13}-c_1\\ c_{14}-c_2\\ c_{15}-c_3\\ ·······\\ c_{N-1}-c_{N-13}\\ c_N-c_{N-12}\end{array}\right)---\left(14\right)\end{array}$

子分步骤D3b，利用季节差分月平均云量序列C_S，进行步骤D1a～D1c 的操作，生成季节差分月平均云量序列C_S的月平均云量预判值，即 c_d(N+1)；c_d(N+2)；......；c_s(N+r)；

子分步骤D3c，进行季节差分还原操作得到最终月平均云量序列预判 值，利用以下公式计算未来12个月的月平均云量序列预判值：

$\left(\begin{array}{c}{c}_{N+1}=c_{s(N-11)},+c_{N-11};\\ c_{N+2}=c_{s(N-\text{10})}+c_{N-10;}\\ ...;\\ c_{N+12}=c_{s\left(N\right)}+c_N;\end{array}\right)---\left(15\right)$

子分步骤D3d，利用D1e步骤中操作计算得到月平均云量序列预判值 的80％置信区间和95％置信区间。

至此，已经结合附图对本实施例进行了详细描述。依据以上描述，本 领域技术人员应当对本发明利用云气候学数据预判目标区域未来月平均 云量的方法有了清楚的认识。以下给出应用本实施例方法进行北京和乌鲁 木齐两地云量预判的具体结果。

实际应用实例1：

以北京所在网格为例，在时间上，选取2007年份进行预判和验证。用 1984年到2006年22年数据预判2007年12个月云量。得出2007年12 个月份的预判值，预判值80％、95％置信区间。与2007年实际值比较， 实际真实值均落在80％置信区间内，平均相对误差3.60％，最大相对误差 5.29％，均方误差3.2199。

表格1北京地区2007年预判值和实际值

实际应用实例2：

以乌鲁木齐所在网格为例，在时间上，选取2009年份进行预判和验 证。用1984年到2008年24年数据预判2009年12个月云量，得出2009 年12个月份的预判值，预判值80％、95％置信区间。与2009年实际值比 较，实际真实值除了一月份(落在95％置信区间)外，均落在80％置信区 间内，平均相对误差2.80％，最大相对误差6.55％，均方误差2.8871。

表格2乌鲁木齐地区2009年预判值与实际值

此外，上述方法并不仅限于实施方式中提到的各种具体算法，本领域 的普通技术人员可对其进行简单地熟知地替换，例如：

(1)步骤C1中参数估计环节提到的高斯牛顿法还可以使用最速梯度 法或阻尼最小二乘法代替；

(2)步骤C1中提到的赤池信息准则(Akaike Information criterion， AIC)还可以使用贝叶斯信息准则(Bayesian Information Criterion，BIC)。

综上所述，本发明提供一种利用云气候学数据预判目标区域未来月平 均云量的方法。该方法基于目标区域历史云量的平稳性和季节周期性特征， 选用不同的预判模型来预判目标区域未来月平均云量，具有可靠性高、应 用广泛的优势，可广泛应用于农业、资源调查、城市规划、地图测绘等领 域光学卫星数据获取的可行性分析。

以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行 了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而 已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修 改、等同替换、改进等，均应包含在本发明的保护范围之内。