CN110276104B - 一种集合气候模式下的分期设计洪水推求方法 - Google Patents

一种集合气候模式下的分期设计洪水推求方法 Download PDF

Info

Publication number
CN110276104B
CN110276104B CN201910433611.8A CN201910433611A CN110276104B CN 110276104 B CN110276104 B CN 110276104B CN 201910433611 A CN201910433611 A CN 201910433611A CN 110276104 B CN110276104 B CN 110276104B
Authority
CN
China
Prior art keywords
flood
climate
function
distribution
gcm
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201910433611.8A
Other languages
English (en)
Other versions
CN110276104A (zh
Inventor
尹家波
张官政
陈杰
顾磊
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Wuhan University WHU
Original Assignee
Wuhan University WHU
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Wuhan University WHU filed Critical Wuhan University WHU
Priority to CN201910433611.8A priority Critical patent/CN110276104B/zh
Publication of CN110276104A publication Critical patent/CN110276104A/zh
Application granted granted Critical
Publication of CN110276104B publication Critical patent/CN110276104B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/10—Geometric CAD
    • G06F30/13—Architectural design, e.g. computer-aided architectural design [CAAD] related to design of buildings, bridges, landscapes, production plants or roads
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A10/00—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE at coastal zones; at river basins
    • Y02A10/40—Controlling or monitoring, e.g. of flood or hurricane; Forecasting, e.g. risk assessment or mapping
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Geometry (AREA)
  • General Physics & Mathematics (AREA)
  • Computer Hardware Design (AREA)
  • Theoretical Computer Science (AREA)
  • Civil Engineering (AREA)
  • Structural Engineering (AREA)
  • Computational Mathematics (AREA)
  • Architecture (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Evolutionary Computation (AREA)
  • General Engineering & Computer Science (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开了一种集合气候模式下的分期设计洪水推求方法，通过建立GCM与VIC的耦合模型，以预测未来气候变化情景下的水库径流过程，通过算数平均法，获得气候变化情景下的平均态径流过程，采用概率变点法分析径流过程划分主汛期和非主汛期，通过Copula函数构建了水库设计洪水峰量最可能组合法计算模型，并将未来气候变化情景下的水库预测径流资料作为设计洪水峰量组合法计算模型的输入，同时考虑多变量同频率组合和最可能组合，计算水库在不同重现期水平下的联合设计值。本发明具有较强的统计基础，为气候变化情景下的分期设计洪水计算提供新途径。

Description

一种集合气候模式下的分期设计洪水推求方法
技术领域
本发明属于水库防洪安全设计技术领域，具体涉及一种集合气候模式下的分期设计洪水推求方法。
背景技术
分期设计洪水是确定分期汛限水位的重要依据，优化分期设计洪水，使之反映汛期洪水演变规律对准确刻画水库面临的防洪风险和提高利效益具有重要意义。分期设计洪水既要符合防洪标准，又要能反映洪水的季节性规律，我国现行的分期设计洪水方法假定分期设计洪水频率等于防洪标准的倒数，达不到防洪标准的要求。为了解决该问题，国内外学者主要提出了全概率公式法和联合分布法这两类技术方案。
全概率公式法假定各分期的年最大洪水相互独立，通过全概率公式建立年最大洪水与分期年最大洪水的关系式，但是该法在非主汛期抽取到的样本容量很少，最终导致该分期频率分析结果不可靠，难以在实际中采用。联合分布法的关键技术在于基于Copula函数的联合分布构建方法，能够充分挖掘主汛期与非主汛期洪水特征的内在相关性，例如：肖义等[1]通过Copula函数构建了分期设计洪水频率与防洪标准应满足的关系函数，推求了满足设计标准的分期洪水等值线；陈璐等[2]建立了基于Copula函数的洪水发生时间和量级的两变量联合分布，提出了一种新的分期设计洪水计算方法。
然而，以上基于Copula函数推求分期设计洪水的文献均得到重现期等值线，在等值线上有无数种洪水组合满足防洪标准，学者们为获得满足实际要求的分期设计洪水值，仍然采用同频率假设来推求分期设计洪水，并不符合分期洪水的内在发生规律。同时，受人类活动和气候变化的双重影响，下面条件发生了动态改变，从而使水库的入库径流规律发生变化，现有方法只考虑当前径流，未有文献对未来气候变化情景下的分期设计洪水展开研究。
参考文献：
[1]肖义，郭生练，刘攀，等.分期设计洪水频率与防洪标准关系研究[J].水科学进展，2008，19(1)：54-60.
[2]陈璐，郭生练，闫宝伟，等.一种新的分期设计洪水计算方法[J].武汉大学学报(工学版)，2010，43(1)：20-24.
发明内容
本发明的目的是克服现有技术存在的不足，基于全球集合气候模式(GCMs)，耦合统计降尺度方法和水文模型，获取未来气候变化情景下的降雨、径流过程，并基于Copula函数构建主汛期和非主汛期洪水的联合分布，同时考虑洪水的多变量同频率组合和最可能组合模式，提供一种集合气候模式下的分期设计洪水推求方法。本发明为解决现有技术中存在的问题采用的技术方案如下：
一种集合气候模式下的分期设计洪水推求方法，包括以下步骤：
步骤1，基于M个全球集合气候模式(GCMs)，通过基于经验分布的偏差校正方法对GCMs的降雨气温进行偏差校正，模拟气候变化情景下的降雨气温过程，并建立GCMs与可变下渗能力水文模型(VIC)的耦合模型，以模拟M组不同气候模式下的径流过程；
步骤2：对于步骤1的M组未来情景径流过程，通过算数平均法，获得气候变化情景下的平均态径流过程；
步骤3：对于步骤2获得的未来气候情景的平均态径流过程Qt，运用概率变点分析方法对汛期进行分期；
步骤4、依据步骤3的汛期分期结果，采用分期最大值取样方法，获得气候变化情景下的主汛期和非主汛期洪水，基于Copula函数建立未来情景的分期设计洪水计算模型，并推求气候变化情景下的多变量同频率组合和最可能分期设计洪水。
所述步骤1中由于GCMs的网格较大，不能满足流域或站点尺度的应用要求，往往需要通过统计降尺度方法建立大尺度气候模式输出变量与小尺度气象要素之间的统计关系，产生站点或流域尺度气候变化情景的过程，这是开展站点或流域尺度气候变化影响评估的重要技术环节。
所述步骤1中采用基于经验分布的偏差校正方法(Quantile-mapping baseddaily translation)作为统计降尺度方法来校正GCMs的降雨气温输出。
所述步骤1进一步包括以下子步骤：
步骤1.1、逐月计算历史时段模拟与实测日降水(包含无降水日)与气温序列经验分布函数：
Figure GDA0003562025860000031
Figure GDA0003562025860000032
其中，d表示月份；GCM表示全球气候模式(Global climate model)；ref表示基准期；ecdf表示累积分布函数；PGCM，ref，d、TGCM，ref，d分别为历史时段各月(1、2、3、...12月)降水、气温经验累积分布函数；
Figure GDA0003562025860000033
和
Figure GDA0003562025860000034
分别表示GCM模式下的降水、气温系列。
步骤1.2、逐月计算经验分布各概率下日降水(观测序列与模拟序列比值)与日气温(观测序列与模拟序列差值)校正因子：
Figure GDA0003562025860000041
Figure GDA0003562025860000042
其中，CFPGCM，ref，d、CFTGCM，ref，d分别为历史时段各月降水、气温校正因子；
步骤1.3、逐月内插并外延根据经验分布计算的日降水与日气温校正因子：
Figure GDA0003562025860000043
Figure GDA0003562025860000044
其中，PGCM，fut，d和TGCM，fut，d分别表示GCM模式下未来情景的降水、气温系列；CFPGCM，fut，d、CFTGCM，fut，d分别为处理后应用于未来时段的降水、气温校正因子；
步骤1.4、将处理后的校正因子应用于未来模拟日降水与日气温：
Figure GDA0003562025860000045
Figure GDA0003562025860000046
其中，
Figure GDA0003562025860000047
分别为降尺度后的未来降水、气温情景；
步骤1.5、采用基于经验分布的偏差校正方法(Quantile-mapping based dailytranslation)作为统计降尺度方法来校正GCMs的输出(降水和气温)作为流域水文模型VIC的输入，模拟预测未来气候情景下的流域径流过程；
进一步地，由于单一GCM模型具有较大的不确定性，本文采用M个基于分位数的偏差校正方法得出的全球气候模式输出(GCM)与可变下渗能力水文(VIC)模型耦合，于是在给定的气候情景下，可以模拟M组降雨、径流系列。
所述步骤2中对于第i个GCM模型，通过耦合统计降尺度和水文模型，可以在步骤1中得到气候变化情景下的径流过程
Figure GDA0003562025860000048
其中t代表时间变量，则可通过算数平均法获得未来情景的平均态径流过程：
Figure GDA0003562025860000051
式中，Qt为气候变化情景下气候模式的平均态径流过程。
所述步骤3进一步包括以下子步骤：
步骤3.1、确定一个阈值
Figure GDA0003562025860000052
即认为
Figure GDA0003562025860000053
是洪灾风险发生与否的临界值点，当
Figure GDA0003562025860000054
时，认为洪水发生；
Figure GDA0003562025860000055
时，认为洪水不发生；阈值选取，根据不同流域的水情特征来选定，例如可选取洪水年最大值系列的50％分位数作为阈值。
步骤3.2、改变分期个数，分别计算最优的变点估计，最后，取式(9)的值开始缓慢变化的两个点作为主汛期和非主汛期洪水分期点；
所述步骤3.2进一步包括如下步骤：
步骤3.2.1、设n1，n2分别代表洪水发生和不发生的次数，则统计量
Figure GDA0003562025860000056
的条件分布服从样本大小分别为n1，n2的斯米尔诺夫分布相同，有检验方法如下：当n1，n2小于40时，对于给定的置信度水平α，从斯米尔诺夫表中可以查到对应于n1，n2临界值Cα，当：
Figure GDA0003562025860000057
成立时，则认为变点存在；
步骤3.2.2、当n1和n2大于40时，如果：
Figure GDA0003562025860000058
成立，则认为变点存在，否则变点不存在，常用的两个值为：C0.005＝1.358，C0.01＝1.628，同样，这里也可以用秩和检验法或游程检验法来进行假设检验；
为了实现汛期时间序列的多变点分析，同样可以仿照均值变点的办法，先给定变点个数以及变点的初始估计，固定两端的点，用累次计数法估计中间的变点，直到估计的变点位置不再发生变化为止。同样，最后得出的计算结果与初始点的选取有关；
设总共有n年降水模型，每个时段事件发生次数xi(i＝1，2，...，l)次，则xi服从二项分布，由此得出变点mj(j＝1，2，...，q)的似然函数：
Figure GDA0003562025860000061
式中，pi表示在分期i事件发生的概率，可用下式估计：
Figure GDA0003562025860000062
利用不同初始值求出来的估计变点，取似然函数值最大的作为最终的分期方式，在实际计算中，为了避免直接计算式(前12式)引起的数值溢出，注意到组合项都相同，同时实际操作中只需要比较大小，因此可仅通过下式计算似然函数：
Figure GDA0003562025860000063
所述步骤4进一步包括以下子步骤：
步骤4.1、基于步骤3的分期结果，采用分期最大值取样方法获得步骤2中得到的未来气候情景下的主汛期和非主汛期洪水样本；
现有的洪水取样方法主要有年最大洪水取样、超定量取样和分期最大值取样等三种方法，其中分期最大值取样方法在分期设计洪水计算中应用较广，属于现有技术；
步骤4.2、构造边缘分布和联合分布函数，令X表示气候变化情景的主汛期洪水；Yi(i＝1，2，...，n)表示气候变化情景的非主汛期洪水，采用P-III(皮尔逊III型)分布函数构建X、Y的边缘分布，分别记为FX(x)，
Figure GDA0003562025860000071
密度函数分别记为fX(x)，
Figure GDA0003562025860000072
采用线性矩法估计P-III分布函数的参数；
采用Copula函数构造气候变化情景的分期洪水的联合分布，令Q(x1，x2，...，xn)为一个n维分布函数，其边缘分布分别为F1(x1)，F2(x2)，...，Fn(xn)，则存在一个n-Copula函数C，使得对任意x∈Rn(x为n维向量，Rn为n维实数空间)：
Q(x1，x2，...，xn)＝Cθ(F1(x1)，F2(x2)，...，Fn(xn)) (15)
式中：θ为Copula函数的相关性参数；
优选地，本发明以Archimedean Copula函数族中的Gumbel-Hougaard(G-H)Copula函数作为联合分布函数(参考式15)，构造主汛期和非主汛期洪水的联合分布；
步骤4.3、基于OR重现期作为防洪标准，基于多变量同频率组合构建主汛期与非主汛期的同频率设计值计算模型；
采用OR重现期作为水库防洪标准的度量指标，其定义如下：
T(x，y1，...，yn)＝1/[1-F(x，y1，...，yn)] (16)
式中：T(x，y1，...，yn)即为联合重现期，以年为单位；F(x，y1，...，yn)为联合分布函数；对于主汛期与非主汛期的同频率组合，应满足下式：
Figure GDA0003562025860000073
式中，T为重现期，
Figure GDA0003562025860000074
为采用G-H copula函数构造的联合分布函数。
Copula函数均具有显式表达式，由于设定式(17)中不同变量的累积概率值相等，故可通过解析方法求解得到同频率设计值；
例如，二维G-H copula函数表达式为：
C(u，v1)＝exp{-[(-ln u)θ+(-ln v1)θ]1/θ} (18)
结合式(17)-(18)，可以求得X和Y1的累积概率分别为：
Figure GDA0003562025860000089
式中，
Figure GDA0003562025860000081
θ为Copula函数的参数；
进而，可通过边缘分布求得同频率设计值组合：
Figure GDA0003562025860000082
Figure GDA0003562025860000083
式中：
Figure GDA0003562025860000084
分别为X和Y1的反函数；
类似地，三维对称型G-H copula函数表达式为：
C(u，v1，v2)＝exp{-[(-ln u)θ+(-ln v1)θ+(-ln v2)θ]1/θ} (22)
式中，
Figure GDA0003562025860000085
θ为Copula函数的参数；
Figure GDA0003562025860000086
则可以求解出三变量的同频率组合模式为：
Figure GDA00035620258600000811
同理，将式(23)得到的累积概率值(u，v1，v2)代入到边缘分布(FX(x)、
Figure GDA00035620258600000810
即可求解；
步骤4.4、基于OR重现期作为防洪标准，构建主汛期与非主汛期的最可能设计值计算模型，并推求分期设计洪水的最可能设计值；
F(x，y1，...，yn)相应的联合概率密度函数可以表示为：
Figure GDA0003562025860000087
式中：
Figure GDA0003562025860000088
为Copula函数的密度函数；fX、
Figure GDA00035620258600000812
分别为X和Yi(i＝1，2，...，n)的概率密度函数；
步骤4.5、求解联合概率密度函数最大的汛期洪水组合，推求气候变化情景的最可能分期设计洪水；
分期设计洪水最可能组合模式是指(X，Y)在满足防洪标准的条件下，f(x，y1，...，yn)取最大值时的两变量联合设计值，通过以下方程求解该问题：
Figure GDA0003562025860000091
函数f(x，y1，...，yn)在定义域内是连续的，故必然存在最大值和最小值，本文提出采用拉格朗日乘数法求解该问题，给定联合重现期TU，构造拉格朗日函数如下：
Figure GDA0003562025860000092
分别对x，y1，...，yn和λ求偏导，并令其为0，得：
Figure GDA0003562025860000093
式中：λ为拉格朗日乘子；c＝c(u，v1，v2，...，vn-1，vn)，
Figure GDA0003562025860000094
Figure GDA0003562025860000095
pT＝1-1/T；fX′，
Figure GDA0003562025860000096
分别为相应概率密度函数的导函数。例如，当构建的主汛期与非主汛期为二维联合分布时，可通过求解下式求得分期设计洪水的最可能组合值，即对应公式(28)中求解出来的x和y1：
Figure GDA0003562025860000101
本发明具有如下优点：
1、本发明耦合全球气候集合模式、基于经验分布的偏差校正方法和水文模型，模拟了气候变化情景下的气象水文过程，充分考虑气候变化对径流季节性特征的影响，采用概率变点分析方法对汛期进行分期，为未来水资源综合管理提供重要的应用支撑，科学合理、贴近工程实际；
2、本发明基于气候变化情景的径流系列，采用分期最大值取样方法获得主汛期和非主汛期的洪水样本，通过多维Copula函数建立气候变化情景下主汛期和非主汛期洪水的联合分布，充分考虑洪水特征的内在相关性，采用多变量同频率组合和最可能组合模式推求了分期洪水设计值，为工程应用提供了重要且可操作性强的参考依据。
本发明基于M个全球气候模式(GCMs)，通过基于经验分布的偏差校正方法对GCMs的降雨气温进行偏差校正，模拟气候变化情景下的降雨气温过程，并建立GCMs与可变下渗能力水文模型(VIC)的耦合模型，模拟未来气候变化情景下的水库径流过程，通过分析气候变化情景下的径流特征开展汛期分期，并通过Copula函数构建基于多变量组合和最可能组合的分期设计洪水计算模型，以适应未来气候变化，更具工程实用性。
附图说明
图1为本发明方法的具体流程图；
图2为主汛期洪水的边缘分布的示意图；
图3为主汛期与非主汛期洪水二维联合分布的示意图。
具体实施方式
下面通过实施例，并结合附图，对本发明的技术方案作进一步具体的说明，一种集合气候模式下的分期设计洪水推求方法，包括以下步骤：
步骤1，基于M个全球集合气候模式(GCMs)，通过基于经验分布的偏差校正方法对GCMs的降雨气温进行偏差校正，模拟气候变化情景下的降雨气温过程，并建立GCMs与可变下渗能力水文模型(VIC)的耦合模型，以模拟M组不同气候模式下的径流过程；
由于GCMs的网格较大，不能满足流域或站点尺度的应用要求，往往需要通过统计降尺度方法建立大尺度气候模式输出变量与小尺度气象要素之间的统计关系，产生站点或流域尺度气候变化情景的过程，这是开展站点或流域尺度气候变化影响评估的重要技术环节；
进一步地，本发明采用基于经验分布的偏差校正方法(Quantile-mapping baseddaily translation)作为统计降尺度方法来校正GCMs的降雨气温输出。
步骤1进一步包括以下子步骤：
步骤1.1、逐月计算历史时段模拟与实测日降水(包含无降水日)与气温序列经验分布函数：
Figure GDA0003562025860000111
Figure GDA0003562025860000112
其中，PGCM，ref，d、TGCM，ref，d分别为历史时段各月(1、2、3、...12月)降水、气温经验累积分布函数；
步骤1.2、逐月计算经验分布各概率下日降水(观测序列与模拟序列比值)与日气温(观测序列与模拟序列差值)校正因子：
Figure GDA0003562025860000121
Figure GDA0003562025860000122
其中，CFPGCM，ref，d、CFTGCM，ref，d分别为历史时段各月降水、气温校正因子；
步骤1.3、逐月内插并外延根据经验分布计算的日降水与日气温校正因子：
Figure GDA0003562025860000123
Figure GDA0003562025860000124
其中，CFPGCM，fut，d、CFTGCM，fut，d分别为处理后应用于未来时段的降水、气温校正因子；
步骤1.4、将处理后的校正因子应用于未来模拟日降水与日气温：
Figure GDA0003562025860000125
Figure GDA0003562025860000126
其中，
Figure GDA0003562025860000127
分别为降尺度后的未来降水、气温情景；
步骤1.5、采用基于分位数的日偏差校正方法(DBC)作为统计降尺度方法来校正GCMs的输出(降水和气温)作为流域水文模型VIC的输入，模拟预测未来气候情景下的流域径流过程。
进一步地，由于单一GCM模型具有较大的不确定性，本文采用M个基于分位数的偏差校正方法得出的全球气候模式输出(GCM)与可变下渗能力水文(VIC)模型耦合，于是在给定的气候情景下，可以模拟M组降雨、径流系列。
步骤2：对于步骤1的M组未来情景径流过程，通过算数平均法，获得气候变化情景下的平均态径流过程；
对于第i个GCM模型，通过耦合统计降尺度和水文模型，可以在步骤1中得到气候变化情景下的径流过程
Figure GDA0003562025860000131
其中t代表时间变量，则可通过算数平均法获得未来情景的平均态径流过程：
Figure GDA0003562025860000132
式中，Qt为气候变化情景下气候模式的平均态径流过程。
步骤3：对于步骤2获得的未来气候情景的平均态径流过程Qt，运用概率变点分析方法对汛期进行分期；
步骤3进一步包括以下子步骤：
步骤3.1、确定一个阈值
Figure GDA0003562025860000133
即认为
Figure GDA0003562025860000134
是洪灾风险发生与否的临界值点，当
Figure GDA0003562025860000135
时，认为洪水发生；
Figure GDA0003562025860000136
时，认为洪水不发生；
步骤3.2、对计算出来的分期进行假设检验；改变分期个数，分别计算最优的变点估计，最后，取式(9)的值开始缓慢变化的两个点作为主汛期和非主汛期洪水分期点；
步骤3.2.1、设n1，n2分别代表洪水发生和不发生的次数，则统计量
Figure GDA0003562025860000137
的条件分布服从样本大小分别为n1，n2的斯米尔诺夫分布相同，有检验方法如下：当n1，n2小于40时，对于给定的置信度水平d，从斯米尔诺夫表中可以查到对应于n1，n2临界值Cα：
Figure GDA0003562025860000138
成立时，则认为变点存在；
步骤3.2.2、当n1和n2大于40时，如果
Figure GDA0003562025860000139
成立，则认为变点存在，否则变点不存在，常用的两个值为：C0.005＝1.358，C0.01＝1.628，同样，这里也可以用秩和检验法或游程检验法来进行假设检验；
累次计数法只能分析单个变点，为了实现汛期时间序列的多变点分析，同样可以仿照均值变点的办法，先给定变点个数以及变点的初始估计，固定两端的点，用累次计数法估计中间的变点，直到估计的变点位置不再发生变化为止。同样，最后得出的计算结果与初始点的选取有关；
设总共有n年降水模型，每个时段事件发生次数xi(i＝1，2，...，l)次，则xi服从二项分布，由此得出变点mj(j＝1，2，...，q)的似然函数：
Figure GDA0003562025860000141
式中，pi表示在分期i事件发生的概率，可用下式估计：
Figure GDA0003562025860000142
利用不同初始值求出来的估计变点，取似然函数值最大的作为最终的分期方式，在实际计算中，为了避免直接计算式(前12式)引起的数值溢出，注意到组合项都相同，同时实际操作中只需要比较大小，因此可仅比较下式的大小
Figure GDA0003562025860000143
步骤4、依据步骤3的汛期分期结果，采用分期最大值取样方法，获得气候变化情景下的主汛期和非主汛期洪水，基于Copula函数建立未来情景的分期设计洪水计算模型，并推求气候变化情景下的多变量同频率组合和最可能分期设计洪水；
步骤4进一步包括以下子步骤：
步骤4.1、基于步骤3的分期结果，采用分期最大值取样方法获得步骤2中得到的未来气候情景下的主汛期和非主汛期洪水样本；
现有的洪水取样方法主要有年最大洪水取样、超定量取样和分期最大值取样等三种方法，其中分期最大值取样方法在分期设计洪水计算中应用较广，属于现有技术；
步骤4.2、构造边缘分布和联合分布函数，令X表示气候变化情景的主汛期洪水；Yi(i＝1，2，...，n)表示气候变化情景的非主汛期洪水，采用P-III分布函数构建X、Y的边缘分布，分别记为FX(x)，
Figure GDA0003562025860000151
密度函数分别记为fX(x)，
Figure GDA0003562025860000152
采用线性矩法估计P-III分布函数的参数。
如图3所示，给出了主汛期洪水的边缘分布的示意图。
采用Copula函数构造气候变化情景的分期洪水的联合分布，令Q(x1，x2，...，xn)为一个n维分布函数，其边缘分布分别为F1(x1)，F2(x2)，...，Fn(xn)，则存在一个n-Copula函数C，使得对任意x∈Rn(x为n维向量，Rn为n维实数空间)：
Q(x1，x2，...，xn)＝Cθ(F1(x1)，F2(x2)，...，Fn(xn)) (15)
式中：θ为Copula函数的相关性参数；
优选地，本发明以Archimedean Copula函数族中的Gumbel-Hougaard(G-H)Copula函数作为联合分布函数，构造主汛期和非主汛期洪水的联合分布；
步骤4.3、基于OR重现期作为防洪标准，基于多变量同频率组合构建主汛期与非主汛期的同频率设计值计算模型；
采用AND重现期作为水库防洪标准的度量指标，其定义如下：
T(x，y1，...，yn)＝1/[1-F(x，y1，...，yn)] (16)
式中：T(x，y1，...，yn)即为联合重现期，以年为单位；F(x，y1，...，yn)为联合分布函数；对于主汛期与非主汛期的同频率组合，应满足下式：
Figure GDA0003562025860000161
Copula函数均具有显式表达式，由于设定式(17)中不同变量的累积概率值相等，故可通过解析方法求解得到同频率设计值；
例如，二维G-H copula函数表达式为：
C(u，v1)＝exp{-[(-ln u)θ+(-ln v1)θ]1/θ} (18)
结合式(17)-(18)，可以求得X和Y1的累积概率分别为：
Figure GDA0003562025860000168
式中，
Figure GDA0003562025860000162
θ为Copula函数的参数；
如图3所示，给出了主汛期和非主汛期二维联合分布的等值线和同频率设计值；
进而，可通过边缘分布求得同频率设计值组合：
Figure GDA0003562025860000163
Figure GDA0003562025860000164
式中：
Figure GDA0003562025860000165
分别为X和y1的反函数；
类似地，三维对称型G-H copula函数表达式为：
C(u，v1，v2)＝exp{-[(-ln u)θ+(-ln v1)θ+(-ln v2)θ]1/θ} (22)
式中，
Figure GDA0003562025860000166
则可以求解出三变量的同频率组合模式为：
Figure GDA0003562025860000167
同理，将式(14)得到的累积概率值代入到边缘分布，即可求解；
步骤4.4、基于OR重现期作为防洪标准，构建主汛期与非主汛期的最可能设计值计算模型，并推求分期设计洪水的最可能设计值；
F(x，y1，...，yn)相应的联合概率密度函数可以表示为：
Figure GDA0003562025860000171
式中：
Figure GDA0003562025860000172
为Copula函数的密度函数；fX、
Figure GDA0003562025860000179
分别为X和Yi(i＝1，2，...，n)的概率密度函数；
步骤4.5、求解联合概率密度函数最大的汛期洪水组合，推求气候变化情景的最可能分期设计洪水；
分期设计洪水最可能组合模式是指(X，Y)在满足防洪标准的条件下，f(x，y1，...，yn)取最大值时的两变量联合设计值。通过以下方程求解该问题：
Figure GDA0003562025860000173
函数f(x，y1，...，yn)在定义域内是连续的，故必然存在最大值和最小值。本文提出采用拉格朗日乘数法求解该问题，给定联合重现期T∪，构造拉格朗日函数如下：
Figure GDA0003562025860000174
分别对x，y1，...，yn和λ求偏导，并令其为0，得：
Figure GDA0003562025860000175
式中：λ为拉格朗日乘子；c＝c(u，v1，v2，...，vn-1，vn)，
Figure GDA0003562025860000176
Figure GDA0003562025860000177
pT＝1-1/T；fX′，
Figure GDA0003562025860000178
分别为相应概率密度函数的导函数；例如，当构建的主汛期与非主汛期为二维联合分布时，可通过求解下式求得分期设计洪水的最可能组合值：
Figure GDA0003562025860000181
本发明的保护范围并不限于上述的实施例，显然，本领域的技术人员可以对本发明进行各种改动和变形而不脱离本发明的范围和精神。倘若这些改动和变形属于本发明权利要求及其等同技术的范围内，则本发明的意图也包含这些改动和变形在内。

Claims (4)

1.一种集合气候模式下的分期设计洪水推求方法，其特征在于，包括以下步骤：
步骤1、基于M个全球集合气候模式(GCMs)，通过基于经验分布的偏差校正方法对GCMs的降雨气温进行偏差校正，模拟气候变化情景下的降雨气温过程，并建立GCMs与可变下渗能力水文模型(VIC)的耦合模型，以模拟M组不同气候模式下的径流过程；
步骤2、对于步骤1的M组未来情景径流过程，通过算数平均法，获得气候变化情景下的平均态径流过程；
步骤3、对于步骤2获得的未来气候情景的平均态径流过程Qt，运用概率变点分析方法对汛期进行分期；
步骤4、依据步骤3的汛期分期结果，采用分期最大值取样方法，获得气候变化情景下的主汛期和非主汛期洪水，基于Copula函数建立未来情景的分期设计洪水计算模型，并推求气候变化情景下的多变量同频率组合和最可能分期设计洪水；
所述步骤1中采用基于经验分布的偏差校正方法作为统计降尺度方法来校正GCMs的降雨气温输出；
所述步骤3进一步包括以下子步骤：
步骤3.1、确定一个阈值
Figure FDA0003594293360000011
即认为
Figure FDA0003594293360000012
是洪灾风险发生与否的临界值点，当
Figure FDA0003594293360000013
时，认为洪水发生；
Figure FDA0003594293360000014
时，认为洪水不发生；其中阈值选取，根据不同流域的水情特征来选定：
步骤3.2、改变分期个数，分别计算最优的变点估计，最后，取式(9)的值开始缓慢变化的两个点作为主汛期和非主汛期洪水分期点；
所述步骤4进一步包括以下子步骤：
步骤4.1、基于步骤3的分期结果，采用分期最大值取样方法获得步骤2中得到的未来气候情景下的主汛期和非主汛期洪水样本；
步骤4.2、构造边缘分布和联合分布函数，令X表示气候变化情景的主汛期洪水；Yi表示气候变化情景的非主汛期洪水，i＝1,2,...,n，采用P-III(皮尔逊III型)分布函数构建X、Y的边缘分布，分别记为FX(x)，
Figure FDA0003594293360000021
密度函数分别记为fX(x)，
Figure FDA0003594293360000022
采用线性矩法估计P-III分布函数的参数；
采用Copula函数构造气候变化情景的分期洪水的联合分布，令Q(x1,x2,...,xn)为一个n维分布函数，其边缘分布分别为F1(x1),F2(x2),...,Fn(xn)，则存在一个n-Copula函数C，使得对任意x∈Rn，x为n维向量，Rn为n维实数空间：
Q(x1,x2,...,xn)＝Cθ(F1(x1),F2(x2),...,Fn(xn)) (15)
式中：θ为Copula函数的相关性参数；
步骤4.3、基于OR重现期作为防洪标准，基于多变量同频率组合构建主汛期与非主汛期的同频率设计值计算模型；
采用OR重现期作为水库防洪标准的度量指标,其定义如下：
T(x,y1,...,yn)＝1/[1-F(x,y1,...,yn)] (16)
式中：T(x,y1,...,yn)即为联合重现期，以年为单位；F(x,y1,...,yn)为联合分布函数；对于主汛期与非主汛期的同频率组合，应满足下式：
Figure FDA0003594293360000023
式中，T为重现期，
Figure FDA0003594293360000024
为采用G-H copula函数构造的联合分布函数；
步骤4.4、基于OR重现期作为防洪标准，构建主汛期与非主汛期的最可能设计值计算模型，并推求分期设计洪水的最可能设计值；
F(x,y1,...,yn)相应的联合概率密度函数可以表示为：
Figure FDA0003594293360000031
式中：
Figure FDA0003594293360000032
为Copula函数的密度函数；fX、
Figure FDA0003594293360000033
分别为X和Yi的概率密度函数，其中i＝1,2,…,n，u＝FX(x)、
Figure FDA0003594293360000034
Figure FDA0003594293360000035
步骤4.5、求解联合概率密度函数最大的汛期洪水组合，推求气候变化情景的最可能分期设计洪水；
分期设计洪水最可能组合模式是指(X，Y)在满足防洪标准的条件下，f(x,y1,…,yn)取最大值时的两变量联合设计值，通过以下方程求解该两变量联合设计值：
Figure FDA0003594293360000036
函数f(x,y1,…,yn)在定义域内是连续的，故必然存在最大值和最小值，本文提出采用拉格朗日乘数法求解分期设计洪水最可能组合值，给定联合重现期TU，构造拉格朗日函数如下：
Figure FDA0003594293360000037
分别对x,y1,…,yn和λ求偏导，并令其为0，得：
Figure FDA0003594293360000041
式中：λ为拉格朗日乘子；c＝c(u，v1，v2，...，vn-1，vn)，
Figure FDA0003594293360000042
Figure FDA0003594293360000043
pT＝1-1/T；fX′，
Figure FDA0003594293360000044
分别为相应概率密度函数的导函数，其中i＝1,2,…,n。
2.如权利要求1所述的一种集合气候模式下的分期设计洪水推求方法，其特征在于，所述步骤1进一步包括以下子步骤：
步骤1.1、逐月计算历史时段模拟与实测日降水与气温序列经验分布函数，日降水包含无降水日：
Figure FDA0003594293360000045
Figure FDA0003594293360000046
其中，d表示月份；GCM表示全球气候模式(Global climate model)；ref表示基准期；ecdf表示累积分布函数；PGCM,ref,d、TGCM,ref,d分别为历史时段各月降水，各月包含1、2、3、…、12月，和气温经验累积分布函数；
Figure FDA0003594293360000047
和
Figure FDA0003594293360000048
分别表示GCM模式下的降水、气温系列；
步骤1.2、逐月计算经验分布各概率下日降水：观测序列与模拟序列比值，与日气温：观测序列与模拟序列差值，校正因子：
Figure FDA0003594293360000049
Figure FDA0003594293360000051
其中，CFPGCM,ref,d、CFTGCM,ref,d分别为历史时段各月降水、气温校正因子；
步骤1.3、逐月内插并外延根据经验分布计算的日降水与日气温校正因子：
Figure FDA0003594293360000052
Figure FDA0003594293360000053
其中，PGCM，fut，d和TGCM，fut，d分别表示GCM模式下未来情景的降水、气温系列；CFPGCM,fut,d、CFTGCM,fut,d分别为处理后应用于未来时段的降水、气温校正因子；
步骤1.4、将处理后的校正因子应用于未来模拟日降水与日气温：
Figure FDA0003594293360000054
Figure FDA0003594293360000055
其中，
Figure FDA0003594293360000056
分别为降尺度后的未来降水、气温情景；
步骤1.5、采用基于经验分布的偏差校正方法(Quantile-mapping based dailytranslation)作为统计降尺度方法来校正GCMs的输出：降水和气温，作为流域水文模型VIC的输入，模拟预测未来气候情景下的流域径流过程。
3.如权利要求1所述的一种集合气候模式下的分期设计洪水推求方法，其特征在于：所述步骤2中对于第i个GCM模型，通过耦合统计降尺度和水文模型，可以在步骤1中得到气候变化情景下的径流过程
Figure FDA0003594293360000057
其中t代表时间变量，则可通过算数平均法获得未来情景的平均态径流过程：
Figure FDA0003594293360000058
式中，Qt为气候变化情景下气候模式的平均态径流过程。
4.如权利要求1所述的一种集合气候模式下的分期设计洪水推求方法，其特征在于，所述步骤3.2进一步包括如下步骤：
步骤3.2.1、设n1，n2分别代表洪水发生和不发生的次数，则统计量
Figure FDA0003594293360000061
的条件分布服从样本大小分别为n1，n2的斯米尔诺夫分布相同，有检验方法如下：当n1，n2小于40时，对于给定的置信度水平α，从斯米尔诺夫表中可以查到对应于n1，n2临界值Ca，当：
Figure FDA0003594293360000062
成立时，则认为变点存在；
步骤3.2.2、当n1和n2大于40时，如果：
Figure FDA0003594293360000063
成立，则认为变点存在，否则变点不存在，常用的两个值为：C0.005＝1.358，C0.01＝1.628，同样，此处可用秩和检验法或游程检验法来进行假设检验；
为实现汛期时间序列的多变点分析，同样可仿照均值变点的办法，先给定变点个数以及变点的初始估计，固定两端的点，用累次计数法估计中间的变点，直到估计的变点位置不再发生变化为止，最后得出的计算结果与初始点的选取有关；
设总共有n年降水模型，每个时段事件发生次数xi次，i＝1,2,…,l，则xi服从二项分布，由此得出变点mj的似然函数，j＝1,2,…,q：
Figure FDA0003594293360000064
式中，pi表示在分期i事件发生的概率，可用下式估计：
Figure FDA0003594293360000065
利用不同初始值求出来的估计变点，取似然函数值最大的作为最终的分期方式，在实际计算中，为了避免直接计算式(1)-(12)引起的数值溢出，注意到组合项都相同，同时实际操作中只需要比较大小，因此可仅通过下式计算似然函数：
Figure FDA0003594293360000071
CN201910433611.8A 2019-05-23 2019-05-23 一种集合气候模式下的分期设计洪水推求方法 Active CN110276104B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910433611.8A CN110276104B (zh) 2019-05-23 2019-05-23 一种集合气候模式下的分期设计洪水推求方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910433611.8A CN110276104B (zh) 2019-05-23 2019-05-23 一种集合气候模式下的分期设计洪水推求方法

Publications (2)

Publication Number Publication Date
CN110276104A CN110276104A (zh) 2019-09-24
CN110276104B true CN110276104B (zh) 2022-06-07

Family

ID=67960205

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910433611.8A Active CN110276104B (zh) 2019-05-23 2019-05-23 一种集合气候模式下的分期设计洪水推求方法

Country Status (1)

Country Link
CN (1) CN110276104B (zh)

Families Citing this family (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111027175B (zh) * 2019-11-06 2021-09-24 中国地质大学（武汉） 基于耦合模型集成模拟的洪水对社会经济影响的评估方法
CN111611692B (zh) * 2020-04-26 2022-08-05 武汉大学 气候变化情景下基于等可靠度的设计洪水推求方法及系统
CN111797129B (zh) * 2020-06-01 2024-01-26 武汉大学 一种气候变化情景下水文旱情评估方法
CN112036683B (zh) * 2020-07-14 2023-02-14 中国电建集团华东勘测设计研究院有限公司 一种适用于未来气候变化情景下的水库防洪风险预估方法
CN112819293B (zh) * 2021-01-14 2023-01-06 中国长江三峡集团有限公司 气候变化影响下水库调度规则的失效预警分析方法
CN112765912B (zh) * 2021-01-26 2022-06-14 武汉大学 基于气候模式集合的洪涝灾害社会经济暴露度的评估方法
CN113158542B (zh) * 2021-01-29 2022-10-04 武汉大学 适用于缺资料地区的多变量设计洪水估计方法
CN117522012A (zh) * 2023-11-02 2024-02-06 长江水利委员会水文局 一种基于季节周期特性的径流场景生成方法

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
GB201405788D0 (en) * 2014-03-31 2014-05-14 Imp Innovations Ltd A computer implemented method of deriving performance from a local model
CN105714729B (zh) * 2016-02-29 2017-10-24 武汉大学 一种自适应气候变化的水库多变量设计洪水推求方法

Also Published As

Publication number Publication date
CN110276104A (zh) 2019-09-24

Similar Documents

Publication Publication Date Title
CN110276104B (zh) 一种集合气候模式下的分期设计洪水推求方法
Yan et al. Comparison of four nonstationary hydrologic design methods for changing environment
Kou et al. Sparse online warped Gaussian process for wind power probabilistic forecasting
Mendoza et al. An intercomparison of approaches for improving operational seasonal streamflow forecasts
Stauffer et al. Ensemble postprocessing of daily precipitation sums over complex terrain using censored high-resolution standardized anomalies
CN110598290B (zh) 考虑气候变化的流域未来水电发电能力预测方法和系统
CN111611692B (zh) 气候变化情景下基于等可靠度的设计洪水推求方法及系统
CN108898250B (zh) 一种基于D藤copula函数的月径流模拟方法
CN107563554B (zh) 一种统计降尺度模型预报因子的筛选方法
Smithers et al. Design rainfall and flood estimation in South Africa
CN102495937A (zh) 一种基于时间序列的预测方法
Hostache et al. Propagation of uncertainties in coupled hydro-meteorological forecasting systems: A stochastic approach for the assessment of the total predictive uncertainty
Chen et al. Copula-based method for stochastic daily streamflow simulation considering lag-2 autocorrelation
Hendrickx et al. Impact of warming climate on water management for the Ariège River basin (France)
Lemaitre-Basset et al. Climate change impact and uncertainty analysis on hydrological extremes in a French Mediterranean catchment
Nasri et al. Non-stationary hydrologic frequency analysis using B-spline quantile regression
CN115238947A (zh) 气候变化下旱涝急转事件的社会经济暴露度预估方法
Zhang et al. A weighted ensemble of regional climate projections for exploring the spatiotemporal evolution of multidimensional drought risks in a changing climate
CN109460526A (zh) 基于Copula函数和信息熵理论相结合的综合式水文站网评估模型
CN109165455A (zh) 基于互信息和vine copula的水文相依结构建模方法
CN114492952B (zh) 一种基于深度学习的短临降水预报方法和装置
Oesting et al. Sampling from a Max‐Stable Process Conditional on a Homogeneous Functional with an Application for Downscaling Climate Data
Wilson et al. Use of meteorological data for improved estimation of risk in capacity adequacy studies
CN109614642A (zh) 一种水文频率分析的模型选择方法
CN106557614B (zh) 一种基于Halphen IB分布的洪水频率分析方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant