冯介玲，李 宁，刘 丽，陈 曦，白 扣

(1. 北京师范大学 地理科学学部 环境演变与自然灾害教育部重点实验室, 北京 100875；2. 北京师范大学 地理科学学部 减灾与应急管理研究院, 北京 100875)

灾害发生频率研究的基本数学手段是概率统计方法，影响灾害的变量通常很多，变量间的相关性会影响灾害的发生频率。为了全面地研究发生频率的统计规律，需要了解多变量之间存在的各种相关关系。在实际灾害的求相关问题中，传统的最小二乘线性回归方法，前提是随机误差项需要符合正态性假设等，当假定条件不满足时，主观勉强求解会使得数据信息发生偏离，从而影响分析结果，因此不适用于非线性关系下的灾害多变量研究[1]。变量间整体和尾部的相关性大小或者正负方向可能不一致，当多变量的相关关系出现厚尾现象时，以传统相关系数为基础的频率分析及风险测度可能会引起对实际风险的高估或低估。因此完整刻画变量之间线性或非线性相关关系，特别是尾部的相关性特征，对于准确分析灾害的发生频率从而评估灾害的风险具有重要意义。相关结构是变量之间完整相关性信息的表征，可以比较全面地描述变量间的相关程度和相关模式，包含变量在不同取值时的相关性，能够刻画各种形式的相关性特征。

Copula函数模型作为一种不受边缘分布形式限制，不基于线性假设的方法，越来越广泛的被运用在自然灾害相关结构的分析中[2]。然而传统上变量之间的相关结构都是利用一种Copula函数来构建的，绝大多数的Copula函数只能刻画一种相关结构特征，但实际上有些数据并不完全是由一种相关结构所决定，它们有可能是几种相关结构的混合[3-4]。所以仅仅用单个Copula函数进行模型拟合可能会出现信息失真的情况[5]。如对于常见的Archimedean Copula函数，Gumbel Copula无法捕捉在变量分布下尾处的相关性特征；Clayton Copula函数无法刻画分布在上尾处的相关结构；Frank Copula只能刻画上下尾结构对称时的相关性[6]。根据单个Copula函数模型局限性的分析，我们利用更适用于灾害分析的混合Copula函数模型，来克服单一Copula函数模型在刻画复杂相关结构时的不足。混合Copula函数结合了其组成部分单个Copula函数的特点，因而能够更为全面地反映变量间的相关信息，在应用上也更具有灵活性，近年来越来越多地受到国内外学者的青睐[7-11]。

本文基于内蒙古中部的强沙尘暴灾害资料，以“相关结构”为手段解决传统相关分析方法产生的评估偏差问题。通过对比单个及混合Copula函数对灾害变量间相关结构的拟合效果，说明混合Copula模型在相关结构分析中的优势，并且探讨了极端灾害事件对应的变量尾部阈值，为提高灾害发生频率及风险水平的评估精度提供技术支撑。

1 研究方法

1.1 Copula函数

1959年Sklar提出，可以将一个联合分布分解为n个边缘分布和一个Copula函数，这个Copula函数描述了这n个变量间的相关结构[12]。在Copula函数中，由于构造和计算的简便灵活性等特点，单参数Archimedean Copula函数得到了最广泛的应用。三种最常用的Archimedean Copula包括Clayton、Frank和Gumbel Copula。

1.2 混合Copula函数

最常见的一种混合(Mixed)Copula模型，即是将多个Copula函数进行线性组合。混合Copula函数的表达式如下[13]：

(1)

混合Copula函数的参数估计方法有极大似然估计法(EM法)、最小二乘法和矩量估计法等方法。从数据处理的简便性，并满足计算精度要求的角度考虑，本研究选择最小二乘法估算权重系数[14]。

2 实例分析

内蒙古自治区位于我国北部边疆，中部地区是我国沙尘暴的频发区之一，也是沙尘暴灾害研究的重点区域，强沙尘暴集中出现在春季。本文选取内蒙古中部(37.61°～46.78°N，105.24°～119.89°E)为研究区，选择有详细强沙尘暴观测资料以及完整气象资料的17个中国地面气象台站(图1)。该区域每年春季的大风日数较多，沙尘天气发生频繁，主要的地表覆盖类型包括草原、半荒漠灌木以及沙漠。同时该区域的云量普遍较少，卫星遥感观测数据质量较高，能够保证遥感监测该区域植被生长过程的准确性。

图1 研究区气象站点位置示意图

本研究选取内蒙古中部1982-1999和2001-2007年25年中165次春季强沙尘暴事件(2000年的风速数据缺失)。根据沙尘暴的致灾机理及已有研究，选取与沙尘暴发生频率相关性较高的植被返青期和春季大风事件作为模型的特征变量[15]。本研究选取的植被返青期是指植被萌发变绿并开始进行光合作用的时间，即土壤暴露期的结束以及植被快速生长的开始；春季大风事件是指3-5月份距地面10 m处风速测量超过沙尘物质起沙风速阈值时(9 m/s)的记录总和。

