水文随机分析课件

2024/3/121水文随机分析

2024/3/122第一章随机过程基础知识随机过程概念随机过程的概率分布及数字特征随机过程的基本分类平稳随机过程泊松过程2024/3/123随机过程概念实际上，常遇到实验过程中随某个参变量变化而变化的随机变量，数学上称该随机变量为随机函数。（随机变量定义：随机事件的实数值函数，有一个基本事件，对应一个实数值，这个实数在一次试验中能否发生，是很难事先确定的）。2024/3/124如南京滁河某站水位或流量，它是随时间而变化的，包括年平均流量，年最大流量，日、月平均流量或水位都是随时间而变。再如南京日、月平均气温值也随时间而变化，当然南京气温值还随空间位置不同而变化。换句话说，参变量不一定总是时间，可以是其他。这些随机变量即为随机函数。特别是：我们常称以时间t为参变量的随机函数为随机过程。当然如果涉及水文现象的随机过程则称为水文随机过程，一般用表示。

2024/3/125t可以是连续的，也可以是离散的。如t为离散的则简称随机序列或时间序列，如年最大等，如t为连续的则仍称随机过程。在给定t情况下，就是一个随机变量，其取值可以是离散的，也可以是连续的。（水文上一般是连续型的，年径流量、年最高水位等；也有离散的，如年降水天数，1，2，…，365）2024/3/126随机过程的分类2024/3/127在研究水文随机现象时，如研究洪水过程，如果把t当作连续的（瞬时过程），理论上讲是最好的，但这样的随机过程建模对资料要求高，工作量很大，实际上几乎难于实现。因此，常根据实际水文现象特性对t作离散化处理，如对大江大河洪水过程，不要求t为连续的，而只要日平均过程，即一年365个数据则可，经验表明：日平均过程可以近似反映长江干流洪水。当然对中小河洪水过程，则不能用日平均流量来反映，而应该取时段长为几个小时的平均流量做离散化（山区河流因流量变化大则应更短，一次洪水过程可用1、2个小时平均流量过程代替洪水过程）。2024/3/128在任意给定一个t值，为随机变量，既然是随机变量，那么就会有概率分布。，对于时刻t1，对于时刻t1和t2

。

，n个时刻的联合分布。

由于研究多维联合分布难度很大，因此常需要研究随机过程的数字特征，一般到2阶即可。二、随机过程的概率分布与数字特征2024/3/129随机过程的主要数字特征

2、方差()

3、自相关系数1、数学期望2024/3/12101、平稳过程与非平稳过程。主要看随机过程的统计特性是否随时间变化分类。如年径流或年降水过程在人类活动影响很小时可以认为是平稳随机过程，但洪水过程不是平稳的。2、独立随机过程与非独立随机过程主要看各时刻状态之间是否相互独立。年最大洪峰流量过程为独立随机过程，而日流量过程则为非独立随机过程。其中有一种特殊过程：Markov（马尔科夫过程）如AR（1）过程，是非独立随机过程里常见的一种，实际上应用此较多。即将来状态与现在有关，而与其前面状态毫无关系。10三、随机过程的基本分类2024/3/1211

定义

如果对于时间t的任意n个数值t1,…,tn和任意实数k，随机过程的n维分布函数满足关系式则称为平稳随机过程。11当n=1

时，不管K取何值，；任何一维分布都是用同分布。四、平稳随机过程2024/3/1212当n=2时，说明二元联合分布仅与这两个随机变量时间间隔有关，与取多少无关。而且还可以证明，。这种平稳过程，称为严平稳过程。由于实际上要求得多元联合分布难度很大，应该说绝大多数情况下是不可能办到的，因此，一般只要求关心随机过程一、二阶矩。当，均与无关，这时可称随机过程为宽平稳过程。今后所说平稳过程，一般都是指宽平稳过程。2024/3/1213②宽平稳过程自相关函数性质≤1，。如

，则≥0。b.，偶函数。③平稳随机过程各态历经性设为平稳的一个样本或一个现实，令为一个样本或现实的平均值。当，则称平稳随机过程具有各态历经性。2024/3/1214事实上，这个概念还是很重要的。在许多实际随机过程特别是水文过程中，仅能有一个样本或一个现实，那么要估计不同t下，(某一个时刻t仅一个数据)则无法进行，这时如果具备历经性，则可用代替。各年年径流随机过程

这些资料仅有一个样本或现实。2024/3/1215以上用代替，满足两个条件n足够大，太小了不能反映实际，误差大，同时满足各态历经性。2024/3/1216五、泊松过程①独立增量过程若随机过程X(t)(t≥0)满足条件a:b:对任意时刻，（任意给定n+1时刻），如果过程增量相互独立，则称X(t)为独立增量过程。可以证明它是Markov过程。②泊松过程泊松分布：n次独立试验中A事件发生了k次2024/3/1217独立增量过程X(t)，若其增量的频率分布为泊松分布，t2>t1≥0,n=0,1,2,…则称X(t)为泊松过程。2024/3/1218当t1=0,t2=t，则，为平均数值（随机变量）

对于这种随机过程可以在水文中描述，（0，t）时间间隔内出现降水次数的概率。（在一段时间内接听电话的次数也可按泊松分布）2024/3/1219例如，某站汛期[0，t]（t=30,30天内）年平均暴雨发生次数为4.8次，即=4.8，这个数值完全可以根据实际观测次数系列求平均得到，则在汛期开始30天内发生n次暴雨的概率。

