基于Copula函数的黑龙江省作物生育期丰枯遭遇研究

李红芳，王会肖，赵茹欣，杨雅雪，郭嘉豪

（北京师范大学水科学研究院，北京100875）

0 引 言

随着粮食产量的减少、水资源的供给不足以及气候变化等问题的出现，分析研究区域农业干旱特征，对于加强对干旱的预防和应急管理，促进区域水资源的合理配置具有重要意义。降雨量和参考作物蒸腾量是灌溉系统中两个具有相关关系的重要水文随机变量，可表征灌区供需水变化关系[1-3]。因此十分有必要通过研究降雨量和参考作物腾发量的二维联合概率分布模型，来揭示自然降水条件下的灌区水资源短缺风险[4]，近年来，多变量联合分析成为水文计算领域的一个研究热点，并证实比单变量分析能更好地描述水文事件的内在规律和分析各个特征属性之间的相互关系[5]。目前，一些学者已将Copula 函数应用到水文水资源领域[6-9]，Copula 函数能够独立于随机变量的边缘分布来反映变量之间的相关性，不限定变量的边缘分布，能够通过边缘分布和相关性结构两部分构建多维联合分布，并且在转换过程中不会改变原始随机变量的信息[10-12]。针对干旱事件的研究往往采用游程理论从长期实测资料中提取干旱指标（如干旱研究中常用的干旱历时、干旱烈度、干旱强度等），然后进行Copula 拟合，从而分析研究区域的干旱分布状况[13]。

以往的方法在系统刻画研究区域干旱特性，分析区域水文干旱、气象干旱特性方面作用显著[14]。但由于作物生育期是作物生长需水的关键时期，且作物生育期缺水对作物产量影响非常严重，而上述方法不能得出作物生育期丰枯遭遇情况，因此，该方法在分析农业干旱方面稍显不足。本文基于黑龙江省作物生育期降雨量和参考作物腾发量之间的天然相关性，收集研究区1980-2017年的气象资料，根据Penman-monteith公式求得黑龙江省生育期参考作物腾发量，通过比较简单易行且已在水文计算中得到成功应用的Copula 方法[15，16]构建作物生育期降雨量和参考作物腾发量联合分布模型，利用Kolmogorov-Smirnov 检验和最小离差平方和（OLS）方法进行拟合检验；利用联合分布模型对其丰枯遭遇频率进行研究，以期为定量评价灌溉系统的干旱风险提供基础依据和科学参考，从而促进粮食增产。

1 数据与方法

1.1 研究区概况与数据来源

黑龙江省位于中国东北部，地理范围为东经121°11'～135°05'，北纬43°26'～53°33'，东西跨越14 个经度，根据黑龙江省耕地和水体的空间分布情况，本文选取黑龙江省富裕和海伦2 个代表性站点进行研究（图1）。黑龙江省属于寒温带与温带大陆性季风气候，全省的主要气候特征为：春季低温干旱，夏季温热多雨，秋季易涝早霜，冬季寒冷漫长，无霜期短，气候地域性差异大。降水的季风性特征明显，具体表现为夏季降水充沛，冬季干燥少雨，年降水量一般处于400～650 mm 之间[17]。空间上，中部山区降水较多，西部和北部较少。全省多年平均气温-5～5℃。年太阳辐射总量在4.4×109～5.0×109J/m3之间[18]。全年约5 000～7 000 h 的风速不小于3 m/s，风能密度为200 W/m3[19]。该地区的耕作制度为一熟制，主要作物有春玉米、大豆、春小麦等，其中玉米产量占全国产量的比例为23%，本文以春玉米生育期为代表进行研究，即5-9月[20]。

本文的日尺度降水、最高气温、最低气温、相对湿度、日照时数和平均风速数据库由中国气象局提供。数据处理主要参照全球观测系统世界气象组织指南（World Meteorological Organization’s guide to the Global Observing System）和中国气象局天气观测技术规范标准，数据下载网站参照http：//cdc.cma.gov.cn。其中日照时数通过公式（参考FAO）转化为日太阳辐射量为计算参考作物腾发量所用。作物生育期数据来自于中国气象局的中国农作物产量资料旬值数据集。本文研究所需数据的时间范围见表1。

表1 研究数据描述Tab.1 Description of study data

1.2 研究方法

1.2.1 参考作物腾发量

1998年FAO推荐的Penman-monteith 公式，所使用的数据包括测站的高程、纬度、风速测量、高度、日最高气温、日最低气温、日平均气温、日平均风速、日平均相对湿度和日照时数等。其具体形式如下：

式中：ET0为参考作物腾发量，mm/d；∆为饱和水汽压—温度曲线上的斜率，kPa/℃；Rn为植物冠层表面净辐射，MJ/(m2·d)；G为土壤热通量，MJ/(m2·d)，逐日计算G=0；γ为湿度计常数，kPa/℃；u2为2 m 高处的风速，m/s；es和ea分别为饱和水汽压和实际水汽压，kPa；T为日平均气温，℃。

