农作物气象灾害风险识别与评估研究 *

赵思健，张　峭，聂　谦，郑茗曦

(1.中国农业科学院农业信息研究所 农业部智能化农业预警技术重点开放实验室，北京 100081; 2.华中科技大学武昌分校 经济管理学院金融学，湖北 武汉 430064)

我国农业生产总体虽呈增长趋势，但农作物产量年际间的波动日益加剧，是普遍存在、不容忽视的问题[1]。农作物产量波动的影响因素大致归结为社会因素、气象因素及其他一些偶然因素，其中气象因素、尤其是气象灾害对作物产量波动的影响最大[2-3]。据统计，自1990年代以来，我国几乎每年都会暴发重大的气象灾害，包括干旱、洪涝、冷害、风雹、台风等。农作物平均每年受灾面积约为46 966 khm2，占总播种面积的30.57%，平均年成灾面积约为24 774 khm2，占播种面积的16.13%，给农作物产量造成很大的影响。为此，亟需对农业生产的气象灾害风险进行有效评估，其中识别造成作物产量波动的关键气象要素，并对气象要素造成产量的异常波动模拟是核心技术。该项技术将能为农业气象灾害的监测预警及农业保险的气象指数产品设计提供重要的支撑。与农业气象灾害风险相关的研究有很多，有代表性的包括杜鹏等较早就提出农业气象灾害风险的体系和特点[4]，李世奎等提出农业气象灾害风险评估体系及模型[5]，李志明等[6]、王丽媛[7]等对农业气象灾害风险评估进行综述，张继权等系统总结农业气象灾害风险评估方法及模型[8]，等。在方法论上，农业气象灾害风险评估方法也层出不穷，有资料统计分析法[9-11]、数学模型法[12-13]、实验机理模拟法[14-16]、作物模型模拟法[17-18]、GIS分析法[19-20]，等。但可以发现，这些研究大都局限于对作物的某种气象灾害风险进行分离式评估，忽略了作物减产损失往往是由于生长期内受到了多种气象灾害的多重影响。因此，农业气象灾害风险评估需要从微观层面上研究作物生长期内各种极端气象事件共同作用造成的作物减产。另一方面，有关气象要素与作物产量之间关系的研究也有很多，它们通常是借助相关分析法[20-23]、因子分析法[24]、多元回归分析法[25-27]、神经网络法[28]等统计方法，直接建立生长期气象要素与作物产量的关系模型，研究目的大都是为作物产量预测服务，在研究视角上非真正意义上的作物产量气象风险研究。本文摆脱直接构建气象要素和作物产量之间关系的传统方法，从作物气象灾害的原理(即气象灾害是气象异常事件、而气象异常事件又导致作物产量异常的原理)出发，通过构建气象灾害指数，识别作物生长期气象灾害指数异常和作物最终产量异常之间的关系来达到评估作物产量的气象灾害风险。由于作物在生长期间面临着多种气象灾害，且各种气象灾害对作物最终产量的影响是综合性的。为解决多种气象灾害对产量的综合影响，采用气象灾害指数异常与产量损失异常之间的关系匹配技术，并通过不断调整异常识别阈值来寻找最优匹配关系，最终形成一套新的作物产量气象灾害风险的识别与评估技术。

1　资料准备

1.1　数据资料

本研究以安徽省来安县水稻和萧县玉米为案例，收集如下2类数据资料开展分析建模工作。

(1)作物单产数据：收集安徽省来安县1980-2015年水稻单产数据，萧县1980-2015年玉米单产数据，数据来源中国农业部县级农村经济统计数据库。

(2)气象要素数据：收集安徽省来安县和萧县的气象站点自1954-2015年气象要素日观测数据，具体的气象要素包括日累计降水、日平均温度、日照时数，数据来源中国气象局。

1.2　生长期资料

通过实地调研，获悉来安县水稻和萧县玉米的种植规律如下：水稻整个生长期从4-9月份，每年4月上旬开始播种，4月中旬开始出苗，4月下旬开始育秧，6月中旬返青，6月下旬开始分蘖，8月上旬开始孕穗，8月中旬开始抽穗，8月下旬开始灌浆，9月上旬开始乳熟，9月中下旬开始成熟收获；玉米整个生长期从6-10月份，每年6月中旬开始播种出苗，6月下旬长出三叶、5月下旬长至七叶，7月中下旬开始拔节，8月开始吐丝，9月上旬开始乳熟，10月上旬开始成熟收获。

