基于多维Copula函数的松嫩平原干旱特征分析

郭嘉豪，王会肖，赵茹欣，宁理科

(1.北京师范大学水科学研究院，北京 100857；2. 中国科学院地理科学与资源研究所生态系统网络观测与模拟重点实验室，>北京 100101； 3. 中国科学院禹城综合试验站，北京 100101)

0 引 言

干旱，作为一种自然灾害，长期困扰着社会经济发展和人民生活。同时由于我国地形、山脉和逐年间不稳定的季风气候造成的水热分布不均进一步导致了我国干旱的频繁发生[1]，且干旱造成的损失逐年增加[2]，因此对干旱进行频率分析十分必要。干旱事件一般拥有多个方面的特征属性，如干旱历时(D0)、干旱烈度(S0)和干旱峰值(P0)等，单一分析单个干旱特征变量则会忽略干旱特征变量间的相互联系及影响，无法全面反映干旱事件的真实性，不能完整地评估干旱事件的发生概率[3]，所以需综合考虑多个变量，才能更好地理解和认识到不同条件下干旱事件发生可能性的大小。

传统的多变量干旱分析计算方法有多元概率函数法和非参数法，如，Beersma和Adri Buishand[4]应用二维正态分布建立了荷兰降雨不足量和莱茵河流最不足量的二维联合分布，分析了荷兰干旱事件的特点；Kim等[5]用二维非参数核密度估计法进行了干旱的特征分析。但传统的多变量分析大多都要求边缘分布具有相同的分布函数，且大多只适用于二维，很难向高维转换[6]。为解决上述问题，国内外学者将Copula函数应用到干旱事件分析计算领域[7]。Copula函数的优点在于不限定变量的边缘分布，能够通过边缘分布和相关性结构两部分构建多维联合分布，同时求解简单[8-9]。陆桂华等[10]采用Copula函数，建立了区域干旱历时和干旱强度的联合分布并计算了联合分布的重现期；闫宝伟等[11]以汉江上游为研究对象，利用标准化降水指数(Standardized Precipitation Index，SPI)，通过Copula函数构建了干旱历时和干旱程度的联合分布，更加全面地反映了干旱特征；周玉良等[12]利用Archimedean Copula函数构建了干旱历时与干旱烈度的联合分布并计算干旱重现期，得到的结果与区域实际受旱状况相符；李计等[13]通过新疆乌鲁木齐和石河子气象站数据并基于Archimedean Copula函数分别构建二维、三维干旱变量的联合分布，表明Copula函数能够描述多维干旱特征。

东北地区是我国主要粮食产区之一，耕地面积占全国18.1%，粮食总产量占全国16.5%[14]，但极易遭受干旱影响；同时2020年由于疫情影响，多国限制了粮食出口，粮食安全更加受到威胁。松嫩平原是东北三大平原之最，然而目前利用Copula函数来研究松嫩平原干旱特征的研究还鲜有报道，且由于Copula函数的“维数诅咒”[7]，大多数研究主要集中在二维Copula函数的应用方面。因此，本文选取松嫩平原为研究区，计算3个月时间尺度的SPI，根据游程理论提取干旱历时、干旱烈度和干旱峰值3个干旱特征变量，最后通过Copula函数构建二维和三维联合分布，推求多维变量的组合重现期并以此分析松嫩平原的干旱特征，以期能为松嫩平原防旱抗旱以及保障粮食安全提供理论基础。

1 研究内容与方法

1.1 研究区概况及数据资料

松嫩平原，位于黑龙江省西南部和吉林省西北部，主要由松花江和嫩江冲积而成，平原面积为21.6 万km2。地理位置为121°38′～128°33′ E，42°12′～49°12′ N。研究区四季气候变化明显，为半干旱半湿润交界地带，全年降水量为530～645 mm，主要集中在夏季，多年平均蒸发量为1 250～1 650 mm[15,16]。

本文选取松嫩平原21个气象站点1961-2015年逐月降水数据(国家气象信息中心，http:∥data.cma.cn/)，各站点分布见图1。

图1 松嫩平原气象站点分布图Fig.1 The location of meteorological stations in Songnen Plain

1.2 研究方法

本文利用游程理论从SPI序列中提取出干旱历时、干旱强度和干旱峰值，并利用韦伯(Weibull, WEI)分布、伽马(GAMMA)分布、指数(Exponential, EXP)分布和对数正态(Logarithmic normal, LOGNOR)分布对干旱历时、干旱强度和干旱峰值进行边缘分布拟合，基于Copula函数建立二维、三维联合分布函数，以赤池信息量准则(AIC)选择最优拟合分布函数，并对研究区干旱的重现期进行研究分析。

