基于小波随机耦合模型在流域流量方面应用.docx
《基于小波随机耦合模型在流域流量方面应用.docx》由会员分享，可在线阅读，更多相关《基于小波随机耦合模型在流域流量方面应用.docx（26页珍藏版）》请在冰豆网上搜索。
基于小波随机耦合模型在流域流量方面应用
基于小波随机耦合模型在流域流量方面的应用
摘要：
小波分析是一种时、频多分辨率分析方法，能对时间序列进行多尺度分析，它不仅可以充分展示水文时间序列的精细结构，而且还便于提取水文时间序列的变化规律。
流域流量既是这样一个复杂的时间序列，其变化与区域气候条件和生态环境密切相关，具有明显的年际同期性和随机性波动，用小波随机耦合模型描述流域流量的动态变化具有重要意义和深远的影响。
关键词：
小波分析;随机模型;耦合预测；地下水埋深
Basedonwaveletrandomcouplingmodelintheriverflowapplications
Abstract:
Waveletanalysisisakindofwhenmultiresolutionanalysismethod,thefrequency,thetimeseriesofmultiscaleanalysis,itnotonlycanfullydisplayinhydrologicaltimeseriesoffinestructure,butalsoeasytoextractthechangeruleofhydrologicaltimeseries.Riverflowissuchacomplextimeseries,thechangeandregionalclimateconditionsandecologicalenvironmentrelatedclosely,haveobviousinterannualTongQiXingandrandomfluctuation,usingwaveletrandomcouplingmodeltodescribethedynamicchangeofflowriverbasinhasimportantsignificanceandfar-reachinginfluence.
Keywords:
Waveletanalysis；Randommodel；Couplingforecast；Groundwaterdepth
1前言
水文随机模拟已经取得了不断的进步和广泛应用。
小波分析在水文水资源学中得到了广泛的研究和应用[1]。
文献[2-4]将小波分析与目前成熟的系统分析方法耦合建立了水文预测模型，结果表明其途径是可行的。
流域流量预测有多种方法，文献中[5-9]分别采用时间序列分析、人工神经网络、灰色理论等方法对流域流量进行了预测，但这些方法都属于单一预测方法，难以获得理想的预测效果。
由于流域流量变化过程属于非平稳随机过程，流量的的大小受到多种复杂因素的影响，因此，要对流量进行准确预测是非常困难的。
近些年发展起来的小波分析方法可以充分展示水文时间序列的精细结构，便于提取水文时间序列的变化规律[10]。
2小波分析的发展概述
2.1小波分析的历史回顾
小波分析[11]是自1986年以来由于Y．Meyer、S．Mallat及I．Daubechies等人的奠基工而迅速发展起来的新兴学科，它是Fourier分析划时代发展的结果。
然而，它的展历史可以追溯到1909年Haar的工作。
从现代小波分析的观点来看，1930年前后许多与小波有关的新方向出现，但以后进展不大，直到1960年Calder6n的研究及20后（1980年）Grossmann与Morlet的研究，后人称为“原子分解”，才有了现代小波分的雏形。
进而，由于1986年以后的工作以及应用的广泛性使这个学科得到了飞速发展。
2.2现代小波分析的发展阶段[11]。
第一阶段为孤立研究与应用阶段。
主要特征是得到了一些特殊的小波，并在一些专业领域得到零散应用。
一个代表性的工作是法国地质学家A．Grossmann和J．Morlet于1984年第一次把“小波”用于处理地质数据，得到了以他们的名字命名的Grossmann-Morlet小波。
另一个代表性的工作是1981年J．Stromberg与他的合作者发现的正交小波基。
这个时期一个有趣的现象是各个领域的专家、学者和工程师在完全不了解别人的研究工作的状况下巧妙地、独立地构造自己的“小波”。
第二阶段为小波分析理论成型期。
这个时期，法国数学家Y．Meyer成功地构造出了具有一定衰减性质的光滑函数
，这个函数的二进伸缩和整数平移产生的函数组
构成空间L2（R）的规范正交基。
这引起了大家广泛关注，因为在此之前，学术界普遍认为不可能具有如此好性质的函数，Meyer也是认为这样好的函数不可能存在，在构造反例的过程中，把这个正交小波构造出来。
随后，P．Lemarie和G．Battle又分别独立地构造出同样“好”的小波。
第三阶段为小波理论和应用全面发展期。
一般认为，它是从20世纪90年代开始的。
它的主要特征是，在上述理论框架下，出现了许多很有价值的应用成果，也解决了长期没有解决的应用问题，相反地，在应用中也提出了许多需要解决的问题，从而推动了小波分析理论的发展。
3小波随机耦合模型建模思路
首先将研究的水文时间序列采用快速小波变换算法（例如ATrous算法）进行小波分解，得到某尺度P下的小波变换序列{W1（t），W2（t），…，Wp（t），Cp（t）}；然后对各小波变换序列的主要成分（随机成分或确定成分）进行识别，对各小波变换序列进行互相关分析，并建立各小波变换序列适宜的数学模型；最后采用小波变换重构算法得到所研究水文时间序列的小波随机耦合模型。
4小波随机耦合模型在诺敏河流量变化预测中的应用
根据诺敏河流域99-07年的实际观测资料，对该流域流量变化进行模拟，并预测未来发展趋势，为有关决策者提供参考依据。
4.1实测月流量序列平稳化处理
由图1可以看出，诺敏河流量序列Ht为一非平稳时间序列，要使Ht转换为平稳时间序列，可以对其进行差分处理．由于Ht在均值水平上不平稳，因此认为只需对其进行一次差分即可达到平稳，即：
式中:
B为后移算子通过差分，使得诺敏河流量序列Ht由非平稳序列转化为平稳序列Xt，见图2。
图1诺敏河原始流量变化曲线
图2差分后诺敏河原始流量变化曲线
4.2差分后逐月诺敏河流量序列小波分解与重构
采用ATrous算法，取尺度数P=2，对差分后年降水标准化序列Xt进行分解，得到小波分解序列{W1（t），W2（t），…，Wp（t），Cp（t）}；将各小波分解序列进行叠加，得到重构序列。
重构过程与差分后序列Xt变化过程完全一致，因此，采用ATrous算法对诺敏河流量序列进行分解是可行的。
图3小波分解与重构
4.3小波变换序列成分识别
采用公式
和公式
计算各小波变换序列的自相关系数和方差谱密度，并绘制自相关图和方差谱密度图，见图和图。
根据结果图分析各小波变换序列的变化特性，可以近似认为W1（t）代表序列的随机项，W2（t）代表序列的周期项，C2（t）代表序列的趋势项。
图4各小波变换序列自相关图
图5小波变换序列方差谱密度图
4.4小波变换序列互相关分析
采用下列公式分别计算各小波变换序列的互相关系数，绘制互相关图，并加绘95％容许限。
由图可以看出，各小波变换序列互相关系数基本上落在95%容许限范围以内，且趋近于0.因此，各小波变换序列互相关性较小，可以认为小波变换序列W1（t）、W2（t）和C2（t）两两独立。
上述互相关分析结果表明，小波变换序列W1（t）、W2（t）和C2（t）成分单一，比序列Xt要简单，因此分析和处理Xt就转嫁为对W1（t）、W2（t）和C2（t）进行处理。
图6各小波变换序列互相关图
4.5建立细节序列W1（t）的随机模型
通过计算机编程计算，小波分解细节序列W1（t）的均值W1t=-0.0394≈0，方差σ2W1=0.1.8573，偏态系数Csw1=0.8684，趋近于0，所以认为序列W1（t）近似于正态分布，不必进行正态性转化。
对序列W1（t）分别进行自相关分析和偏相关分析，自相关图具有拖尾性，而偏相关图具有截尾性，所以初步判定模型形式为AR（p）模型[12-13]。
判定模型阶数为1，属于AR
（1）模型。
对AR
（1）模型参数进行计算，建立如下自回归模型:
W1m（t）=0.5557*W1（t-1）+0.3563*W1（t-2）-0.3035*W1（t-12）
采用BIC准则对AR（p）模型的阶数进行进一步识别。
当p=1时，BIC达到最小值。
这说明初步确定的模型阶数为1阶是合适的。
采用自相关系数综合检验法检验残差项εt是否为独立序列。
经过计算，统计Q=5.3570，n服从（21-12）的卡方分布，卡方值为16.919，所以残差是相互独立的假定是可以接受的。
对独立随机序列εt的正态性进行检验:
均值ε-t=-0.03943≈0，方差σ2ε=0.9252，偏态系数Cs（ε）=0.0081，趋近于0，因此近似认为服从正态分布。
因此，细节序列W1（t）随机项模型为:
W1m（t）=0.5557*W1（t-1）+0.3563*W1（t-2）-0.3035*W1（t-12）
4.6建立细节序列W2（t）和背景序列C2（t）的自回归模型
细节序列W2（t）和背景序列C2（t）为确定成分，可以借助于水文学中的自回归模型（略去随机变量εt）对序列W2（t）和C2（t）进行描述。
分别对序列W2（t）和C2（t）进行自相关分析和偏相关分析，并采用BIC准测，判定所需建立的W2（t）和C2（t）自回归模型阶数分别为4阶和3阶。
对于序列W2（t），选定k=1，2，3，4（偏相关系数超出95%容许限范围），则略去随机变量εt的AR（4）模型为:
W2m（t）=0.6240*W2（t-1）+1.5676*W2（t-2）-0.8358*W2（t-3）-1.1532*W2（t-4）；
对于序列C2（t），选定k=1，2，3（偏相关系数超出95%容许限范围），则略去随机变量εt的AR（3）模型为:
C2m（t）=0.7618*C2（t-1）+0.8558*C2（t-2）-0.6245*C2（t-3）；
4.7小波随机耦合模型的组合与拟合
Xm（t）=W1m（t）+W2m（t）+C2m（t）;
Hm（t）=Xm（t）+H（t-12）;
诺敏河流量拟合模型如图7所示
图7诺敏河流量小波随机耦合模型拟合曲线
4.8小波随机耦合模型精度检验
采用所建小波随机耦合模型的拟合数据进行拟合效果检验，得出后验差比值C=0.0721，小误差频率P=1，相对均方差E=0.0034，拟合准确率E2=1，达到一级标准。
结论
1）利用小波随机耦合模型来研究诺敏河流量的动态变化，能合理的模拟真实的流域情况，小波随机耦合模型原理简单，计算快捷，具有一定的先进性和使用性。
2）对于实际数据的拟合精度还存在一定的误差，从而发现理论的研究方法同样有其适用的数据特点。
参考文献
[1]张凯.关于天津市水资源价值与水价格的研究[D].天津:
天津大学管理学院.2006.
[2]姚雨田.天津市水资源优化配置与节水研究[J].天津建设科技，1996，（4）.
[3]韩海燕，于苏俊.基于多目标优化模型的荣县水资源优化配置[J].安徽农业科学，2010，38（19）:
200-202.
[4]张成凤.基于遗传算法的榆林市水资源优化配置的研究[D].陕西:
西北农林科技大学，2008:
21-29
[5]刘东，马永胜，付强.三江平原井灌水稻区地下水动态变化规律研究[J].灌溉排水学报，2006，25（5）:
42-46.
[6]付强，刘建禹，王立昆，等.基于人工神经网络的井灌水稻区地下水位预测[J].东北农业大学学报，2002，33
（2）:
152-159.
[7]NayakPC，SatyajiRaoYR，SudheerKP.Groundwaterlevelforecastinginashallowaquiferusingartificialneuralnetworkapproach[J].WaterResourcesManagement，2006，20
（1）:
77-90.
[8]付强，王立坤，梁川，等.GM（1，1）模型在预测地下水位动态变化中的应用[J].农机化研究，2001，（4）:
79-81，85.
[9]陈南祥，张海丰，李松海.基于混沌时间序列的地下水位多步预测模型[J].地球科学与环境学报，2007，29
（1）:
66-69.
[10]刘东，付强.小波随机耦合模型在三江平原井灌区地下水位预测中的应用[J].数学的实践与认识，2008，38（16）：
54-61.
[11]程正兴，杨守志.小波分析的理论算法、进展和应用[M].北京：
国防工业出版社，2007.
[12]付强．农业水土资源系统分析与综合评价[M]．北京：
中国水利水电出版社，2005．228—235．
[13]丁晶，刘权授．随机水文学[M]．北京：
中国水利水电出版社，1997．54—59．
源程序代码：
clc
t=1:
108;
Q（t）=[945.2676.2829.9832.3841.412761181.7136017762307201416082335160719711343971.6927.71051.3781.61839286716761655...
928.41040.2794.2580.3529.6516.1482.4674.2535.42497424451797798.4658.3766.4521.5721.61226959.91661.72435280716121941...
304.1205.7203.2180.5808.1870.463363179753675495523599740716117632427151687126310901212883.3866.9687.3709.5...
199531562638304928963093253333103816338919031443664.9433.4653.46885274.35197258121804721451228502113...
1031818175018481314962.47551201788.6402.8252.9266.8];
figure
plot（t，Q（t））;
title（'原始流量曲线'）;
t=13:
108;
Q1（t）=Q（t）-Q（t-12）;
figure
plot（t，Q1（t））;
标准化处理
a=mean（Q1）;
fort=13:
108
X（t）=（Q1（t）-a）/std（Q1，0，2）;
end
figure
t=13:
108;
plot（t，X（t））;
-----------------------------------------------
分解与重构
t=13:
108;
C0（t）=X（t）;
%fort=97:
106
%C0（t）=2*C0（t-1）-C0（t-2）;
%end
h=[1/161/43/81/41/16];
fort=13:
96
C1（t）=0;
fork=1:
5
C1（t）=C1（t）+h（k）*C0（t+2*k）;
end
W1（t）=C0（t）-C1（t）;
end
forj=1:
20
C1（96+j）=C1（97-j）;
end
%fort=97:
116
%C11（t）=C1（t-1）;
%C12（t）=2*C1（t-1）-C1（t-2）;
%C1（t）=（C11（t）+C12（t））/2;
%end
fort=13:
96
C2（t）=0;
fork=1:
5
C2（t）=C2（t）+h（k）*C1（t+4*k）;
end
W2（t）=C1（t）-C2（t）;
end
fort=13:
96
Y（t）=C2（t）+W1（t）+W2（t）;
end
t=13:
96;
figure
plot（t，W1（t））;
title（'分解曲线（a）'）;
t=13:
96;
figure
plot（t，W2（t））;
title（'分解曲线（b）'）;
t=13:
96;
figure
plot（t，C2（t））;
title（'分解曲线（c）'）;
t=13:
96;
figure
plot（t，Y（t））;
title（'重构曲线（d）'）;
%检验
fort=13:
96
X1（t）=Y（t）-X（t）;
end
%--------------------------------------------
%求W1（t）自相关系数
t=13:
96;
X2=mean（W1（t））;
X3=0;
fort=13:
96
X3=X3+（W1（t）-X2）^2;
end
fork=1:
21
X4（k）=0;
fort=13:
（96-k）
X4（k）=X4（k）+（W1（t）-X2）*（W1（t+k）-X2）;
end
r1（k）=X4（k）/X3;
end
%求W2（t）自相关系数
t=13:
96;
X5=mean（W2（t））;
X6=0;
fort=13:
96
X6=X6+（W2（t）-X5）^2;
end
fork=1:
21
X7（k）=0;
fort=13:
（96-k）
X7（k）=X7（k）+（W2（t）-X5）*（W2（t+k）-X5）;
end
r2（k）=X7（k）/X6;
end
%求C2（t）自相关系数
t=13:
96;
X8=mean（C2（t））;
X9=0;
fort=13:
96
X9=X9+（C2（t）-X8）^2;
end
fork=1:
21
X10（k）=0;
fort=13:
（96-k）
X10（k）=X10（k）+（C2（t）-X8）*（C2（t+k）-X8）;
end
r3（k）=X10（k）/X9;
end
fork=1:
21
y1（k）=（-1+1.96*sqrt（84-k-1））/（84-k）;
y2（k）=（-1-1.96*sqrt（84-k-1））/（84-k）;
end
figure
k=1:
21;
plot（k，r1（k），'-*'，k，r2（k），'-d'，k，r3（k），'-o'，k，y1（k），k，y2（k））;refline（0，0）;
legend（'r1（W1）'，'r2（W2）'，'r3（C2）'）;
%title（'各小波变换序列自相关图'）;
%--------------------------------------------------
%求W1（t）方差谱密度
k=1:
21;
D（k）=0.5+0.5*cos（pi*k/21）;
forj=1:
21
X11（j）=0;
fork=1:
21
X11（j）=X11（j）+D（k）*r1（k）*cos（2*pi*j*k/42）;
end
Sf1（j）=2*（1+2*X11（j））;
end
%求W2（t）方差谱密度
forj=1:
21
X12（j）=0;
fork=1:
21
X12（j）=X12（j）+D（k）*r2（k）*cos（2*pi*j*k/42）;
end
Sf2（j）=2*（1+2*X12（j））;
end
%求C2（t）方差谱密度
forj=1:
21
X13（j）=0;
fork=1:
21
X13（j）=X13（j）+D（k）*r3（k）*cos（2*pi*j*k/42）;
end
Sf3（j）=2*（1+2*X13（j））;
end
figure
j=1:
21;
plot（j，Sf1（j），'-*'，j，Sf2（j），'-d'，j，Sf3（j），'-o'）;
legend（'r1（W1）'，'r2（W2）'，'r3（C2）'）;
title（'·方差谱密度图'）;
%---------------------------------------------------
%W1（与W2（t）互相关分析
fork=1:
21
X14（k）=0;
fort=13:
（96-k）
X14（k）=X14（k）+（W1（t）-X2）*（W2（t+k）-X5）;
end
r4（k）=X14（k）/sqrt（X3*X6）;
end
%W2（t）与C2（t）互相关分析
fork=1:
21
X15（k）=0;
fort=13:
（96-k）
X15（k）=X15（k）+（W2（t）-X5）*（C2（t+k）-X8）;
end
r5（k）=X15（k）/sqrt（X6*X9）;
end
%C2（t）与W1（t）互相关分析
fork=1:
21
X16（k）=0;
fort=13:
（96-k）
X16（k）=X16（k）+（C2（t）-X8）*（W1（t+k）-X2）;
end
r6（k）=X16（k）/sqrt（X9*X3）;
end
fork=1:
21
r01（k）=1.9841/sqrt（1.9841^2+21-2）;
r02（k）=-1.9841/sqrt（1.9841^2+21-2）;
end
figure
k=1:
21;
plot（k，r4（k），'-*'，k，r5（k），'-d'，k，r6（k），'-o'，k，r01（k），k，r02（k））;
legend（'r4（W1，W2）'，'r5（W2，C2）'，'r6（C2，W1）'）;
%title（'各小波变换序列互相关图'）;
%----------------------------------------------------------------
%建立细节序列W1（t）的随机模型
t=13:
96;
X2=mean（W1（t））;
FW1=var（W1（t））;
CsW1=0;
fort=13:
96
CsW1=CsW1+（W1（t）-X2）^3/（（84-3）*FW1^（3/4））;
end
CsW1;
m1=eye（21）;
fori=1:
21
forj=1:
21
ifim1（i，j）=r1（j-i）;
end
ifi>j
m1（i，j）=r1（i-j）;
end
end
end
m1;
i=1:
21;
m2=r1（i）';
p1=inv（m1）*m2;
%绘制W1（t）偏相关系数图
fork=1:
21
y3（k）=1.96/sqrt（84）;
y4（k）=-1.96/sqrt（84）;
end
k=1:
21;
figure;
plot（k，p1（k），'-*'，k，y3（k），k，y4（k））;
re
升级会员 