2　分析与建模

2.1　气象灾害指数构建

根据历史资料及实地调研，安徽省易发生的农业气象灾害主要包括干旱、洪涝、高温、冷害和寡照(连阴雨)五种。依据气象灾害的定义，结合收集气象站点的日值气象要素数据，计算这五种气象灾害指数，具体如表1所示。

由于作物生长需要经历一个较长的时期，在生长期间很有可能受到同一种气象灾害的多次影响。为了衡量同一种气象灾害在作物生长期间内对作物生产的影响，可以选取作物生长期间内所有发生该种气象灾害事件对应指数的极值或者累计值来综合代表同一类气象灾害的程度，其中极值代表最严重的一次气象灾害程度，累计值则代表多次气象灾害的累计程度，如表2所示。

表1　气象灾害指数表达与解释

表2　气象灾害指数极值与累计值方案

2.2　单产去趋势

影响作物产量的因素很多，主要有三类，①包括种子、化肥、农技水平、经营管理、政策、病虫害控制等在内的非自然因素；②是气象灾害、病虫害等自然因素；③是随机因素，即：

Y=Ytd+Yw+Yδ。

(1)

式中：Y为作物的实际单产(t/hm2)，Ytd为作物的趋势单产(t/hm2)，Yw统称为作物的气象单产(t/hm2)、Yδ称为随机波动单产(t/hm2)。

趋势单产反映了一定历史时期社会发展水平,可通过把年序作为自变量,以某种函数关系去逼近和模拟。趋势单产的拟合方法主要包括回归分析法、滑动平均法和直线滑动平均法。其中，直线滑动平均法是一种线性回归法与滑动平均法相结合的模拟方法，具有较好的逼近和模拟效果[29-30]。

去除趋势单产后，剩余部分的气象单产和随机波动单产统称为波动单产。波动单产是衡量作物产量风险的重要部分，通常用单产波动率为量度，即：

(2)

式中：Vyt表示作物的单产波动率，正值表示增产，负值则表示减产；t表示年份。

2.3　定量关系的优化匹配

灾害学中，气象灾害的实质是极端的天气现象，换句话说是某种气象灾害指数出现异常。所谓异常，就是偏离正常的程度，或者说是偏离正常水平的程度。偏离程度越大，异常越严重、越极端，灾度也就越大。同样，作物产量在天气异常的情况下导致损失。因此，要衡量气象灾害对作物带来的产量风险就可以从建立气象灾害指数异常与产量损失之间的定量关系角度出发进行评估。

其中，作物的产量损失可以通过判断作物单产波动率偏离历史最好增产率的程度来识别，具体表达如下：

dyt=|Vyt-max(Vyt)|。

(3)

式中：dyt表示作物的单产损失率；Vyt表示作物的单产波动率；max(Vyt)表示作物历史的最大增产率；t表示年份。

气象灾害指数采用了极值和累计值两种方案，因此每年都会有两套气象灾害指数数据集，即：

(4)

和

(5)

式中：wm,t表示极值型气象灾害指数集；wc,t表示累计型气象灾害指数集；t表示年份。由于不同的气象灾害指数之间存在着量纲的差异，为此可对各气象灾害指数进行归一化处理，即：

(6)

通常，轻微的气象异常通常构不成气象灾害，只有气象异常达到一定程度才可能诱发气象灾害，才有可能对作物生长产生影响。因此，可以在气象灾害指数上增设一个阈值，超出该阈值的极端异常才判断为气象灾害，才有可能导致作物产量损失。这个阈值并非固定，不同类型气象灾害的阈值有差异，不同地区的阈值有差异，不同类型作物因生理脆弱性的不同对应的气象灾害阈值也有差异。因此，在建模时，可以通过不断微调气象要素阈值来不断匹配气象要素异常与产量损失之间线性关系，再通过寻找最大R2找出最优匹配关系，最终确定各气象灾害指数阈值及相应的关系模型(经过反归一化处理)，即：