这样就可以知道发生不同次数暴雨的概率，如发生3-5次概率50%。这对于防汛决策是有意义的。当然这是假定符合泊松过程为前提。如果要检验它是否正确，需要有大量资料。如果可找到全国各地汛期暴雨资料，可以分析我们在汛期内发生暴雨次数是否符合这一泊松过程规律。2024/3/1220第二章水文时间序列的组成分析概述趋势项分析处理与检验跳跃成分分析处理与检验周期成分的描述与提取2024/3/1221概述趋势S（t)2024/3/1222对系列组成成分分析的目的随机模拟，包括年月径流系列，洪水系列(n年)，如长江中下游地区洪水模拟可以计算三峡水库防洪效益的计算误差。

预测（外延），西北干旱区出口径流作2年预测，上海水情中长期预测，晋江水情预测。

2024/3/1223一般的水文随机序列，非周期[Tt（趋势），Ct（突变），跳跃Kt等]周期[简单或复合周期成分]随机成分

当，则仅为随机成分，这时要作预测？比有确定性成分时难度大，往往不易预测准确。反之，如S（t)=0,则系列仅有确定性成分，这时候一般易于进行预测，但如果规律性发生变化，则预测会出问题，如树与小孩生长例子。2024/3/1224当，则，这时只要趋势及周期等成分分析的比较好，则预测就比较可靠。当然如果这种确定性成分在今后时期内不是按现有资料变化规律变化，那么预测起来精度也是不能保证的。要注意的是有些情况下，2024/3/1225第二节趋势项分析处理与检验对于一个时间序列，随着时间增长，呈现出系统而连续的增加或减少的变化。这种有规则变化称为趋势。这种趋势可分整体趋势与局部趋势，往往是由于人为或自然原因造成，而不是随机抽样波动或观测资料误差所致。如气候因素年际变化，若有某种明显趋势，那么年降水量序列和年径流序列则有可能出现相应的趋势。在一定时期内湖泊中泥沙逐渐淤积，年平均水位就有升高的趋势；如沿河逐年提高防洪堤，年最大洪峰流量就有增大趋势等。流域内，灌溉面积不断增加，流域蒸发量有增加趋势，当然径流量就有减少趋势。为了排除趋势成分应从物理成因和统计分析两个方面着手进行。2024/3/1226查明趋势现象及其产生原因，然后使用数学方法加予描述，进而加于排除。趋势变化可以是线性和非线性的，常用多项式来描述。趋势项剩余项为系数，一般实际上先用简单线性模型来描述如何定？2024/3/1227利用现有观测数据，使得残差最小为原则。

,来优化确定系数。如果是一元线性，仅有两个参数。优化计算在许多方面都有应用，如计算机优化适线，回归分析，最小二乘法、流域模型参数优化计算等2024/3/1228分离趋势方法还有差分法，如：见下一个片子趋势项目检查（验）方法：①移动平均法（把原系列从变幅大→变幅小，容易看出是否有趋势）若观测值为，若移动平均区段为h，例h=3则移动平均值其中为权重系，，当相等则是简单的算术平均，当不等时，则为加权平均。

2024/3/12292024/3/1230得到新系列，（比原来少2项）一般通式：移动平均值y为h为奇数时，为系数，当则为算术平均，还可加权取值，

2024/3/1231h为偶数，移动平均后仅有n-h个数据。H的选择一般根据时间序列的周期来选择，即区段长度等于周期长度，如对月平均径流量时间序列，由于存在年周期，h=12

2024/3/1232②Kendall秩次相关检验（水资源综合规划中使用此方法分析降水和径流变化趋势）

对于序列，先确定所有对偶值中的出现个数（设为P），顺序的子集是：

如果按顺序前进的值全部大于前一个值，是一种上升趋势，则，系为等差级数，总和为2024/3/1233如果序列全部倒过来，即由大到小排列，则P=0，为下降趋势，（说明P数值不能太多，也不能太少）P的数学期望（值）构造统计量当

2024/3/1234：无趋势，当，接受无趋势，：无趋势

，拒绝原假设。

例12个数据，看是否有趋势，Xt509490475482513535498540550560545530789864531010故趋势明显，有上升趋势。2024/3/1235③Spearman秩次相关检验分析序列Xt与时序t的相关关系，在运算时，Xt用其秩次Rt（即把序列Xt从大到小排列时，Xt所对应的原来序列中的序号），t仍为时序秩次相关系数：n为序列长度，，显然如秩次与t序号相近时则d小，秩次相关系数大，接近1，趋势显著。若反过来，则达最大，r接近-1。r是否异于0，用t检验。：无趋势。2024/3/1236

R=-0.471T=-3.12

趋势显著

也可以建立线性回归方程检验线性趋势是否明显，当然还要以肉眼观看其变化趋势，这是最直观的。

例子：数据与上面同t123456789101112数据509490475482513535498540550560545590排序560550545540535530513509498490482475Rt109118612517243Dt=Rt-t978416-2-7-2-8-7-19d2814964161364494644981

R=-0.471T=-3.12

2024/3/1237第三节：跳跃成分分析处理与检验跳跃是指水文系列急剧变化的一种形式，当水文序列从一种状态过渡到另一种状态时表现出来。平稳过程（序列）跳跃大小跳跃一般也出现序列均值、方差与自相关系数等参数之中，实际上多在均值中寻找跳跃。2024/3/1238跳跃是由于人为或天然原因造成，如修筑水库前的坝下年最大流量序列与修建水库后经过水库调节后的年最大流量序列，就是人为引起的跳跃，修造后均值与方差变小（这是事实）。又因为修建水库增加水的面积，蒸发增加，可能引起下游年径流量均值的跳跃。如尼罗河阿斯旺坝断面年径流系列。