本文利用归一化植被指数(NDVI)来提取内蒙古中部的植被返青期[16-18]。其中NDVI产品来自于GIMMS(Global Inventory Monitoring and Modeling studies)工作组基于NOAA(The National Oceanic and Atmospheric Administration)卫星上AVHRR(Advanced Very High Resolution Radiometer, AVHRR)资料。研究中统计大风事件的风速数据来自于中国气象局中国气象科学共享服务网(http://data.cma.gov.cn)发布的全球地面天气资料定时数据集，沙尘天气数据来自中国强沙尘暴序列及其支撑数据集。

图2 特征变量频率直方图

特征变量均值标准差偏度峰度返青期161.642937.87000.15705.7801大风事件44.500032.00480.86181.9638

2.1 灾害变量的特征检验

对研究区的植被返青期和春季大风事件数据进行描述性统计分析，图2显示了返青期和大风事件的频率直方图以及拟合的正态分布曲线，样本的描述性统计数据见表1。由图2和表1可知，与正态分布相比，返青期在均值附近与尾部的概率值更大，峰度更高，尾部更厚，具有尖峰态厚尾分布特征，说明极端值出现的频率高。大风事件具有明显的非对称性，存在一定的右偏，有较多的极端高值。由变量的频率直方图形态及偏度、峰度值初步判断数据不服从正态分布，再利用Matlab软件分别对变量进行Jarque-Bera正态性检验，返回值都为1，则认为数据拒绝服从正态分布的假设，表示返青期和大风事件数据都具有非正态性。因此，变量不满足传统线性相关分析的前提假设，考虑用Copula函数进行相关结构分析。

2.2 边缘分布选择估计

令返青期和大风事件的边缘分布分别为FG(g)和FS(s)。利用Easyfit软件对各变量进行多种概率分布类型的拟合及参数估计，再根据Anderson-Darling检验结果，选取最优的一种概率分布。通过0.01水平下的假设检验得到的单变量最优概率分布类型以及参数结果见表2。

2.3 Copula模型构建

利用极大似然法，选取三种常见的单个Archimedean Copula函数对返青期和大风事件的Copula模型进行拟合，得到单个Copula函数的参数并作为混合Copula函数的迭代初值，应用OLS法，进行混合Copula函数的参数估计，得到混合Copula函数的模型表达式如下：

CMixed(u,v|θ,λ,α)=0.699×CClayton(u,v|θ=0.8525)+0.290×CFrank(u,v|λ=2.5813)+0.011×CGumbel(u,v|α=1.3558)。

(2)

通过计算拟合优度评价指标进行优度检验，结果如表3所示。比较表3中各Copula函数的RMSE，AIC，Bias以及OLS值，根据拟合优度指标越小模型拟合度越好的择优原则，可知混合Copula函数的拟合度优于单个Copula函数，返青期和大风事件之间的相关结构可以用混合Copula函数进行更准确的描述。

从表3可以看出，混合Copula函数中的三个成分Copula函数的相关参数均为正数，说明返青期和大风事件之间具有正的相关结构，在任一概率分布处都为正相关。由于混合Copula函数中Clayton Copula和Gumbel Copula函数的权重系数均不等于0，表明返青期和大风事件之间可能存在非对称的尾部相关结构，而其中Clayton Copula的权重系数最大，表明返青期与大风事件两个特征变量的下尾相关性表现较为明显，即当一个变量的取值变小，另一变量跟着变小的概率比较大。

图3显示了两个特征变量基于混合Copula函数的相关结构模型。从图中可以看出，拟合的混合Copula函数在上尾和下尾的频率都相对较大，符合变量有较多极端值的特征。上尾频率大的特点表示在发生强沙尘暴灾害的情况下，同时出现返青期越晚和大风事件越多的概率越大。两个特征变量在上下尾部表现出一定的非对称结构，尤其在极端下尾呈现出较强的正相关关系。返青期与大风事件之间的尾部相关性能被混合Copula函数清晰地定性刻画出来，为了更准确的分析相关结构特别是尾部相关关系，需要进一步进行基于相关结构的定量测度。

图3 返青期和大风事件基于混合Copula函数的相关结构模型

特征变量边缘分布类型分布函数参数返青期Log-LogisticFG(gα,β,γ)=(1+(βx-γ)α)-1α=64.273β=1308.6γ=-1147.4大风事件WeibullFD(dα,β)=1-exp(-(xβ)α)α=1.2773β=47.884

表3 单个及混合Copula函数的参数估计和拟合优度检验值

2.4 基于混合Copula函数的相关结构分析

根据构建的混合Copula函数模型来分析返青期和大风事件的相关结构，分别用Kendall秩相关系数τ和尾部相关系数λ来表征总体和尾部的相关性[19]。

通过计算得到基于混合Copula函数的Kendall秩相关系数为0.29，表明从总体来看，二者存在一定的正相关关系。

计算单个Copula函数及混合Copula函数的尾部相关系数，结果如表4所示。

表4 尾部相关系数值

由表4可知，基于混合Copula函数的下尾相关系数为0.31，表明返青期和大风事件同时出现极小值的概率为31%，而上尾相关系数为0.0965，表明返青期和大风事件同时出现极大值的概率为9.65%。尾部相关系数的结果说明相关结构确实存在上下尾的差异。

为了更好地分析返青期和大风事件的相关结构，以及更加准确地观察到尾部的相关性和拟合效果，对样本聚集的位置，即联合分布中对角线 上的频率变化进行进一步分析，绘制出相应的频率分布图，四种Copula函数拟合情况如图4所示。

综合表4和图4可以看出，Frank Copula函数不能捕捉到变量之间的非对称的尾部相关结构，而Clayton Copula函数和Gumbel Copula函数都只能反映出变量间复杂尾部相关结构的单个侧面。Clayton Copula 函数会高估变量同时取得极端小值的概率，低估变量同时取得极端大值的概率，而Gumbel Copula函数会低估了变量同时取得极端小值的概率，高估变量同时取得极端大值的概率。而混合Copula函数既能准确捕捉到变量取极端值时的相关关系，也可以准确地刻画变量同时出现极小值时的相关性略强于同时出现极大值时的相关性这样的一种非对称的尾部相关关系，能够反映返青期提前伴随大风事件减少的可能性更高这一特征，所以混合Copula函数在刻画变量间复杂的相关结构上具有独特的优势。

基于混合Copula函数的相关结构可以更加准确地反映出返青期和大风事件在不同取值范围上有不同的相关性表现。极端事件往往存在于概率分布的尾部，根据混合Copula函数相关结构的结果，可以进一步探索变量的尾部阈值，即极端事件相关性较高时对应的变量取值。

2.5 极端事件对应尾部阈值

灾害变量的相关结构通常比较复杂，为了掌握组成整体的各部分相互关联的方式，需要找到相关性高低分化的节点，来探讨不同部分的相关性情况。因此获取尾部阈值，识别极端事件同时出现的可能性较大时对应的变量取值，对于全面分析相关结构非常重要。

最优分割法是一种常用的有序样本聚类方法，其基本原理是通过对样本进行不同的分段划分，找出使得段内部各样本间的差异最小，而段与段之间的差异最大的分割方法，来获得样本的最优分段结果。研究表明，利用最优分割法结合相关结构可以很好地划分变量的上下尾部阈值[20]。由图4显示，两个特征变量的经验联合分布样本聚集的主对角线频率具有下尾、中部和上尾的三段特征，所以本文选择最优三分割方法计算上下尾的分割点及对应的变量原始值，即极端事件同时出现时对应的变量尾部阈值，结果如表5所示。

表5 最优三分割结果

由相关结构和最优三分割结果可知，返青期的下尾阈值为第90 d，上尾阈值为第223 d；大风事件的下尾阈值为9次，上尾阈值为100次。在发生强沙尘暴的情况下，当植被返青期早于90d和春季大风事件小于9次，或者当植被返青期晚于223d和春季大风事件大于100次时，两者的正相关性就会显著提高。如果研究区某一站点观测到返青期晚于223d时，由于和发生较多的大风事件的正相关性显著提高，因此需要特别警惕发生强沙尘暴灾害的风险。

3 结论与讨论

本文基于Copula函数理论，针对内蒙古中部1982-1999和2001-2007年25年间的春季强沙尘暴事件，分别建立植被返青期和春季大风事件的单一和混合Copula函数模型，比较不同Copula函数的拟合效果，并利用最优Copula函数进行相关结构的分析，得到以下结论：

(1)根据拟合优度评价结果，混合Copula函数的拟合精度优于单个Copula函数，适合于构建植被返青期和春季大风事件两个特征变量的相关结构模型；

(2)由返青期和大风事件相关结构的秩相关系数可知，二者有一定的正相关关系。混合Copula函数反映的相关结构表明，两个特征变量具有非对称的尾部关系，且在极端下尾呈现出较强的正相关结构，表明极端低值发生时，返青期和大风事件二者具有更高的相关性，其中一个变量取值变小时对另一变量的影响更大。相关结构可以更清晰地说明不同变量取值阶段的相关情况，增加对尾部的研究能够更全面的分析变量的相关性。

(3)对发生强沙尘暴灾害的内蒙古中部地区而言，根据上下尾部的阈值探索可知，植被返青期早于90d和春季大风事件小于9次，以及植被返青期晚于223d和春季大风事件大于100次存在较高的正相关关系。两者同时取得尾部阈值时对应发生沙尘暴的可能性更高，需要引起关注。

相关结构模型可以作为评估自然灾害发生频率及灾害风险的重要工具，混合Copula函数能够有效改善相关结构的拟合精度，提高尾部相关性分析的准确性，更好地探究尾部极端值处的变量相关性变化情况。关注混合Copula相关结构模型在灾害评估中的应用，有利于提高风险分析的准确度，为优化风险评估模型提供理论依据。