(7)

式中：fop(dY,W)表示最终寻找的最优气象灾害指数与产量损失关系模型；dY是产量损失；Wg,m(c)是第g种气象灾害指数；a0是线性回归拟合常数项；ag是气象灾害指数的拟合系数；sWg,m(c)是气象灾害指数的异常判断阈值；R2是相关系数。

2.4　气象灾害风险模拟与评估

利用上述气象灾害指数异常与产量损失的关系模型，结合气象灾害指数的概率分布，便可以对其风险进行评估。利用历史气象要素的观测数据，根据气象灾害指数的定义，可以获得各种气象灾害指数的历史样本，即：

wg,m(c)={wg,m(c),1945,…,wg,m(c),t，…,wg,m(c),2015}。

(8)

式中，wg,m(c)表示第g种气象灾害指数的历史样本集，wg,m(c),t表示该气象灾害指数第t年的样本。利用历史样本集，采用参数估计法模拟获得该气象灾害指数的概率分布函数f(wg,m(c))和分布函数F(wg,m(c))(可借助科学计算软件Easyfit实现)。

通过假设重现期的方式可以模拟气象灾害指数发生异常。为此，利用该气象灾害指数的分布函数模拟出重现期下T0的气象灾害指数值，具体表达如下：

(9)

式中，wg,m(c),T0表示重现期T0下的该气象灾害指数值(极值或累计值)。假设其他气象灾害指数不发生异常，将代入公式7计算出该气象灾害在重现期下的作物产量损失dyg,T0，再代入公式3计算出实际的产量减产率Dyg,T0，即：

(10)

假设一组不同的重现期{T0,T1,…,Tn}，按照公式(11)便可计算出一组各重现期对应的气象灾害指数和作物产量损失，即：

{(T0,wg,m(c),T0,Δyg,T0),(T1,wg,m(c),T1,Δyg,T1),

…,(Tn,wg,m(c),Tn,Δyg,Tn)}。

(11)

将上述信息绘制图表上，可获得作物在该种气象灾害影响下的产量风险曲线。

为衡量该种气象灾害的风险，可计算该气象灾害造成的作物期望减产率Vg，即：

(12)

3　结果分析

3.1　作物产量年际损失分析

利用来安县水稻产量时间序列数据和萧县玉米产量时间序列数据，利用直线滑动平均法对两个地区的作物产量波动进行去趋势分析，获得两个地区的单产波动率年际变化情况，并以历史最佳增产率为准，计算两个地区单产损失率(图1)。从图1中看出，萧县玉米单产波动较来安县水稻剧烈。其中，来安县1994年、2003年水稻产量损失较为严重，萧县1982年、1987年、1988年、1994年、2003年玉米产量损失较为严重。

图1　作物单产波动率和损失率

图2　气象灾害指数年际变化

3.2　气象灾害指数年际变化分析

利用来安县和萧县气象站的气象要素时间序列数据，依据气象灾害指数构建的原理，计算出两个地区的气象灾害指数(通过后续的匹配结果，来安县累计型气象灾害指数匹配度较好，萧县极值型气象灾害指数匹配度较好)，并进行归一化处理，最终获得两个地区气象灾害指数的年际变化曲线(图2)。从图2中可以看出，萧县总体气象灾害指数水平较来安县要高。以归一化值0.6为阈值，高于阈值的判断为突出年份的话，来安县暴雨指数在1991年，干旱指数在1981至1984年、1986年、1991年、1994至1996年、2000年、2012年至2013年，高温指数在2013年，寡照指数在1980年、1999年、2002年至2003年、2008年、2011年、2014年较为突出；萧县暴雨指数在1982年、2000年、2006年，干旱指数在1989年、1998年、2001年，高温指数在1988年、1994年、2013年，冷害指数在1992年、寡照指数在2003年、2008年较为突出。

3.3　最优定量关系模型

在产量年际损失率和气象灾害指数(极值和累计值)年际变化的基础上，采用异常优化匹配的方法建立作物气象灾害易损性关系模型。