1.2.1 干旱特征变量的识别

图2 干旱特征变量的定义与提取示意图Fig.2 Schematic diagram of definition and extraction of drought feature variables

1.2.2 干旱特征变量的边缘分布和联合分布模型

采用WEI分布、GAMMA分布、EXP分布和LOGNOR分布分别拟合干旱历时、干旱烈度和干旱峰值的边缘分布，采用极大似然法估算参数，再利用K-S检验(Kolmogorov-Smirnov test)从中选择最优边缘分布[20]。

WEI分布函数为：

(1)

式中：a为比例参数；b为形状参数。

GAMMA分布函数为：

(2)

式中：α为形状参数；β为尺度参数。

EXP分布函数为：

f(x)=λexp(-λx)

(3)

式中，λ>0。

LOGNOR分布函数为：

(4)

式中：μ为对数均值；σ为对数标准差，其中x>0。

Copula函数不限定变量的边缘分布，通过Copula模型，可以将k个任意的边际分布连接起来，形成一个多变量联合分布概率模型。Copula函数的理论基础是Sklar’s 定理，随机变量X1、X2、…、Xn连续，F1(x1)、F2(x2)、…、Fn(xn)是其边缘分布函数，F为n维联合概率分布函数，则存在Copula函数C[0 1]n→[0 1]，使得：

F(x1,x2,…,xn)=C[F1(x1),F2(x2),…,Fn(xn) ]

(5)

若边缘分布函数F1(x1),F2(x2),…,Fn(xn)是连续的，则Copula函数唯一确定[1]。

水文领域常用的Copula函数中的3种：Frank Copula、Clayton Copula和Gumbel-Hougarrd Copula函数(表1)[21]，本文利用这3种函数构建多维干旱特征变量的联合分布。拟合优度的检验选用赤池信息量准则(AIC)，以AIC最小值作为最优拟合的标准，公式为[22]：

AIC=Nln(MSE)+2k=

(6)

式中：N是样本个数；xc为经验频率；x0为理论频率；k为参数个数；MSE为均方误差。

设N为观测序列样本对数，按x升序排序所得数据系列为(x1,y1)，(x2,y2)，…，(xn,yn),则二维经验概率采用下式[23]计算：

P(xi,yi)=P(X≤xi,Y≤yi)=

(7)

式中：No.of代表样本升序排列后所得的样本序号。

其中，三维经验概率计算方法与二维经验概率计算公式类似[9]，具体如下：

P(xi,yi,zi)=P(X≤xi,Y≤yi,Z≤zi)=

(8)

已知干旱历时、干旱烈度和干旱峰值的边缘分布函数和联合分布函数，根据重现期理论，边缘分布为u的单变量重现期计算公式[20]为：

(9)

式中：N为资料的序列长度；m为干旱事件的次数。

多维变量联合分布组合重现期包括联合重现期和同现重现期，边缘分布为u,v的二维分布联合重现期Ta和同现重现期To。计算公式为：

(10)

(11)

边缘分布为u,v和w的三维联合分布联合重现期Ta和同现重现期To。计算公式为：

(12)

To=

(13)

2 多维干旱变量联合分布在松嫩平原的应用

2.1 干旱特征变量边缘分布

2.1.1 相关性分析

一个矩阵也许很复杂，线性代数中的一个常用的方法，便是把矩阵分解成一些更简单的矩阵的乘积，通过研究这些更简单的矩阵，来得到原始矩阵的一些有用的性质。矩阵的UV分解是线性代数中分解一个矩阵的经典的方法，我们在这里不做详细地解释，但是我们推荐对此感兴趣的读者去参看参考文献[2]中的第二章。

采用Pearson和Kendall相关系数计算3个干旱变量的两两相关性，结果列于表2。结果表明，变量均呈现出较强的正相关特性；其中干旱烈度和干旱历时的相关性最强，Pearson相关系数达到0.90，Kendall相关系数为0.76；干旱烈度和干旱峰值的相关性次之，Pearson和Kendall相关系数分别为0.57和0.70。

表2 干旱事件三变量Pearson 和 Kendall 相关系数Tab.2 Pearson and Kendall correlation coefficients of three variables of drought events (SPI<0)

2.1.2 边缘分布函数的推求

松嫩平原3个干旱特征变量边缘分布的拟合结果如表3所示。

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

由表3可知，在显著性水平α=0.05的条件下，松嫩平原干旱历时边缘分布符合WEI分布，通过K-S检验；干旱烈度和干旱峰值的边缘分布符合LOGNOR分布。通过相关性分析可知，3个变量存在较高的相关性，可利用Copula函数进行分析。