为从水文序列中排除跳跃成分，也和趋势分析一样，应先查明跳跃现象产生原因，并进行数学描述，再加排除。2024/3/1239跳跃成分是否存在的检验：检验方法：分成两个样本，假设前面样本分布，后面样本分布

若拒绝，认为总体发生显著变化，跳跃显著；若接受，认为总体发不发显著变化，跳跃不明显。关键，如何定：

①调查流域自然地理条件变化，确定因自然或人为原因使序列发生显著变化的时间②用时序累积值相关曲线法确定③统计推断确定，也可以目估看变化2024/3/1240时序累积值曲线法：设确定序列，参证序列（无跳跃和趋势突变），两个序列的累积值分别为点绘关系图（如右图）如果研究序列跳跃不显著，则为一条通过原点的直线，否则为折线，转折点即为。图中，1956年为研究序列的一个跳跃点，当然选择参证序列时，参证序列不应包含有暂态（趋势、跳跃、突变等成分2024/3/1241分布一致性检验，①秩和检验法假定前后两个分布密度为和，从总体中取样本长度分别为，方法：将两个样本所有数据依小→大排列并统一编号，规定每个数据在排列中所对应的序数称为该数的秩，对于相同的数值，则用它的序数的平均值作秩。现记容量小的样本各数值的秩之和为W（统计量），秩的检验就是对W作检验（W太大或太小都证明总体前后不一致）：当时，统计量W近似于正态分布，2024/3/1242于是可用U检验~N（0，1）——小样本容量——大样本容量接受，否则拒绝。2024/3/1243例子：对所有数据由小→大排，所排序号为该数据的秩，把数据容量小的，如本例把个数据对应秩累加起来为W，再求U。W太大或太小,表明有跳跃成份。t12345678910111213Xt250210230275220245221265247220250205215t14151617181920212223Xt231202206209218204209219202214Xt*202202204205206209209210214215218219220wi22345778910111213Xt*220221230231245247250250265275wi141516171819202122232024/3/1244

②游程检验法有一观测值序列：11，9，7，12，14，15，16，10，13，假设分割点,将它们由小到大的顺序进行排列。7，9，10，11，12，13，14，15，16把属于的记为A，的记为B

，这样得到一个新的序列。，把每一连续出现同一字母的称为1个游程，每个游程所含元素的个数为游程长，例如上式中，有2个A游程，游程长度为2；有2个B游程，游程长度分别为1和4，全部游程数为4个。2024/3/1245

当游程出现个数较期望的游程少时，就比较趋向于拒绝两个样本来自同一总体的假设。因为此时长的游程较多,表明个别样本中的元素有较大聚集现象，因此,我们认为不服从同一总体，这是游程检验的指导思想。2024/3/1246游程个数检验法：当游程总个数K服从正态分布跳跃不显著2024/3/1247K太多、太少都说明有跳跃上例计算游程数K=6U=2.77

拒绝有跳跃成分。以上秩和检验与游程检验均属于非参数假设检验，其他还有参数检验，如方差、均值齐性检验等。当经过成因分析和统计推断，水文系列中有趋势或跳跃明显时，可用适当方法加以描述，再从序列中排除掉，剩余的就是具有原始状态或一致条件。2024/3/1248突变成分检测：由于人为或自然原因，水文序列中可能出现突变。这一变化一过就又恢复原状，这可看作是跳跃的一个特殊情况，如由于塌方拦截江河，形成水库，以后又溃坝，这就引起流量突变，但临时水坝冲毁，又恢复原来状态。2024/3/1249第四节：周期成分的描述与提取一、周期成因周期地球绕太阳公转：对于月径流、旬径流、日平均径流存在年周期，周期分别是12个月、36旬和365天。对月、旬、日降水量、蒸发量一样存在着年周期。地球自转：气温、蒸发量等存在24小时周期。近似周期月球绕地球旋转影响，潮水位过程出现周期，但频率不可通约。太阳黑子影响等影响，存在若干年为一个周期（近似）可能不均匀变化周期。大气环流2024/3/1250如黄河上游和松花江相邻两年年径流之和与前一年太阳黑子数有对应关系，因为太阳黑子有一定循环周期，因而年径流多年变化中也可能存在一定循环周期。再如1960年代分析结果，长江汉口站最大流量多年变化有55年主要周期，1989-90年代作分析，长江宜昌站100年资料（1881-1980）分析汛期流量存在15年主周期。2024/3/1251二、周期成分描述及谱分析

存在简单周期成分，（同正弦或余弦表示）—振幅—角频率当时，周期，∴周期对于月径流序列，应该存在T=12这一情况。—相位这个正弦函数，均值等于0，方差这个有周期过程的自相关函数表现出周期性质或特点。

具有复合周期的周期成分，可用以下公式（傅立叶级数表示）。2024/3/1252为基本周期，这个式子说明序列Xt可由l个谐波（周期）线性叠加而成，而且不同谐波周期与基本周期是倍比关系。

——为Xt

平均值。2024/3/1253那么如何根据已有实测序列求，下面先介绍谱分析技术方法---方差线谱。满足一定条件时，进行傅立叶级数展开，用上面式子表示。l为谐波总个数（n为偶数，）谐波振幅。角频率（为基本角频率，令f=1/n）2024/3/1254分析一下：在已知前提下，关系一一对应，不同角频率下谐波振幅，反映谐波大小。可以证明：为对应谐波方差大小。反映不同频率谐波所占比重，愈大，谐波周期含量愈大，愈显著。方差谱密度2024/3/1255代表任一点上对应的方差密度为方便定义，方差谱密度函数可以证明实用时