针对来安县，通过回归系数R2比较，累计值方案(R2=0.994)较极值方案(R2=0.949)优越，因此选择累计值方案获得关系模型，具体如下：

(15)

式中：洪涝灾害指数的异常阈值为546.9 mm，干旱灾害指数的异常阈值127.3 d，高温灾害指数的异常阈值为18.8℃，冷害灾害指数的异常阈值为5.91℃，寡照灾害指数的异常阈值为173.9 h。

针对萧县，通过回归系统R2比较，极值方案(R2=0.953)较累计值方案(R2=0.902)优越，因此选择极值方案获得关系模型，具体如下：

(16)

式中：洪涝灾害指数的异常阈值为334.16 mm，干旱灾害指数的异常阈值为31.8 d，高温灾害指数的异常阈值为5.72℃，寡照灾害指数的异常阈值为131.04 h。

3.4　历史数据回测分析

利用关系模型，对两县的历史数据进行回测分析，回测结果如图3所示。从图3可以看出，来安县与萧县的历史数据回测结果都比较好。来安县的实际单产损失与模型计算的气象灾害单产损失匹配较好的年份包括1980年、1981年、1988年、1991年、1994年、1999年、2001年、2002年、2003年、2006年、2008年、2010年、2011年、2013年和2014年，其中1994年和2003年属于实际产量损失极值的年份，模型计算单产损失与实际损失匹配度很高；萧县的实际单产损失与模型计算的气象灾害单产损失匹配较好的年份包括1980年、1982年、1987年、1988年、1989年、1993年、1994年、1996年、1997年、1998年、2003年、2005年、2008年、2009年、2010年、2012年、2013年和2014年，其中1982年、1987年、1994年和2003年属于实际损失极值的年份，模型计算单产损失与实际损失匹配度很高，这些充分说明了关系模型对于模拟气象灾害产量损失效果很好。当然，由于模型采用了阈值匹配的方案，对于存在那些实际产量损失、但气象灾害指数不异常的负相关性年份会被剔除在模型外，因此会出现一些年份存在单产损失，但模型计算不出现灾害损失的情况。例如，来安县1982年、1992年，萧县1999年、2007年，存在一定的单产损失，但气象灾害指数却不异常，模型并不进行匹配。这些年份的产量损失有可能是由病虫害等其他自然因素造成，或是单产本身存在的随机损失波动。

图3　历史数据回测

地区气象灾害指数最佳概率分布参数来安县洪涝指数(累计)Wakeby分布a=053658，b=12661，g=001671，d=061687，x=-002475干旱指数(累计)Wakeby分布a=24124，b=54723，g=01746，d=-02320，x=-00160高温指数(累计)Wakeby分布a=0，b=0，g=02114，d=009729，x=-002303冷害指数(累计)Logistic分布s=008451，m=007356寡照指数(累计)Wakeby分布a=24639，b=22469，g=03851，d=-03048，x=-005003萧县洪涝指数(极值)GenLogistic分布k=01763，s=00849，m=02920干旱指数(极值)GenExtremeValue分布k=00802，W=01550，m=01890高温指数(极值)JohnsonSB分布g=01219，d=06707，l=12542，x=-00178冷害指数(极值)GumbelMax分布s=01517，m=-00159寡照指数(极值)GumbelMax分布s=01451，m=01672

3.5　气象灾害风险模拟

利用来安县和萧县历史62年(1954年至2015年)的气象要素(降水、温度和日照)日观测数据，结合各种气象灾害指数(极值型和累积型)表达形式，构建历史气象灾害指数的时间序列数据样本，利用Easyfit软件对各种气象灾害指数进行概率分布拟合，通过KS检验、AD检验和卡方检验选取各气象灾害指数的最佳概率分布，如表3所示。

