基于小波随机耦合模型在流域流量方面应用.docx

1、基于小波随机耦合模型在流域流量方面应用基于小波随机耦合模型在流域流量方面的应用摘要：小波分析是一种时、频多分辨率分析方法，能对时间序列进行多尺度分析，它不仅可以充分展示水文时间序列的精细结构，而且还便于提取水文时间序列的变化规律。流域流量既是这样一个复杂的时间序列，其变化与区域气候条件和生态环境密切相关，具有明显的年际同期性和随机性波动，用小波随机耦合模型描述流域流量的动态变化具有重要意义和深远的影响。关键词：小波分析 ; 随机模型 ; 耦合预测；地下水埋深Based on wavelet random coupling model in the river flow applications

2、Abstract: Wavelet analysis is a kind of when multiresolution analysis method, the frequency, the time series of multiscale analysis, it not only can fully display in hydrological time series of fine structure, but also easy to extract the change rule of hydrological time series. River flow is such a

3、 complex time series, the change and regional climate conditions and ecological environment related closely, have obvious interannual TongQiXing and random fluctuation, using wavelet random coupling model to describe the dynamic change of flow river basin has important significance and far-reaching

4、influence.Key words: Wavelet analysis；Random model；Coupling forecast；Groundwater depth1 前言水文随机模拟已经取得了不断的进步和广泛应用。小波分析在水文水资源学中得到了广泛的研究和应用1。文献2-4将小波分析与目前成熟的系统分析方法耦合建立了水文预测模型，结果表明其途径是可行的。流域流量预测有多种方法，文献中5-9分别采用时间序列分析、人工神经网络、灰色理论等方法对流域流量进行了预测，但这些方法都属于单一预测方法，难以获得理想的预测效果。由于流域流量变化过程属于非平稳随机过程，流量的的大小受到多种复杂因素的影

5、响，因此，要对流量进行准确预测是非常困难的。近些年发展起来的小波分析方法可以充分展示水文时间序列的精细结构，便于提取水文时间序列的变化规律10。2小波分析的发展概述2.1 小波分析的历史回顾小波分析11是自1986年以来由于YMeyer、SMallat及IDaubechies等人的奠基工而迅速发展起来的新兴学科，它是Fourier分析划时代发展的结果。然而，它的展历史可以追溯到1909年Haar的工作。从现代小波分析的观点来看，1930年前后许多与小波有关的新方向出现，但以后进展不大，直到1960年Calder6n的研究及20后(1980年)Grossmann与Morlet的研究，后人称为“原

6、子分解”，才有了现代小波分的雏形。进而，由于1986年以后的工作以及应用的广泛性使这个学科得到了飞速发展。2.2 现代小波分析的发展阶段11。第一阶段为孤立研究与应用阶段。主要特征是得到了一些特殊的小波，并在一些专业领域得到零散应用。一个代表性的工作是法国地质学家AGrossmann和JMorlet于1984年第一次把“小波”用于处理地质数据，得到了以他们的名字命名的Grossmann-Morlet小波。另一个代表性的工作是1981年JStromberg与他的合作者发现的正交小波基。这个时期一个有趣的现象是各个领域的专家、学者和工程师在完全不了解别人的研究工作的状况下巧妙地、独立地构造自己的“

7、小波”。第二阶段为小波分析理论成型期。这个时期，法国数学家YMeyer成功地构造出了具有一定衰减性质的光滑函数，这个函数的二进伸缩和整数平移产生的函数组构成空间L2(R)的规范正交基。这引起了大家广泛关注，因为在此之前，学术界普遍认为不可能具有如此好性质的函数，Meyer也是认为这样好的函数不可能存在，在构造反例的过程中，把这个正交小波构造出来。随后，PLemarie和GBattle又分别独立地构造出同样“好”的小波。第三阶段为小波理论和应用全面发展期。一般认为，它是从20世纪90年代开始的。它的主要特征是，在上述理论框架下，出现了许多很有价值的应用成果，也解决了长期没有解决的应用问题，相反地