，（对称函数）从图形中要以看出哪一个频率下谱密度大，即该周期谐波分量比较大。可以根据实测样本估计出来（仍是离散化）2024/3/1256计算时m取值：，但自相关系数求相关系数（纠偏系数，对估计量而言）从估计结果可以知道哪些周期谐波比重大，哪些小不显著。2024/3/1257三、用简单分波法分析系列周期（近似）

—随机水文系列，—趋势，跳跃，突变—不含非周期确定性成分序列（趋势项或跳跃项）简单分波法寻找后序列中周期项，把（一般即可）。一般针对年、月径流，年、月降水量序列2024/3/12581、分析思路：①先分离第一周期，得出余波，计算余波系列均方差。②对余波A1分离第二周期，得出余波，计算余波系列均方差（显然会愈来愈少）。③对余波分离第三周期（如果存在），得余波（如果已无第二周期那么就不要做第③

），计算。④再对余波分离第四周期（如存在），得，计算（如果已无第三周期就不再做第④

）。2024/3/12592、分析周期用途：①了解序列未来变化规律②用于预测未来年份或月份的水文特征如对年径流序列分析下来有2个周期：第一个7年周期，其周期值已知，P1674.2P2658.2P3426.2P4520.4703.1P5751.4P6907.5P7第二个3年周期，Q1Q2Q392.7-22.9-78.6此时，可以用于预测，先假设无趋势项，即T=02024/3/1260(预则值)，作为第n+1年的周期项取值，再加T(n+1)（如果有趋势的话）。则：2024/3/12613、第一周期分离与提取要对n个数据分析周期，显然根据这些数据分析出来同期不可能超过（偶数），（奇数）。∴周期长度为：对于任意l年周期成分是否显著，。①计算整个系列均值②把现有时间序列分组（分成m行，l组，各组数据至多m个）如不为整数，取行。2024/3/1262按l年分组即周期长为l年情况对于水文时间序列数据，按时间顺序，从左→右，从上→下排列。（先取l=2作周期显著性检验，再取l=3，再按次序来作）2024/3/1263③计算每一个组计算均值，，及，，每一组数据个数，最多为。④计算组间离差平方和，

或计算组内离差平方和，⑤：不存在长度为l的周期成分，：当显著性水年，可查F表，当由样本计算则原假设成立，当，则拒绝原假设即存在明显的长度为l的周期成分。2024/3/1264如何理解？假如一组数据完全按周期排列，即每7年重复同一套数据（7个数据不同差异较大），在这种情况下，按分组，计算大小则F很大，存在周期但是如按分组（打乱了），小大F很小，周期不显著实际做的时候要计算，或2024/3/1265年数据显然在时，周期显著，其他情况下不显著，有时可能出同几个F计算值超过，这时应选哪个？选F值最大的。2024/3/1266⑥如何提取出这个7年周期值就是在时，计算的作为周期为7年周期成分，⑦求余波⑧求余波均方差2024/3/12674、第二周期分离与提取这时分析序列的周期成分，不再是原序列，而是对余波进行，余波与原水文序列差别主要在，一般水文序列不为负值，但余波有正有负，正负相抵消。分析周期方法与第一周期分析完全一样，只不过这时分组组数l应不包括第一周期长度（如第一周期长度7年，则这时对l

=7可不进行周期显著性检验）。2024/3/1268求余波求余波（残差）均方差，比第一余波均方差要小。假如经过检验后发现时明显存在周期，这时对余波序列，求时分组平均值，如2024/3/12695、第三、四周期分离与提取

与第二周期完全一样，不再叙述。当然，分析周期成分个数愈多，最后余波均方差会愈小，但由于实际上资料较短，本身存在误差，不一定愈多愈好。（特别用于预报时）原则上，为了提高周期分析精度，应使资料尽可能长一些。当增加新的资料，则应重新分析周期成分。（哪怕只增加一年或一次）另外显著性水平大小选取对周期成分判断也会有影响，所做假设检验前提是系列符合正态分布，事实上也未必成立。这样分析出来的周期成分会犯错误。对周期成分还可以通过小波分析，最大熵估计等方法。2024/3/1270四、平稳随机成分的分析

当前面介绍非周期成分及周期成分被提取之后，所得序列随机成分原序列非周期周期一般可认为是平稳随机过程。对于这样一个平稳随机成分，一般情况下，还可以分成相依随机成分+独立成分。对于这种平稳随机成分，可以用不同随机模型：线性、非线性（模型加以描述）。即平稳相依成分纯随机——不同t之间相互独立，平稳独立随机过程2024/3/1271即方差为为常数，但不同时刻t时相互独立。在水文计算中，主要研究的就是平稳独立随机序列，如年最大洪峰或洪量序列、年最大t天暴雨量序列。当中无相依成分，则为纯随机序列。对于纯随机序列，其分布线型（我国一般）：对数正态，P-III等，有时为了简化，要用正态分布，当然对年最高水位还采用极值分布。如何检验一个随机序列是平稳独立随机序列？2024/3/1272一般是通过计算自相关系数作检验根据不同k做计算，看是否落在上下置信项之间。若都落在置信项之间，则认为是独立的，若不在之间，那么就是有相关。事实上，我们处理一个无确定性成分的年径流序列，它一般并不是独立，而是有一定相关关系，但年最大洪峰、洪量序列一般可认为是独立平稳序列。因此把年径流序列当作纯随机成分作频率计算会带来一些误差。2024/3/1273对于一个平稳随机序列，若不是独立的，则应该把该序列分成两部分：相依部分与纯随机部分。对相依部分（成分）可以用下面介绍的平稳随机模型加以描述。有些同学未学习过随机水文学课程，故这里把几种常见纯随机变量的随机数生成复习一下。2024/3/1274①[0，1]均匀分布生成，在FORTRAN语言下，调用RANDOM(u)即可生成一个[0，1]均匀分布随机数，如需生成100个[0，1]均匀分布随机数u。