倘若要对某一种气象灾害风险进行刻画，可在假设其他气象灾害指数不发生异常的前提下，利用该气象灾害指数概率分布，通过假设重现期的方式模拟各种气象灾害异常，并利用最优关系模型计算出该异常下作物产量减产率。如假设一组重现期(2年至200年)，通过上述步骤便可获得一组重现期下作物该气象灾害的减产率，绘制成作物的气象灾害风险曲线，如图4所示(以来安县的洪涝灾害和萧县的干旱灾害为例)。图4可以看出，来安县水稻洪涝灾害风险较大，200年一遇的减产风险可以达到28%，而萧县玉米干旱风险在200年一遇的减产风险才达到4.1%；然而，萧县玉米干旱风险的发生条件要低，5年一遇就可以导致减产，但减产率较低；来安县玉米洪涝灾害的发生条件较高，20年一遇才可能导致减产，但减产率却较高。因此，萧县玉米干旱灾害属于高频率但低损失的风险，而来安县水稻洪涝则属于低频率但高损失的风险。

3.6　气象灾害风险量度与排序

通过对各类气象灾害风险模拟后，计算获得两个县各种气象灾害造成的作物单产期望减产率作为一个统一的风险度量值，如图4所示。图中可以看出，来安县水稻的洪涝灾害风险最大，单产期望减产率是4.00%；其次是高温和寡照风险，它们的单产期望减产率是1.39%和1.01%；再次是冷害和干旱，它们的单产期望减产率是0.76%和0.22%。萧县玉米的寡照和干旱风险较大，它们的单产期望减产率是3.91%和2.64%；再次是高温和洪涝，它们的单产期望减产率是1.03%和0.88%；冷害基本没有影响。

图4　作物气象灾害风险曲线

图5　各类气象灾害单产期望减产率

4　讨论

本文以安徽省来安县水稻和萧县玉米为例开展作物气象灾害风险识别与评估研究，在如下几个方面有所突破。

(1)将气象灾害表达成气象要素的指数形式，通过指数可以快速地识别气象灾害的发生和程度。由于区域之间的差异，灾害指数对作物生长的影响有可能是极端效应、也有可能是积累效应，因此气象灾害指数构建时存在极值型指数和累计型指数两种方案，通过评价选择较好的一种方案。

(2)作物生长期间可能会受到多种气象灾害的综合影响，仅考虑一种气象灾害对产量损失进行建模，是不科学的。本研究将作物生长期中潜在影响的所有气象灾害与产量损失之间关系进行建模，科学量化出多种气象灾害对作物最终产量的综合影响，实际效果较显著。

(3)摆脱直接建立气象灾害与产量关系的传统方法，立足气象灾害源于气象异常、气象异常导致产量损失的原理，通过设置气象灾害指数异常判断阈值，结合数据优化匹配的思路来构建气象灾害指数与作物产量损失，为建立气象与产量的关系模型提供新的思路。

(4)造成作物产量波动包括政策等非自然因素，气象灾害和病虫害等自然因素和随机波动因素。本研究可有效地分离非自然因素带来的趋势产量、自然因素带来的自然波动产量和随机因素造成的随机波动产量，为全面认识作物产量波动开展风险评估提供了技术支撑。

(5)通过异常优化匹配构建的气象灾害指数与产量损失的关系模型，可以有效地模拟与评估出各种气象灾害对作物产量损失的影响，并对各种气象灾害风险进行量度与排序，为农业气象灾害风险评估和农业气象指数保险等应用领域提供了技术支撑。

与此同时，研究本身也存在一些不足，需要在接下来的工作中加以改进与完善，具体包括：① 每种气象灾害指数都有极值型和累计型两种方案，不同气象灾害指数之间很有可能存在极值型和累计型方案的组合，再与作物产量进行优化匹配；② 作物生长期中存在关键生育期，关键生育期内的气象灾害对作物产量波动影响较显著，因此建模时还需要从作物生育机理出发对气象灾害指数进行筛选和构造。

参考文献：

[1]林燕，于冷. 中国粮食产量波动分析[J]. 吉林农业大学学报，2002，28( 3) : 346-350.

[2]张宇. 近40年来我国粮食产量变化特征初步分析[J]. 中国农业气象，1995，16 (3) : 1-4.

[3]王玉斌，蒋俊朋，王晓志，等.中国粮食产量波动影响因素实证分析[J].北京农学院学报, 2007, 22 ( 4) : 38-41.

[4]杜鹏，李世奎.农业气象灾害风险分析初探[J].地理学报，1998，63(3)：202-208.