8、，在应用中也提出了许多需要解决的问题，从而推动了小波分析理论的发展。3 小波随机耦合模型建模思路首先将研究的水文时间序列采用快速小波变换算法(例如A Trous算法)进行小波分解，得到某尺度P下的小波变换序列W1(t)，W2(t) ，Wp(t)，Cp(t)；然后对各小波变换序列的主要成分(随机成分或确定成分)进行识别，对各小波变换序列进行互相关分析，并建立各小波变换序列适宜的数学模型；最后采用小波变换重构算法得到所研究水文时间序列的小波随机耦合模型。4 小波随机耦合模型在诺敏河流量变化预测中的应用根据诺敏河流域99-07年的实际观测资料，对该流域流量变化进行模拟，并预测未来发展趋势，为有关决策

9、者提供参考依据。4.1 实测月流量序列平稳化处理由图1可以看出，诺敏河流量序列Ht为一非平稳时间序列， 要使Ht转换为平稳时间序列，可以对其进行差分处理由于Ht在均值水平上不平稳，因此认为只需对其进行一次差分即可达到平稳，即：式中:B为后移算子通过差分，使得诺敏河流量序列Ht由非平稳序列转化为平稳序列Xt， 见图2。 图1 诺敏河原始流量变化曲线图2 差分后诺敏河原始流量变化曲线4.2 差分后逐月诺敏河流量序列小波分解与重构采用A Trous算法，取尺度数P=2，对差分后年降水标准化序列Xt进行分解，得到小波分解序列W1(t)，W2(t) ，Wp(t)，Cp(t)；将各小波分解序列进行叠加，得

10、到重构序列。重构过程与差分后序列Xt变化过程完全一致，因此，采用A Trous算法对诺敏河流量序列进行分解是可行的。图3 小波分解与重构4.3 小波变换序列成分识别采用公式和公式计算各小波变换序列的自相关系数和方差谱密度，并绘制自相关图和方差谱密度图，见图和图。根据结果图分析各小波变换序列的变化特性，可以近似认为W1(t)代表序列的随机项，W2(t)代表序列的周期项，C2(t)代表序列的趋势项。图4 各小波变换序列自相关图图5 小波变换序列方差谱密度图4.4 小波变换序列互相关分析采用下列公式分别计算各小波变换序列的互相关系数，绘制互相关图，并加绘95容许限。由图可以看出，各小波变换序列互相关

11、系数基本上落在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.03940，方差2W1= 0.1.8573，偏态系数Csw1=0.8684，趋近于0，所以认为序列W1(t)近似于正态分布，不必进行正态性转化。对序

12、列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，所以残差是

13、相互独立的假定是可以接受的。对独立随机序列t的正态性进行检验:均值-t=- 0.039430，方差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

14、准测，判定所需建立的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

15、(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张凯.关于天

16、津市水资源价值与水价格的研究D.天津:天津大学管理学院.2006.2姚雨田.天津市水资源优化配置与节水研究J.天津建设科技，1996，(4).3韩海燕，于苏俊.基于多目标优化模型的荣县水资源优化配置J.安徽农业科学，2010，38(19): 200-202.4张成凤.基于遗传算法的榆林市水资源优化配置的研究D.陕西:西北农林科技大学，2008:21-295刘东，马永胜，付强.三江平原井灌水稻区地下水动态变化规律研究J.灌溉排水学报，2006，25(5):42-46.6付强，刘建禹，王立昆，等.基于人工神经网络的井灌水稻区地下水位预测J.东北农业大学学报，2002，33(2): 152-159.

17、7 Nayak PC，Satyaji Rao YR，Sudheer KP. Groundwater level forecasting in a shallow aquifer using artificialneural network approachJ.Water Resources Management，2006，20(1):77-90.8付强，王立坤，梁川，等. GM(1，1)模型在预测地下水位动态变化中的应用J.农机化研究，2001，(4):79-81，85.9陈南祥，张海丰，李松海.基于混沌时间序列的地下水位多步预测模型J.地球科学与环境学报，2007，29(1):66-69.1

18、0刘东，付强.小波随机耦合模型在三江平原井灌区地下水位预测中的应用J.数学的实践与认识，2008，38（16）：54-61.11程正兴，杨守志.小波分析的理论算法、进展和应用M.北京：国防工业出版社，2007.12付强农业水土资源系统分析与综合评价M北京：中国水利水电出版社，200522823513丁晶，刘权授随机水文学M北京：中国水利水电出版社，19975459源程序代码：clct=1:108;Q(t)=945.2 676.2 829.9 832.3 841.4 1276 1181.7 1360 1776 2307 2014 1608 2335 1607 1971 1343 971.6 92

19、7.7 1051.3 781.6 1839 2867 1676 1655 . 928.4 1040.2 794.2 580.3 529.6 516.1 482.4 674.2 535.42 4974 2445 1797 798.4 658.3 766.4 521.5 721.6 1226 959.9 1661.7 2435 2807 1612 1941 . 304.1 205.7 203.2 180.5 808.1 870.4 6336 3179 7536 7549 5523 5997 4071 6117 6324 2715 1687 1263 1090 1212 883.3 866.9 68

20、7.3 709.5 . 1995 3156 2638 3049 2896 3093 2533 3310 3816 3389 1903 1443 664.9 433.4 653.4 688 5274.3 5197 2581 2180 4721 4512 2850 2113 . 1031 818 1750 1848 1314 962.4 755 1201 788.6 402.8 252.9 266.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)