DIMENSIONA（100）

DO10I=1，10010CALLRANDOM（A(I)）

WRITE(*,*)ASTOPENDVB语言用内部函数RND生成。2024/3/1275②正态分布生成生成出这样生成……是相互独立随机数

，即符合均值为0，方差1正态分布。对于一般正态分布，把转换成……2024/3/1276③对数正态分布生成X----三参数对数正态生成生成出这样生成出来符合对数正态分布，其平均值为，均方差，偏态系数。2024/3/1277生成un+1,un+2,un+3NoYesz’=x+yz=z’/+a0生成u1,u2,…,unP-III分布随机数的生成EX,Cv,Cs2024/3/1278举例：EX=103.4,CV=0.482,CS=1.621,生成的随机数有:u1=0.86515，u2=0.69186，u3=0.41686,u4=0.86122要求生成一个P-III分布随机数.(1)

n=1p=0.522(2)生成y

(3)生成x,由u2，u3，u4

成立2024/3/1279

不同P情况下，抽样效率不一样

P愈小，效率愈高。

(4)问题:如何生成第二个P-III随机数?2024/3/1280第三章水文随机模拟第一节概述第二节线性平稳水文序列模型（ARMA(p,q)）第三节分数高斯噪声模型第四节季节性随机模型第五节散粒噪声模型第六节多变量模型2024/3/1281第一节概述两种模拟水文模拟：确定模型结构+估计参数→短期水文预报

水文随机模拟：通过随机数生成，得到很长水文模拟系列，可以充分利用已有序列信息。主要应用:理论研究（参数估计方法比较），风险与可靠性分析。

2024/3/1282主要简要介绍：ARMA(p,q)及其简化模型AR(p)，季节性AR(p)，多站AR(1)。2024/3/1283不同水文变量采用模型

对于平稳的年径流，年降水量过程，一般采用ARMA(p,q)模型，特别是AR(P)模型加以描述。但如有长持续性特征，应用分数高斯噪声模型。对于月、日旬径流或降水过程：也可以做平稳化变换（中心化和标准化变换）变成平稳过程，但如果以上变换后仍是不平稳的，则采有非平稳模型季节性模型，如月径流或日径流过程，可用季节性AR(1)，AR(2)模型描述，解集模型。洪水过程：对于特大江河一般可采用平稳模型。（长江干流，汉江，长江上游地区，珠江流域），对于一般河流，可采用季节性AR(P)模型、散粒噪声、解集模型等（山区性河流涨落明显）2024/3/1284水文随机模拟一般程序实测资料选样与审查（有时包括调查雨洪资料）水文随机序列组成分析与处理随机模型选择模型参数估计结束水文随机序列生成模型残差项独立计算检验有时还有正态性检验模拟序列检验包括Box—Cox，对数变换等(有些模型需要序列是正态)有时还要定阶不合理合理模型检验不好2024/3/1285第二节线性平稳水文序列模型（ARMA(p,q)）一、模型形式一般形式（直接形式）（离均差） 研究变量Xt，中心化变量,标准化变量

。 自回归模型AR(P)：

是p个参数，也是参数，另一个是的方差，这个模型中与相互独立，但与有关，p为阶数。2024/3/1286

显然这个模型告诉大家的是，t时刻值与前面p个时刻取值有关且是线性之和，这在水文中有现实意义，如年、月径流量大小是与前面时刻取值大小有关的。 表示残差，不同时刻t的残差要求独立同分布（但不一定是正态分布）。这一条件要求较高，所以水文变量不一定都能满足这一要求，即白噪声的要求，也是上面一般工作程序中所要做的检验内容之一。如果出现这种情况，应该把残差(假如用表示)，用以下形式表示：,q为阶数这个模型称为自回归滑动平均模型，ARMA(p,q),共p+q+2个参数.2024/3/1287 该模型能表征许多水文变量的变化特征，因此使用比较多。当，自回归模型。当，滑动平均模型。对于用表示 对于用表示

为的均方差。这三种形式都可使用。2024/3/1288传递形式 对AR(1)模型

,不失一般性，令

逐次用代入上式 这说明由无数个独立白噪声线性叠加而成的，是其权重，时间距离愈远，权重愈小。2024/3/1289如用Green函数表示，（对AR(1)模型而言）对模型（一般ARMA模型）传递形式：

2024/3/1290逆转形式 即把用表示形式 如对MA(1)模型

2024/3/1291也有逆转形式这两种形式在作水文时间序列预测时都用到。2024/3/1292二、自回归模型定阶与参数估计AR(P)

[对公式可用乘以方程两边求数学期望， 经过整理可得出]矩法参数估计：（P个方程，P个未知数）

2024/3/1293有了样本那么可以估计出参数。(k阶自相关系数)这样方程中参数均可估计出