[5]李世奎, 霍治国, 王素艳，等.农业气象灾害风险评估体系及模型研究[J].自然灾害学报，2004, 13(1): 77-87.

[6]李志明，臧俊岭，焦仁庆.农业气象灾害风险评估研究综述[J].现代农业科技，2009(14): 269-270.

[7]王丽媛，于飞.农业气象灾害风险分析及区划研究进展[J].贵州农业科学，2011,39(11): 84-88.

[8]张继权，李宁.主要气象灾害风险评价与管理的数量化方法及其应用[M].北京：北京师范大学出版社，2007.

[9]王春乙，娄秀荣，王建林.中国农业气象灾害对作物产量的影响[J].自然灾害学报，2007,16(5): 37-43.

[10] 王光宇，安徽省农业气象灾害粮食损失率评估[J].农业系统科学与综合研究，2008,24(3)：330-332.

[11] 陈怀亮，邓伟，张雪芬，等.河南小麦生产农业气象灾害风险分析与区划[J].自然灾害学报，2006,15(1)：135-143.

[12] 王积全，李维德.基于信息扩散理论的干旱区农业旱灾风险分析——以甘肃省民勤县为例[J].中国农业气象，2007,28(4)：440-442.

[13] 潘护林.山东省农业气象灾害灾情统计特征与灰色关联分析[J].西北师范大学学报(自然科学版)，2008,44(5)：94-98.

[14] 刘玲，郭建平，高素华.低温、干旱并发对玉米影响的评估研究[J].气象，2006,32(4)：116-120.

[15] 高素华，刘玲.低温、干旱胁迫对抽雄期玉米叶片光化效率和光合作用速率的影像[J].气象，2007,33(4):88-91.

[16] 陈振林，张建平，王春乙，等.应用WOFOST模型模拟低温与干旱对玉米产量的综合影响[J].中国农业气象，2007,28(4)：440-442.

[17] 董姝娜, 庞泽源, 张继权, 等.基于CERES Maize模型的吉林西部玉米干旱脆弱性曲线研究[J].灾害学, 2014,29(3):115-119.

[18] 贾慧聪, 王静爱, 潘东华, 等.基于EPIC模型的黄淮海夏玉米旱灾风险评价[J].地理学报, 2011,66(5):643-651.

[19] 罗培.GIS支持下的气象灾害风险评估模型[J].自然灾害学报，2007,16(1)：40-43.

[20] 赵霞，王平，龚亚丽，等.基于GIS的内蒙古中部区域洪水灾害风险评估[J].北京师范大学学报(自然科学版)，2007,43(6)：666-669.

[21] 杨晓华，曹天堂，杨小利.影响平凉金果产量的生态气象因子及其变化趋势[J].安徽农业科学，2009，37(27): 13033-13036.

[22] 罗丽华，陈桂华，胡英，等.气象因素与早稻产量因子的相关性分析[J].自然资源学报，2010，25(10): 1718-1726.

[23] 陈彦虎，顾宁，刘静，等.宁夏石嘴山市枸杞产量与气象条件的关系[J].安徽农业科学，2012，40(27): 13508-13511.

[24] F M Chmieleki, W.Kohn.Impact of weather on yield components of winter rye over 30 years[J].Agricultural and Forest Meteorology.2000, 102: 253-261.

[25] 彭九慧，杨梅，陆倩，等.玉米干旱等级划分及气象产量评估方法[J].环境科学与技术,2012,35(12):137-139.

[26] 范磊，郑国清，赵全志，等.河南省水稻气象灾害因子分析及主要灾害发生区划研究[J].2013，34(6): 162-168.

[27] M G Pereira, L.Caramelo, C.Gouveia.Assessment of weather-related risk on chestnut productivity[J].Nat.Hazards Earth Syst.Sci., 2011, 11: 2729-2739

[28] 乔加新，周森鑫，马季.基于BP神经网络的农业气象产量预报系统[J].微计算机信息, 2009(35): 44-46.

[29] 白林，段东霞，万 忠.广东水稻生产风险评估研究[J].西南农业学报，2014，27(4) : 1762-1767.

[30] 王克.农作物单产分布对农业保险费率厘定的影响[D].北京：中国农业科学院, 2008.