21、;for t=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);% for t=97:106% C0(t)=2*C0(t-1)-C0(t-2);% end h=1/16 1/4 3/8 1/4 1/16;for t=13:96 C1(t)=0; for k=1:5 C1(t)=C1(t)+h(k)*C0(t+2*k); end W1(t)=C0(t)-C1(t);end for j=1:20 C1(96+j)=C1(97-j);end % for t=97:

22、116 % C11(t)=C1(t-1);% C12(t)=2*C1(t-1)-C1(t-2);% C1(t)=(C11(t)+C12(t)/2;% end for t=13:96 C2(t)=0; for k=1:5 C2(t)=C2(t)+h(k)*C1(t+4*k); end W2(t)=C1(t)-C2(t);end for t=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

23、; figure plot(t，C2(t); title(分解曲线(c); t=13:96; figure plot(t，Y(t); title(重构曲线(d);%检验for t=13:96 X1(t)=Y(t)-X(t);end% -% 求W1(t)自相关系数t=13:96;X2=mean(W1(t);X3=0;for t=13:96 X3=X3+(W1(t)-X2)2;endfor k=1:21 X4(k)=0; for t=13:(96-k) X4(k)=X4(k)+(W1(t)-X2)*(W1(t+k)-X2); end r1(k)=X4(k)/X3; end%求W2(t)自相关系数t

24、=13:96;X5=mean(W2(t);X6=0;for t=13:96 X6=X6+(W2(t)-X5)2;endfor k=1:21 X7(k)=0; for t=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;for t=13:96 X9=X9+(C2(t)-X8)2;endfor k=1:21 X10(k)=0; for t=13:(96-k) X10(k)=X10(k)+(C2(t)-X8)*(C2(t+k)-X8)

25、; end r3(k)=X10(k)/X9;endfor k=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);for j=1:21

26、X11(j)=0; for k=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)方差谱密度for j=1:21 X12(j)=0; for k=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)方差谱密度for j=1:21 X13(j)=0; for k=1:21 X13(j)=X13(j)+D(k)*r3(k)*cos(2*pi*j*k/42); end S

27、f3(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)互相关分析for k=1:21 X14(k)=0; for t=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)互相关分析for k=1:21 X15(k)=0; for t=13:

28、(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)互相关分析for k=1:21 X16(k)=0; for t=13:(96-k) X16(k)=X16(k)+(C2(t)-X8)*(W1(t+k)-X2); end r6(k)=X16(k)/sqrt(X9*X3);end for k=1:21 r01(k)=1.9841/sqrt(1.98412+21-2); r02(k)=-1.9841/sqrt(1.98412+21-2);endfigurek=1:21;

29、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;for t=13:96 CsW1=CsW1+(W1(t)-X2)3/(84-3)*FW1(3/4);endCsW1;m1=eye(21);for i=1:21 for j=1:21 if ij m1(i，j)=r1(i-j); end endendm1;i=1:21;m2=r1(i);p1=inv(m1)*m2; %绘制W1(t)偏相关系数图 for k=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