不是真正的相容估计，模型定阶AR(P)2024/3/1294BIC定出的阶数P一般要比AIC小，是相合估计量，一致性估计量P上限值一般取为或。对于以上参数估计公式k阶自相关系数当P=1，2024/3/1295 对AR（1）模型，当为正数，按指数衰减，一直不为0（理论上）。因此，如果实际上计算出来自相关系数，若能符合上述要求（指数衰减要求）的话，那么可以认为用AR（1）模型来描述。AR（P）平稳性条件，方程根如都在单位园内，即，则是平稳的，否则模型不平稳。2024/3/1296三、一阶自回归模型建模与模拟一阶自回归模型，一阶马尔科夫模型最常用。参数估计2024/3/1297平稳性条件，要求在单位圆内即，一般都能满足。自相关系数有正有负按指数衰减①k∴不单单是图①情况下才可以用AR（1）模型加予描述2024/3/1298随机模拟：正态残差若用,则要模拟水文时间序列，必须生成白噪声序列：由[0,1]均匀分布随机数，生成生成时从t=1开始，或任意原序列中的数值2024/3/1299要生成水文序列，如年径流等

一般为了消除初值影响，前50~100次去掉，甚至可以更多一点，显然这样生成序列符合正态分布。问题将会出现，即出现负值，这是不符水文变化特性的。….2024/3/12100 解决方法： ①把原序列先变成正态序列模拟正态序列

同上述方程建模，包括定阶P及估计参数，这样即可生成很长y的模拟序列

至于完全可以利用，转换在统计书上都有（随机水文学）,上一章已经讲过。

2024/3/12101②把残差项改成,分布均值0，方差1，其偏态系数为一阶自相关系数,

偏态系数

这样生成序列只能说近似符合P-III分布分布用舍选抽样生成2024/3/12102四、二阶自回归模型建模与模拟AR（2）r平稳域要求(主要要求)模型形式2024/3/12103随机模拟先假定为原序列中最后两个数值,从t=2开始去掉前50~100次,这是正态过程2024/3/12104非正态过程与AR（1）不一样,残差P-III分布t-1t-12024/3/12105五、滑动平均模型 水文上用得少（中心化,即减去均值后） 该过程肯定是平稳的，但是否可逆，需要满足以下条件（其参数估计在参数估计中做介绍） 若方程根都在单位园内，则MA(q)系列满足可逆性条件。

2024/3/12106

自相关系数是截尾的模型估计系数

2024/3/12107六、ARMA(p,q)模型主要特点：比AR(p)和MA(q)能更好地反映水文变量在时序上变化统计特性，具有更大弹性。在达到一定的要求下，较AR(p)和MA(q)具有更少参数。但ARMA(p,q)模型参数估计较复杂。2024/3/12108

可逆平稳条件：

根要求在单位圆内则平稳。

根要求在单位圆内则可逆。

2024/3/12109

参数估计①对系列(中心化系列),按自回归模型估计参数(矩法)②从原始系列yt中减去具有参数的AR（P）系列③用矩法滑动平均模型估算参数参数矩法估计较粗，精度不高也可以用最小二乘法估计t-12024/3/12110

不同参数，使得

所对应的值，即为最小二乘估计。（该方法精度较高）

ARMA（1,1）为较常用不做详细介绍（有兴趣自己看书）

实例：用AR（p）模拟年径流过程用AR（p）模型模拟大渡河下游水文站铜街子控制站的年径流过程，F=7.64万km2，资料1937～1979年，铜街子站附近拟建一个水电站，现需对未来年径流变化作分析，资料是一致，可靠及有代表性（无趋势，无明显周期存在）

2024/3/12111

从年径流过程线上看，无趋势存在,从计算的自相关系数看,有明显相依性，考虑用AR(p)来加以模拟，用修正纠编矩法计算

考虑到系列偏态特征明显,建模时应考虑其偏态特性。

先对序列建立AR（1）模型

2024/3/12112是平稳的模型t=1生成随机数t=2生成随机数,，代入模型得…2024/3/12113

当然在生成年径流序列之前还应检验模型是否符合独立性前提条件

对做独立性检验，计算出

2024/3/12114作独立性检验，查，具备零相关。

∴符合建模前提条件，即残差具有独立性。年径流模拟成果检验

AR(1)2-3阶相关系数差别较大14890.1280.3330.4250.1810.075长系列14910.1280.3310.420.3020.206实测值14890.1260.2880.3720.1180.008短系列AR(2)14890.1290.3370.4260.3040.172长系列148914890.1250.3540.2180.076短系列2024/3/12115用AIC（p）定阶进一步从理论上证明P=2是最佳的。实用性检验：短序列检验即生成K个43年资料，然后计算各个43年序列的。对K个43年系列求所列参数, 的平均值与均方差。以正负一个均方差为检验可接受区间。P0123AIC451.7445.3446448.2BIC2024/3/12116第三节分数高斯噪声模型 用ARMA（p，q）仅能模拟短持续水文过程，当水文随机过程表现出长持续性（连续丰水年组或枯水年组交替出现）在相关系数图缓慢衰减时，此时应使用分数高斯噪声模型。 （类似滑动平均）独立的正态分布。M—记忆长度，取值较大。2024/3/12117Hurst系数K计算若有一个系列2024/3/12118 即把每年与平均值之差累加起来，相当于水库按平均流量泄水，水库累计水量。相当于调节库容（极差）