2.2 基于Copula函数的多维干旱变量联合分布

目前Copula函数在干旱研究中的应用，大多集中在干旱历时和干旱烈度的二维联合分布模型的构建，但干旱特征变量还包括干旱峰值等[24]，同时干旱是一个多变量事件，二维的联合分布无法全面反映干旱的特征。干旱峰值能反映干旱事件中最严重的状态，因此构建干旱历时、干旱烈度和干旱峰值的三维联合分布模型可以更加全面反映干旱的特征。

利用极大似然法计算Copula函数的参数θ，通过函数计算理论频率与经验频率之间的均方误差(MSE)，以最小AIC为选择最优拟合函数的标准。结果如表4所示，对于干旱历时和干旱峰值以及3变量的分布，Frank Copula具有最小的AIC，拟合优度最好；干旱历时和干旱烈度中，Gumbel Copula的AIC最小，具有最好的拟合优度；干旱烈度和干旱峰值中，Clayton Copula函数为最优拟合函数。

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

根据上述评估得到的不同变量组合所服从的最优Copula函数，以经验累积概率为横坐标，相应的理论累积概率为纵坐标绘制Q-Q图(图3)，发现对应的点均匀分布在y=x的两侧，表明所选函数整体上拟合效果良好。

图3 联合分布经验累积概率与理论累积概率对比Fig.3 Comparison of cumulative experience probability and theoretical cumulative probability

分别以干旱历时、干旱烈度和干旱峰值为横纵坐标，绘制两两间干旱特征变量二维联合分布图(图4)。由图4可知，当干旱历时、干旱烈度和干旱峰值都相应增大时，二者的累积概率也都呈相应的增大趋势；其中，当干旱特征变量较小时，其两两对应的累积概率的增长速率相对稳定，但当干旱历时大于8，干旱烈度大于7和干旱峰值大于2时，其两两相应的联合概率增长速率相对减缓。可见，同时满足干旱长历时、强烈度、高峰值的干旱事件的发生概率较小。当然，及时采取措施是避免干旱缓慢发展为旱灾的必要手段。

图4 松嫩平原干旱特征二维联合分布图Fig.4 Two-dimensional joint distribution map of drought characteristics in the Songnen Plain

2.3 松嫩平原干旱重现期分析

在水文分析中，确定某一稀遇事件或重大灾害发生的概率，具有极其重要的现实意义。在工程中常用重现期来度量风险的大小。所谓重现期就是指某随机变量的取值在长时期内平均多长时间出现一次[1]。重现期在水文中被视为一个通用的标准，它提供了一种简单实用的方法来评估水文风险。然而，联合重现期和同现重现期在干旱的多变量频率分析中，没有得到很好的应用。

分别绘制在二维联合概率为0.1～0.9的干旱历时、干旱烈度和干旱峰值两两变量的联合重现期(图5)，其单位为a。由图5可知，随着干旱特征变量不断增加，两变量的联合重现期也随之增大，且联合重现期均随着两干旱特征变量的增加呈现先陡后缓的增加趋势。干旱特征变量间的变化也会对联合重现期产生影响，由图5(a)和图5(b)可以看出，在干旱历时保持不变的情况下，干旱峰值所导致联合重现期增加的影响要大于干旱烈度，例如在当干旱历时为4个月时，干旱峰值由1到2，联合重现期从4.6 a增加到了10.8 a，而干旱烈度从1到2，联合重现期仅由4.1 a增加到了5.4 a；从图5(b)和图5(c)可得，干旱烈度对于联合重现期的影响要略大于干旱历时；由图5(a)和图5(c)可以看出，当干旱烈度保持不变时，同等程度的干旱峰值变化导致的重现期增加值大于同等程度干旱历时变化引起的联合重现期增加值。综上可知，松嫩平原干旱峰值的变化对联合重现期的影响要大于干旱历时和干旱烈度。

图5 松嫩平原干旱二维联合重现期等值线图Fig.5 The contour map of the 2D joint recurrence period of the Songnen Plain

分别计算单变量重现期为5、10、15、30、50、100、300 a条件下3个干旱特征变量的边缘分布函数值。根据边缘分布的反函数推求对应的3个干旱特征变量值，并利用该干旱特征变量值计算最优拟合二维和三维Copula函数值。将Copula函数值带入式(10)～(13)中，求出相应的多维组合重现期，结果如表5所示。

表5 干旱特征值及其对应的重现期Tab.5 Drought characteristic value and its corresponding return period