1.2.2 二维Copula 函数

Copula 方法的理论基石是Sklar 定理：设X和Y分别为降雨量、参考作物蒸腾量，其联合分布函数为F(x,y)，边缘分布函数分别为Fx(x)和Fy(y)，则存在唯一联合函数Copula 函数Cθ(u,v)，满足：

Copula 理论联合分布的建立具体包括以下几个步骤：①确定春玉米生育期降雨量和参考作物腾发量的边缘分布；②基于生育期降雨量和参考作物腾发量数据，计算Kendall 秩相关系数τ；③根据τ与θ的关系计算θ；④根据评价指标选取Copula 函数，建立联合分布；⑤根据所建分布模型进行丰枯遭遇频率分析。水文领域常用的Copula 函数中的3 种：Frank Copula、Clayton Copula 和Gumbel-Hougarrd Copula 函数[21，22]，本文利用这3 种函数构建二维干旱特征变量的联合分布（表2）。

表2 Copula函数类型Tab.2 Copula function types

1.2.3 拟合检验和拟合优度评价指标

本文采用Kolmogorov-Simirnov(K－S)检验对Copula 函数进行拟合检验，采用离差平方和(OLS)最小准则对Copula 函数进行拟合优度评价，其中K－S检验统计量D和离差平方和OLS的定义如下：

式中：ck为联合观测值样本(xk,yk)的Copula 值；mk为联合观测值样本中满足条件x≤xk且y≤yk的联合观测值的个数。

式中：Pi和Pei分别为联分布的理论频率和经验频率。

式中：Num表示联合观测值对小于等于值对(xi,yi)的个数；N为联合观测值的总数。

1.2.4 联合分布概率与重现期

以pf＝37.5%，pk＝62.5%作为对富裕和尚志春玉米生育期降雨量与参考作物腾发量进行丰枯划分的频率[23-25]，运用已构建的二维Copula 联合分布模型，对二者在年际间的丰枯遭遇情况进行分析，两个站点的丰枯遭遇理论上可分为9 种情况[15]，其中，X表示降雨量序列，Y表示参考作物腾发量序列：

丰丰型：p1=p(X≥xpf,Y≥ypf)；丰平型：p2=p(X≥xpf,ypk＜Y＜ypf)；丰枯型：p3=p(X≥xpf,Y≤ypk)；

平丰型：p4=p(xpk＜X＜xpf,Y≥ypf)；平平型：p5=p(xpk＜X＜xpf,ypk＜Y＜ypf)；平枯型：p6=p(xpk＜X＜xpf,Y≤ypk)；

枯丰型：p7=p(X≤xpk,Y≥ypf)；枯平型：p8=p(X≤xpk,ypk＜Y＜ypf)；枯枯型：p9=p(X≤xpk,Y≤ypk)。

考虑到本文主要研究的是黑龙江省春玉米生育期丰枯遭遇风险研究，所以重点关注P-ET0枯丰型的两种联合分布概率及其重现期，即：

式中：P(X≤x or Y≥y)指X不超过某特定值、Y超过某特定值这两个事件中至少有一个事件发生的联合分布概率；P(X≤x,Y≥y)指X不超过某特定值且Y超过某特定值时事件发生的联合分布概率。Tx,y为联合重现期，表示两事件至少有一个发生的重现期；T'x,y为同现重现期，表示两事件同时发生的重现期。

2 结果与分析

2.1 P、ET0边缘分布与丰枯划分结果

图2为富裕和海伦1980-2017年春玉米生育期降雨量与参考作物腾发量系列，整体上，富裕P-ET0差值波动幅度大于海伦。其中，富裕的参考作物腾发量均大于降雨量，仅1998年和2014年例外；海伦1980-2017年的参考作物腾发量波动比较稳定。

假设春玉米生育期降雨量和生育期参考作物腾发量服从Gamma 分布，Gamma 分布与P-III 型分布同属于Γ 类分布，是P-III的一种特殊形式，在我国降水研究中适用性强[15]。取KS检验的显著性水平α=0.05,n=38 时，对应的分位点D0=0.215 44，富裕的P和ET0对应的D分别为0.131、0.081，海伦的P和ET0对应的D分别为0.153、0.089，均小于D0，且对应的PK-S值均大于0.05，说明假设成立（表3），且由Gamma 分布可得两个站点春玉米生育期P和ET0的丰枯划分数值（表4）。

2.2 Copula函数评价结果

表3 单变量边缘分布参数估计拟合检验Tab.3 Univariate marginal distribution parameter estimation fit test

表4 春玉米生育期降雨量和ET0的丰枯划分值Tab.4 Precipitation and ET0 during the growth period of spring maize