Hurst建议：Vn为Xt的标准差对于ARMA模型，n无穷大时,K=0.52024/3/12119 利用以上k计算公式，对800个系列（河川径流，降水量，气温，树木年轮等）计算k的数值。 发现 但对于ARMA(p、q)，当 这与Hurst计算实际系列不一致，出现矛盾，这种现象叫赫斯特现象。（原因：样本序列不够长或非平稳式序列或具有很长相关结构）。2024/3/12120

实际工作中，当n逐渐加大，k较快趋于0.5，可认为水文序列为短相关结构模型，可用ARMA（p，q）模拟；而当n逐渐加大，k较慢趋于0.5或不趋于0.5，则可把水文序列当作长持续模型。 由于该模型缺乏物理基础，因此实际上未得到普遍应用。2024/3/12121第四节季节性随机模型一、季节性AR（1）模型是一种比较实用模型，除了模拟洪水进程外，还可以用于模拟月径流过程。1、模型形式，第年份，

t为截口， 当表示月径流序列时，形式：

2024/3/12122—第t个截口水文序列均值—均值0，方差1，的P—Ⅲ分布随机数（独立）—为待定系数，不随而变，但随t而变。∴称为季节性模型，因为模型参数在年内不同月份间是变化的，随t而变。2024/3/121232、参数估计

为第t个截口的均方差，为第t-1个截口的均方差，为第t-1与第t个截口的相关系数，Cst为第t个截口偏态系数。2024/3/121243、水文序列的模拟逐年进行模拟。每年生成步骤：先第一截口即x1,1生成，可用纯随机模型生成。如已知，则按P-Ⅲ生成方法生成。第二截口生成：

其余截口类似于第2截口，直至最后一个截口这是生成第一年，对于第2年生成从头开始生成第一截口仍采用纯随机模型生成。一直到第k年。

2024/3/121254、实例（含如何检验）对长江支流乌江武隆站1954-1983年共30年洪水过程建立了季节性AR（1）模型，每年取最大20d平均流量作为其洪水过程（用日平均流量过程代表洪水过程，该流域平均）属大江大河。检验时，残差检验：独立性；截口特性检验：截口均值、Cv、Cs、最大值、最小值及一阶自相关系数，还有时段量特性检验，全部属于实用性检验。生成500个30年资料，对每个截口先计算检验参数数，再统计500组各检验参数的样本平均值与均方差2024/3/12126

经过计算绝大部分能通过一个方差检验（这个应该说较容易通过，因为建模时中考虑了这些截口特性参数）

实际上水文上更关心的是时段量，如最大洪量的平均值，的检验

2024/3/12127时段量长度19.989.800.835实测模拟0.3360.4250.067实测模拟9.989.800.835实测模拟750.8051.164.1530.3670.4120.0600.8820.7200.525326.5725.932.1780.3510.4270.06407020.9120.593从该表时段量检验结果中可以看出，大部分能通过一个均方差检验，但还是有偏差。事实上对于模拟洪水过程，还应打印出其过程线，看是否符合水文实际情况（凭肉眼看看，特别是看特大洪水过程）。

2024/3/12128是否相关独立除了不同截口间的独立性，还要看各截口的自相关如何，如

该模型参数较多，在样本长度n较小时，可能导致模型不稳定，误差大。如果考虑平稳模型则参数大大减少（作标准化变换后序列相关结构是平稳的）。季节性AR（1），还可应用于模拟月径流过程，缺点未直接考虑更高阶自相关，如等自相关关系，参数较多。

残差特性的独立性检验2024/3/12129二、解集模型（相关解集Disaggregationmodel） 形式：由ValenciaSchaake提出

2024/3/12130残差为：P-Ⅲ分布，均值0，方差1，（独立同分布）（X,Y都是离均系列）该模型把总量X分解成各个分量A反映不同分量所占比重，B表示不同分量之间相关程度。2024/3/12131当X表示年径流量，年径流等于不同月份月径流量之和。表示月径流量，当模拟洪水过程时，X表示每年洪水过程总量（如长江洪水模拟时，X表示180天洪水总量，而y表示每天流量值，显然）2024/3/12132参数估计，用矩法：已知资料设2024/3/12133写出通式：对对2024/3/12134偏态分布的考虑：2024/3/12135为了求残差偏态系数，此时应把分解成（水文序列非负定矩阵一般解满足）

随机模拟时，先生成X（可以用纯随机模型或AR（p）模拟）通过生成重复上述步骤可生成所需水文过程2024/3/12136优点：①结构简单，概念清楚；②总量分成分量后，各分量之间之和严格等于总量，即保持水量平衡；③应用较广。不足：①模型参数太多；②相关结构不一致。

实例对宜昌站180天洪水过程采用解集模型模拟效果较为理想（约100年资料），残差为

P-Ⅲ分布。

2024/3/12137第五节多变量模型实际上常会碰到模拟几个站点水文序列的情况，由于上下游之间存在相关关系，因此不能分别模拟，而应同时模拟。2024/3/12138多变量AR（1）与AR（2）模型：2024/3/12139不仅与有关，还与其他站时刻值有关。到底是什么分布？一般是正态分布。另一种对数正态变换后作变换生成，至于如何进行参数估计可自行查阅有关文献。2024/3/12140二阶模型：为正态分布。

先生成第一截口分布随机数，在通过建立第二截口与第一截口关系的生成公式：2024/3/12141长江三峡多站洪水模拟AR（3）及一种特殊变换:2024/3/12142第四章水文随机序列的预报概述平稳线性最少方差预报AR（p）序列预报门限自回归模型2024/3/12143第一节概述模机水文模型应用：最简单一个例子：对中心化（离均差后）系列建立了一阶自回归模型：（1）