由表5分析可知，随着单变量重现期的增加，联合重现期和同现重现期都在增加，且联合重现期的增长速率没有同现重现期增长速率大，当单变量重现期大于15 a时，两者差异更加明显。单变量重现期(T)，联合重现期(Ta)和同现重现期(To)三者之间有明显的大小关系，为Ta15 a时，D0-S0-P0的To最大。

从表5可知，发生干旱历时达4个月的重现期大约在10～30 a，表明松嫩平原的干旱一般都发生在一个季节内，较少发生跨季节性的干旱事件，要做好季节内防旱措施。发生特旱(P0>2)的重现期约140 a，表明松嫩平原的干旱情势相对较好，特旱事件较少发生，相对多发生轻旱和中旱。以单变量重现期为100 a为例，其干旱历时为7.06个月，干旱烈度为8.91，干旱峰值为2.44，该地区发生干旱历时为7.06个月且干旱烈度为8.91的重现期大约为770.47 a；发生干旱历时为7.06个月且干旱峰值为2.44的重现期大约721.25 a；发生干旱烈度为8.91且干旱峰值为2.44的重现期大约为128.38 a；发生干旱历时为7.06个月，干旱烈度为8.91且干旱峰值为2.44的重现期大约为1 409.59 a 。说明松嫩平原发生长历时、强烈度和极度干旱事件的概率相对较小。这可能与北方地区近些年来，降水增加和蒸发相对减少的原因有关[2]，同时由于温度上升，轻度和中度干旱依旧会发生，对松嫩平原造成一定危害。

3 结 语

松嫩平原所处地区是水汽辐射区，同时易受孟加拉湾和日本海的水汽带影响，易发生干旱事件，也是我国重要的商品粮基地。干旱特征及其发生概率等对松嫩平原的粮食安全有较大影响。目前关于松嫩平原的研究多集中于时空变化和干旱影响因素方面[16,25]，较少有学者对松嫩平原的重现期有一个整体的研究，本文利用多维Copula函数分析了松嫩平原的干旱重现期，评价了松嫩平原的干旱风险，可为松嫩平原防旱抗旱及其干旱研究提供一定的基础。本研究利用多干旱特征变量对重现期进行分析，发现在单变量重现期一定时，其相应同现重现期和联合重现期可作为其取值范围，同时通过两变量干旱特征变量的分析可知，当干旱特征变量较小时，其联合概率增长稳定，但干旱特征变量到达一定程度后，其联合概率增长较为缓慢，这与前人的研究一致[20,21]；同时发现松嫩平原较少发生严重的干旱事件，这也与前人研究相符[16,26]。

本文利用泰森多边形法则得到的面降雨量计算松嫩平原1961-2015年3个月时间尺度的SPI，基于游程理论提取研究区的干旱历时、干旱烈度和干旱峰值3个干旱特征变量，利用Copula函数构建二维和三维的干旱特征联合分布，分析研究区的组合重现期特征。主要结论如下：

(1)干旱历时、干旱烈度和干旱峰值两两呈现较强的正相关，WEI分布对干旱历时的拟合最好，LOGNOR对干旱烈度和干旱峰值的拟合最好。

(2)Gumbel Copula函数对干旱历时和干旱烈度二维干旱特征变量的拟合度最好，Clayton Copula函数对干旱烈度和干旱峰值二维干旱特征变量的拟合度最好，Frank Copula函数对干旱历时和干旱峰值以及3维干旱特征变量的拟合度最优。干旱历时大于8个月，干旱烈度大于7和干旱峰值大于2时，两两的累积联合概率增长速率变缓。

(3)组合重现期可视为单变量重现期的取值范围，其中联合重现期可视为相应单变量重现期的下阈值，同现重现期可视为相应单变量重现期的上阈值；两变量干旱特征变量的联合重现期整体上呈现先陡后缓的趋势，其中干旱峰值对联合重现期的影响最大；在单变量重现期一定的情况下，两变量重现期中，干旱历时和干旱烈度的同现重现期最大，联合重现期最小，其相关性最强。

(4)松嫩平原发生长历时强烈度高峰值的干旱事件的概率虽然较低，同时该地区长历时的干旱事件，其烈度较强，但峰值不高，但发生短历时、中度干旱的概率较高。

干旱是一种复杂现象，受多种因素干扰，干旱指标仅能在一定程度上量化干旱特征，不能完全表征干旱全部特征。本文采用SPI指数来识别干旱，但SPI在计算过程中，仅考虑了降水这一单一因素，未考虑温度、蒸发和土壤含水量等与干旱有关的因素，故在表征干旱特征时可能存在一定的误差；此外，本文以松嫩平原整体情况为研究重点，未考虑松嫩平原的空间差异性，所以今后可以以不同站点进行相应的分析，也可以对干旱特征的更高维进行研究。