富裕和海伦春玉米生育期降雨量和参考作物腾发量的Kendall 秩相关系数τ分别为-0.391、-0.511，且对应的PKendall小于0.05，说明春玉米生育期降雨量和参考作物腾发量之间具有显著的负相关性，由τ与θ的关系（表2）可知，只有Frank Copula 函数可作为春玉米生育期降雨量和参考作物腾发量的联结函数（表5），其计算公式为：

式中：u、v为降雨量和参考作物腾发量的边际分布函数。

表5 二维Copula函数参数值及拟合优度检验Tab.5 Two-dimensional Copula function parameter value and goodness-of-fit test

对Frank Copula 函数进行K-S 检验，富裕和海伦的检验统计量D分别为0.104、0.092，均小于D0，OLS分别为0.039、0.041，且Frank Copula函数得出的理论频率能够较好地与经验频率拟合，两个站点经验频率与理论频率的R2均在0.9 以上（图3），说明由Frank Copula 函数建立的联合分布模型对春玉米生育期降雨量和参考作物腾发量的联合分布情况拟合较好，联合分布理论频率分布模型比较合理。两个站点春玉米生育期降雨量与参考作物腾发联合分布模型分别为：

（1）富裕：

（2）海伦：

2.3 P-ET0丰枯遭遇联合分布频率分析

在9 种丰枯遭遇组合中，P-ET0枯丰型组合是最不利于灌溉调度的情况，图4分别绘制了富裕和海伦春玉米生育期PET0联合分布GX,Y(x,y)、G'X,Y(x,y)等值线，以及两站点1980-2017年春玉米生育期降雨量和参考作物腾发量在其中的分布，GX,Y(x,y)表示P枯、ET0丰这两个事件至少有一个事件发生的联合分布概率；G'X,Y(x,y)表示P枯、ET0丰这两个事件同时发生的联合分布概率。整体上，富裕和海伦联合概率值GX,Y(x,y)、G'X,Y(x,y)变化趋势是相同的，都是随着P的增大和ET0的减小而增大。局部看，富裕[图4（a）]1980-2017年春玉米生育期P-ET0枯丰型中实测值(P≤351.71mmor ET0≥540.01mm)的联合分布理论概率约为0.4，有13年的P-ET0在该概率范围内，占总年份的34.21%，其中7年集中在2000-2010年；(P≤351.71mm,ET0≥540.01mm)的联合分布概率约为0.3，有11年的P-ET0在该概率范围内，占总年份的28.95%，其中6年集中在2000-2010年。海伦[见图4（b）]1980-2017年春玉米生育期P-ET0枯丰型中实测值(P≤437.14 mmor ET0≥507.54 mm)的联合分布理论概率约为0.4，有15年的P-ET0在该概率范围内，占总年份的39.47%，其中1980-1989年和 2000-2010年分别占有5年；(P≤437.14 mm,ET0≥507.54 mm)的联合分布概率约为0.3，有9年的P-ET0在该概率范围内，占总年份的23.68%，其中有3年在1980-1989年，有4年在2000-2007年。

综合富裕和海伦的P-ET0枯丰型联合分布概率，由式（7）计算得出黑龙江省春玉米生育期P枯、ET0丰这两个事件至少有一个事件发生的联合重现期为2.50年，P枯、ET0丰这两个事件同时发生的同现重现期为3.33年。说明在天然来水条件下黑龙江省供水和需水处于不协调状况的概率较高，事件发生的间隔时间较短，不能满足作物需水的风险较大。在农业生产中，需要进行灌溉活动、科学调配水资源。

3 结 语

随着粮食产量的减少、水资源的供给不足以及气候变化等问题的出现，分析作物生育期丰枯遭遇特征，对指导农业防旱抗旱有重要作用。本文以黑龙江省富裕和海伦1980-2017年的春玉米生育期（5-9月）降雨量和参考作物腾发量为研究对象，运用Copula函数构建P-ET0丰枯遭遇联合分布模型，并对其进行分析得出以下结论：

（1）黑龙江省春玉米生育期的降雨量与参考作物腾发量均服从Gamma分布，由Frank Copula函数建立的联合分布模型对春玉米生育期降雨量和参考作物腾发量的联合分布情况拟合较好，其经验频率与理论频率的R2大于0.9。

（2）黑龙江省的联合概率值GX,Y(x,y)、G'X,Y(x,y)均随着P的增大和ET0的减小而呈增大趋势。在天然来水条件下，黑龙江省供水不能满足作物需水的风险较大，且间隔时间较短，春玉米生育期P枯、ET0丰的联合概率值GX,Y(x,y)、G'X,Y(x,y)分别为0.4、0.3，联合重现期和同现重现期分别为年2.50、3.33年。

（3）针对水资源供需不协调、P-ET0枯丰遭遇风险较大且重现期间隔短的情况，可通过降雨量和ET0的不同遭遇组合情况对其干旱风险进行定性描述和定量评估，从而为制定和调整灌溉规划、科学调配水资源提供一定参考依据。