—水文序列已知1984年2024/3/12144对（1）式取数学期望得：（预报公式）其实，这是一个期望预报，把随机变量取为0值（平均值）因此一步预报误差：已知水文序列作预报。①建立水文随机模型（假如无趋势及周期）②给出预报公式和预报方法。③利用现在及过去观测值作预测，必要时作实时修正。④区间预报，误差分析，对模型作评定及检验。2024/3/12145第二节:平稳线性最小方差预报对正态平稳系列作l步预报问题：当前时刻k和过去时刻为已知，需对时刻随机变量是，可记为。所谓平稳线性最少方差预报定义为：∴∴水文变量2024/3/12146一、ARMA（p、q）序列差分形式公式：离均化ARMA（p、q）两端取数学期望（条件,l>p,q)）当MA(0,q)模型时，如当为AR（p）模型时，将在后面作介绍。2024/3/12147二、ARMA（p、q）序列传递形式预报公式：把上式(2)式中t用代替2024/3/12148上式两边取条件数学期望该公式2024/3/12149三、ARMR（p、q）递推预报公式（可用于实时修正）对（2）式作一些变换，把k改为k+12024/3/12150事实上以上公式含义：要做k+1时刻l步预报等于在k时刻作l+1步预报值+（k+1时刻实测值与在k时刻对k+1时刻预测值之差）*权重系数。相当于对时刻作预测时，用误差来校正预测值，把这个称为实时校正（现时校正）。预测如显然在时刻预测比用预测效果要好，用了误差校正。2024/3/12151预报误差（评定与检验）[合格率（允许相对误差范围内,如20%）合格率>85%，甲等；70%-84%，乙等；60%-69%，丙等]（枯季径流允许误差30%，每日径流20%）2024/3/12152如区间估计。2024/3/12153第三节AR（p）序列预报（一般ARMA(p,q)模型有此结果）一步预报公式先估计模型参数2024/3/12154再作一步预报，又由预报如此递推可以作l步预测，当然还可以作实时校正预报。实例：某站中心化模型已建立（正态分布假定）求预测1983-1985平均流量。2024/3/12155①以1982年为k时刻进行预报区间预报2024/3/12156实时校正：以1983年为k+1进行实时校正预报2024/3/12157第四节:门限自回归模型(1978年提出H.Tong）（ThresholdAutoregressiveModel）

观察：①②

这样生成不是线性模型生成随机过程，若生成数据100年，用线性模型拟合效果较差即残差相关比较大，不一定独立，正态。2024/3/12158模型形式（一般）一般l取2～3，不宜取太多了，否则太复杂。这个模型实际说明了随机序列是分段线性的，即每个区间内可用自回归模型来描述，但这些模型序数在不同区间是不一样的。2024/3/12159建模：已知当用传统线性模型建模时发现，独立性差，或预测精度较差，则可以考虑用门限自回归模型来建模，（当然还应考虑是否在成因上做分段线性合理性）据以上模型采用以上已知数据，来辩认模型参数步骤如下：1：令。分成段，如确定各段自回归阶数，L为逐段自回归允许最大阶数，这里都认为各段有相同L值。2024/3/12160把实数轴分成l个区间一般是预报值取为这样可以把取值区间分l成个。对于已知观察分析动态数据中在不同区间个数（以排在开始统计），如L=3，d=1，从开始统计在不同区间个数。2024/3/12161记第一个区间N1个第二个区间N2个第l个区间Nl个设动态数据中有：对于第一段线性自回归模型具体形式已知：2024/3/12162用最小二乘法可求出可以想象，给定以上一组系数则可求出一组残差及其残差平方和，即可得第一段自回归模型，其他段求法一致。显然不同2024/3/121632、对于固定RSS(k1)为第1个模型残差平方和（最小二乘法估计参数）显然k1愈大，则残差愈小，但AIC考虑惩罚因子则第二项大，故找到。其他段的自回归模型阶数一样估计。2024/3/121643、固定

达到最少的即为门限值4、估计值AIC（）一般5、2024/3/12165不同那么d初值取多少？若效果不好，可先进行常用2024/3/12166第五章方差分析单因素方差分析多因素方差分析（双因素方差分析）2024/3/12167单因素方差分析（引例）某农科所为了比较四种不同的肥料对农作物产量的影响，进行了下面试验。他们选择采用一块肥沃程度比较均匀的土地，分成十六块。为了减少土地肥沃程度的影响，作了以上试验安排：A1表示第一种肥料A3表示第三种肥料A2表示第二种肥料A4表示第四种肥料2024/3/12168每种肥均施在由上到下不同层次的土地中，并可得出16块土地农作物产量通过上表，来推断不同肥料对农作物产量有无影响？如果影响显著，施哪一种肥料产量最高，而哪一种肥料最好？2024/3/12169在统计学中，考虑因素称为因子，因子不同状态称为水平，上例中一个因子是肥料，有四个不同水平，A1，A2，A3，A4

。由表又可知，肥料不同可能产生产量不同。因此，采用不同肥料所得产量不可以当作来自同一总体的样本，而应看作4个不同总体.抽取的样本长度为4的样本。通常假定这里要求假定方差相等，即方差齐性，那么要知道四种不同肥料对农作物产量的影响有关显著差别可归结为作检验：2024/3/1217012r总体参数结果因子一、数学模型①将作分解②求2024/3/12171其中反映了水平对试验结果的影响，称作水平的主效应，总结方差分析模型如下：2024/3/12172二、与估计三、统计检验误差平方和组内差不