一种帕尔默干旱指数计算方法及系统

摘要

本发明涉及一种帕尔默干旱指数计算方法及系统，该方法包括：采集多个气象站点的气象数据，气象数据包括纬度、年份、温度、降水量和可利用水量；根据多个纬度确定各气象站点；根据年份确定各气象站点记录数据的开始时间和结束时间；根据各气象站点和各气象站点记录数据的开始时间和结束时间，按照第一列为维度、第二列为年份的方式依次将各气象站点的纬度、年份、温度和降水量保存于第一文件；将各可利用水量记录为一列存储于第二文件，各可利用水量与第一文件中的维度和年份相对应；根据第一文件和第二文件记录的数据，采用MATLAB计算各气象站点的帕尔默干旱指数。本发明提高了帕尔默干旱指数的计算效率。

说明书

技术领域

本发明涉及干旱指标计算领域，特别是涉及一种帕尔默干旱指数计算方法及系统。

背景技术

干旱指标是研究干旱问题的基础和关键，合理的干旱指标才能够较精准的监测和量化干旱，反映干旱的变化规律和趋势。国内外许多专家和学者致力于对干旱指标的研究，提出了很多干旱指标，并为提高其在不同区域的适用性进行了改进。Palmer(1965)综合考虑了降水量、土壤含水量、径流量和潜在蒸散量，基于水分平衡原理创建了帕尔默干旱指标(Palmer Drought Severity index，PDSI)，该指标物理机理较为清晰，PDSI能够很好的监测长期的发展和变化，因此被广泛用于干旱研究中。

PDSI在美国已经成为普适性的干旱指标，NOAA(National Oceanic andAtmospheric Administration，美国国家海洋和大气管理局)和USDA(United StatesDepartment of Agriculture，美国农业部)联合发布每周每月的PDSI预报图。美国农业部更是每天发布实时的干旱状态分布图，如果想了解感兴趣的区域，只需要点击该地区即可进行细分，像使用地图一样方便。本身而言，PDSI计算过程十分复杂，而后又根据不同地区的适用性情况进行了一系列的改进，衍生出各种不同版本的PDSI指标。

而目前计算该指标主要方法有两种，一种是美国内布拉斯加大学林肯分校提出的计算程序scpdsi.exe，该程序是基于C++编写，能够运用Windows系统进行单站点干旱指标运算，函数输入参数为4个或6个或7个，降水量和温度格式需要遵照要求输入。但该方法主要是针对于单个站点进行运算，对于中国甚至是更大范围的干旱监测与运算而言，单个站点逐一计算的效率过慢，且scpdsi.exe内部代码未知，无法根据不同地区不同下垫面情景进行计算过程的修正，不利于大范围干旱指标的计算和修改。另一类是基于Fortran语言的干旱指标计算，计算代码可以更改，但对于代码编程能力要求很高，需要拥有较强的代码编程基础以及数学功底，便于理清指标的计算过程和步骤，进而修正代码。

发明内容

本发明的目的是提供一种帕尔默干旱指数计算方法及系统，提高了帕尔默干旱指数的计算效率。

为实现上述目的，本发明提供了如下方案：

一种帕尔默干旱指数计算方法，所示方法包括：

采集多个气象站点的气象数据，所述气象数据包括纬度、年份、温度、降水量和可利用水量；

根据多个所述纬度确定各气象站点；

根据所述年份确定各所述气象站点记录数据的开始时间和结束时间；

根据各气象站点和各所述气象站点记录数据的开始时间和结束时间，按照第一列为维度、第二列为年份的方式依次将各气象站点的纬度、年份、温度和降水量保存于第一文件；将各可利用水量记录为一列存储于第二文件，各所述可利用水量与所述第一文件中的维度和年份相对应；

根据所述第一文件和所述第二文件记录的数据，采用MATLAB计算各所述气象站点的帕尔默干旱指数。

可选地，所述根据所述第一文件和所述第二文件记录的数据，采用MATLAB计算各所述气象站点的帕尔默干旱指数，具体包括：

根据各气象站点的所述气象数据计算各气象站点的水分平衡分量；

根据各气象站点的水分平衡分量计算各气象站点的水分亏缺量；

根据各气象站点的水分亏缺量计算各气象站点的水分亏缺指标；

根据各气象站点的水分亏缺指标计算各气象站点的帕尔默干旱指数。

可选地，所述根据各气象站点的所述气象数据计算各气象站点的水分平衡分量，具体包括：

采用Hamon算法或Thornthwaite法计算各气象站点的潜在蒸散发。

可选地，所述根据各气象站点的所述气象数据计算各气象站点的水分平衡分量，具体还包括：

根据各气象站点的降水量和潜在蒸散发，采用双层土壤模型计算各气象站点的实际蒸散发量、实际土壤水分补充量、实际产流量、实际土壤失水量、潜在土壤水分补充量、潜在产流量和潜在土壤失水量。

可选地，所述根据各气象站点的水分亏缺指标计算各气象站点的帕尔默干旱指数，具体包括：

根据各气象站点逐次干旱的历时时间和逐次干旱的累积水分亏缺指标进行拟合，获得拟合方程；

根据所述拟合方程获得所述帕尔默干旱指数的计算公式。

本发明还公开了一种帕尔默干旱指数计算系统，所示系统包括：

数据采集模块，用于采集多个气象站点的气象数据，所述气象数据包括纬度、年份、温度、降水量和可利用水量；

气象站点确定模块，用于根据多个所述纬度确定各气象站点；

时间确定模块，用于根据所述年份确定各所述气象站点记录数据的开始时间和结束时间；

文件记录模块，用于根据各气象站点和各所述气象站点记录数据的开始时间和结束时间，按照第一列为维度、第二列为年份的方式依次将各气象站点的纬度、年份、温度和降水量保存于第一文件；将各可利用水量记录为一列存储于第二文件，各所述可利用水量与所述第一文件中的维度和年份相对应；

帕尔默干旱指数计算模块，用于根据所述第一文件和所述第二文件记录的数据，采用MATLAB计算各所述气象站点的帕尔默干旱指数。

可选地，所述帕尔默干旱指数计算模块，具体包括：

水分平衡分量计算单元，用于根据各气象站点的所述气象数据计算各气象站点的水分平衡分量；

水分亏缺量计算单元，用于根据各气象站点的水分平衡分量计算各气象站点的水分亏缺量；

水分亏缺指标计算单元，用于根据各气象站点的水分亏缺量计算各气象站点的水分亏缺指标；

帕尔默干旱指数计算单元，用于根据各气象站点的水分亏缺指标计算各气象站点的帕尔默干旱指数。

可选地，所述水分平衡分量计算单元，具体包括：

潜在蒸散发计算子单元，用于采用Hamon算法或Thornthwaite法计算各气象站点的潜在蒸散发。

可选地，所述水分平衡分量计算单元，具体还包括：

水分平衡分量计算子单元，用于根据各气象站点的降水量和潜在蒸散发，采用双层土壤模型计算各气象站点的实际蒸散发量、实际土壤水分补充量、实际产流量、实际土壤失水量、潜在土壤水分补充量、潜在产流量和潜在土壤失水量。

可选地，所述帕尔默干旱指数计算单元，具体包括：

拟合方程获取子单元，用于根据各气象站点逐次干旱的历时时间和逐次干旱的累积水分亏缺指标进行拟合，获得拟合方程；

帕尔默干旱指数的计算公式确定子单元，用于根据所述拟合方程获得所述帕尔默干旱指数的计算公式。

根据本发明提供的具体实施例，本发明公开了以下技术效果：

本发明公开了一种帕尔默干旱指数计算方法及系统，按照第一列为维度、第二列为年份的方式依次将各气象站点的纬度、年份、温度和降水量保存于第一文件，将各可利用水量记录为一列存储于第二文件，各所述可利用水量与所述第一文件中的维度和年份相对应，根据所述第一文件和所述第二文件记录的数据，采用MATLAB计算各所述气象站点的帕尔默干旱指数，实现同时计算多个气象站点的帕尔默干旱指数，提高了帕尔默干旱指数的计算效率。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图1为本发明一种帕尔默干旱指数计算方法流程示意图；

图2为本发明一种帕尔默干旱指数计算系统结构示意图；

图3为本发明一种帕尔默干旱指数计算详细流程示意图；

图4为本发明MATLAB上GUI界面示意图；

图5为本发明帕尔默干旱指标的密度分布图、散点图及相关系数；

图6为本发明帕尔默干旱指数自动化计算系统解压后文件示意图；

图7为本发明设置Matlab的工作目录示意图；

图8为本发明AWC.txt文档示意图；

图9为本发明temperature and precipitation.txt文档示意图；

图10为本发明AWC数据预期数据格式预览示意图；

图11为本发明temperature and precipitation数据预期数据格式预览示意图；

图12为本发明运行成功界面；

图13为本发明Matlab当前文件夹示意图；

图14为本发明Palmer.txt结果示意图；

图15为本发明SPI、SPEI、SCPDSI及MPDSI的时间序列分布图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本发明的目的是提供一种帕尔默干旱指数计算方法及系统，提高了帕尔默干旱指数的计算效率。

为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。

图1为本发明一种帕尔默干旱指数计算方法流程示意图，如图1所示，一种帕尔默干旱指数计算方法包括以下步骤：

步骤101：采集多个气象站点的气象数据，所述气象数据包括纬度、年份、温度、降水量和可利用水量。

步骤102：根据多个所述纬度确定各气象站点。

步骤103：根据所述年份确定各所述气象站点记录数据的开始时间和结束时间。

步骤104：根据各气象站点和各所述气象站点记录数据的开始时间和结束时间，按照第一列为维度、第二列为年份的方式依次将各气象站点的纬度、年份、温度和降水量保存于第一文件；将各可利用水量记录为一列存储于第二文件，各所述可利用水量与所述第一文件中的维度和年份相对应。

所述根据所述第一文件和所述第二文件记录的数据，采用MATLAB计算各所述气象站点的帕尔默干旱指数，具体包括：

根据各气象站点的所述气象数据计算各气象站点的水分平衡分量；

根据各气象站点的水分平衡分量计算各气象站点的水分亏缺量；

根据各气象站点的水分亏缺量计算各气象站点的水分亏缺指标；

根据各气象站点的水分亏缺指标计算各气象站点的帕尔默干旱指数。

所述根据各气象站点的所述气象数据计算各气象站点的水分平衡分量，具体包括：

采用Hamon算法或Thornthwaite法计算各气象站点的潜在蒸散发。

根据各气象站点的降水量和潜在蒸散发，采用双层土壤模型计算各气象站点的实际蒸散发量、实际土壤水分补充量、实际产流量、实际土壤失水量、潜在土壤水分补充量、潜在产流量和潜在土壤失水量。

步骤105：根据所述第一文件和所述第二文件记录的数据，采用MATLAB计算各所述气象站点的帕尔默干旱指数，具体包括：

根据各气象站点逐次干旱的历时时间和逐次干旱的累积水分亏缺指标进行拟合，获得拟合方程；

根据所述拟合方程获得所述帕尔默干旱指数的计算公式。

下面以具体实施例说明本发明一种帕尔默干旱指数计算方法。

Step1：对输入数据进行整理

首先根据帕尔默干旱指数计算的要求，准备好输入数据，具体包括计算帕尔默干旱指数所需要的相应气象站点的的纬度、年份、温度、降水和可利用水量(Available WaterCapacity)。将纬度、年份、温度和降水数据存储于temperature and precipitation.txt文件(第一文件)中，规则为：第一列是气象站点的纬度(单位为度)，第二列是年份(格式为四位数)，第三列是温度(华氏摄氏度°F，°F＝℃*1.8+32)，第四个是降水(英寸，inch＝mm/25.4)，且要注意去掉表头。将可利用水量(Available Water Capacity，AWC)存为AWC.txt文件(第二文件)，注意可利用水量存储为一列，且要和temperature andprecipitation.txt文件中的对应气象站点对上，目的是为了后面Step2能够根据纬度和年份自动提取某个气象站点的数据，进行帕尔默干旱指数的自动化计算。

Step2：统计站点个数并分割数据

读取整理好的数据，首先根据气象站点的纬度信息统计参与PDSI计算的气象站点数，再根据气象站点的年份数据，记录每个气象站点的数据记录的开始和结束时间。由此对各个气象站点的所有气温和降水量进行分割。假设气象站点数为100，年份为1951-2010年，则第1-60行为气象站点1的所有运算数据，第61-120行为气象站点2的所有运算数据，以此类推。而后根据行列信息对可利用水量数据同样分割，以保证和气象数据的年份和气象站点匹配。

Step3：潜在蒸散发(PET)的计算

为使用者提供了两种计算方式，第一种是Hamon算法，第二种是Thornthwaite法，计算潜在蒸散发。

Hamon算法计算过程如下，

lat_rad＝(lat_list(j)*2*pi)/360

TLAT＝-tan(lat_rad)

PHI＝4093*sin(((2*(pi())/365)*J)-1.405)

w＝acos(TLAT*PHI)

D＝24*w/pi

DL＝D/12

PET_mm(k,i)＝(2.1*(DL(i)^2)*0.6108*exp((17.27*5*(T(k,i)-32)/9)/(237.3+5*(T(k,i)-32)/9))/(5*(T(k,i)-32)/9+273.3))*days_mo(i)。

其中，j表示是第j个站点；lat_rad表示是气象站点的纬度(弧度制)；lat_list表示气象站点的纬度(弧度制)，TLAT表示纬度(弧度制)的正切值，PHI为太阳偏斜角；w表示该站点的日落时角，D表示该气象站点的昼长(小时)，DL为该气象站点的昼长正态化之后的结果，PET_mm表示潜在蒸散发量(以毫米为单位)，k表示站点数和年份，即第k行；i表示月份(1,2,3,4....12)；days_mo表示12个月的日数，分别为31,28,31,30,31,30,31,31,30,31,30,31；pi表示圆周率。

Thornthwaite算法计算过程如下：

PHI＝0.4093*sin(((2*(pi())/365)*J)-1.405)

w＝acos(TLAT*PHI)

D＝24*w/pi

DL＝D/12

DM＝days_mo/30

h＝((T./5))^1.514

H＝sum(h,2)

a＝((6.75*10^-7)*H^3-(7.71*10^-5)*H^2+(1.792*10^-2)*H+0.49239)。

温度范围小于80华氏摄氏度时

PET_mm(k,i)＝16*DL(i)*DM(i)*(((10*T(k,i))/H(k,1))^a(k,1))。

温度范围大于80华氏摄氏度时，

PET_mm(k,i)＝(-415.85+32.24*T(k,i)-0.43*(T(k,i)^2))*DL(i)*DM(i)。

其中，DM表示对月时间长度正态化的结果，h表示k年第i个月的热指数，H表示所有年的温度的年热指数，a表示第k年的PET系数。

Step4：水分平衡分量的计算

ET、R、RO、L及PET、PR、PRO、PL统称为水分平衡分量，ET、R、RO、L分别代表实际蒸散发量(Evapotranspiration，ET)、实际土壤水分补充量(Recharge，R)、实际产流量(Runoff，RO)、实际土壤失水量(Loss，L)；PET、PR、PRO、PL分别代表潜在蒸散发(PotentialEvaporation，PET)、潜在土壤水分补充量(Potential Recharge，PR)、潜在产流量(Potential Runoff，PRO)、潜在土壤失水量(Potential Loss，PL)。其中，除PET根据Thornthwaite算法计算而得，其余水分平衡分量均利用实际降水量P和潜在蒸散量PET根据双层土壤模型计算而得，这也是基于步骤三的结果进行的。

双层土壤模型是指将土壤划分为两层，则土壤中的有效持水量(Available WaterCapacity，缩写为AWC)通常由两部分组成，上层土壤的有效持水量(AWC

AWC＝AWC

假设每月月初的上层和下层土壤初始含水量分别为Ss和Su，第一个月默认Ss＝AWC

PR＝AWC-(Ss+Su)。

土壤的潜在产量产流量PRO即为土壤总含水量：

PRO＝Ss+Su。

根据双层土壤模型，在降水不足以满足蒸散发的需求时，上层土壤水分可完全蒸发以补充水分亏缺，而下层土壤则遵循一定的规律丧失一部分水量。因此，潜在土壤失水量PL由两部分组成，如下所示：

PL＝PLs+PLu

PLs＝min(PE，Ss)

PLu＝(PE-PLs)·Su/AWC，当PLu≤Su。

实际的水分平衡分量ET、R、RO、L则根据降水量P和潜在蒸散量PET计算而得：

当PE≤P时，

ET＝PET

L＝0

R＝0

RO＝P-ET-PR。

当PE＞P时，

ET＝P+L

RO＝0

L＝L_s+L_u

R＝min(PET-P，PR)。

其中，

L_s＝min(PE-P，Ss)

L_u＝(PE-PLs)·Su/AWC，PLu≤Su。

然后基于以上步骤计算的各水分平衡分量，计算相应的水分平衡分量系数及气候适宜状态下各分量，具体步骤如下：

α

β

γ

δ

(ET)＝α

R＝β

(RO)＝γ

L＝δ

i表示一年中的不同月份，i＝1,2,3…12；α

Step5：气候适宜降水量和水分亏缺量的计算

基于Step4的结果计算气候适宜降水量P

P

并结合实际降水量计算出该月的水分亏缺量d。

d＝P-P

其中，d为水分亏缺量；P为实际降水量；P

Step6：水分亏缺指标Z的计算

水分亏缺量d被用于表征当前月份实际降水量与气候适宜状态下的降水量的差，但对于不同月份或者不同地区而言，其气候适宜状态下的降水量并不相同，这就导致了同样的d值在不同地区或不同月份的所反映的干湿程度并不一致(如在干旱区和湿润区或者在雨季和干季，同样的水分亏缺量所对应的干旱强度不同)。引进了一个修正系数即为K，修正系数K能够反应一个地区的水分供需关系。根据不同地区、不同月份将水分亏缺量d修正为水分亏缺指标Z。

水分需求量＝(PE)+(RO)+R；

水分供应量＝P+L；

K＝水分需求量/水分供应量；

Z＝dK。

Step7：回溯计算过程

为考虑前期水分盈缺的影响，对每次干旱的干旱历时及逐次干旱的累积Z值进行统计，绘制散点图并做拟合线，对拟合方程进行推导得到PDSI计算公式。

PDSI

当-0.5≤PDSI≤0.5时，判定处于正常状态，即一次干旱事件结束。因此可将问题转化为要得到多少水分才能使当前的干旱结束，回归正常状态。

-0.5＝0.897*PDSI

将PDSI

U_w＝Z+0.15。

同理，对于湿润状态的解除需要失去一定量的水分，定义一个有效增干量U_d，如下所示：

0.5＝0.897*PDSI

U_d＝Z-0.15。

使当前干(湿)期结束的概率P_e可以看做是得到(失去)的水分与确切的结束当前干(湿)期所需要得到(失去)的水分的百分比。

设定有效增湿(干)累积量V

若当前干期结束，

V

V

若当前湿期结束，

V

V

则当前干(湿)期结束的概率P_e可由下式得到

Pe＝V

当Pe达到100％时，则开始一个回算过程，以此来决定PDSI值。在前面计算PDSI值时，设置了三个指标来计算，分别为初始湿期的PDSI，初始干期的PDSI及当前确立干(湿)期的PDSI，实际的PDSI就是依据Pe从这三个中挑选出来的。

在干期或湿期交替时PDSI值不能当时确定，要在返回或该时期结束时回溯推求，时间一般较长。持续的干期或湿期时，其结束概率PPe＝0，当0

具体来看该帕尔默干旱指数(PDSI)自动化计算系统，指标计算过程大致如图3所示。

为整个系统制作了一个GUI界面，如图4所示，在MATLAB界面输入“guide”，然后创建一个GUI，包含两部分，一个是figure(图形)，一个是代码文件。GUI相当于是程序的启动器，这样即使用户不会运行代码或者修改代码，也可以做到在这个界面上，通过点击按钮进行干旱指标的计算。实现了操作界面设计简洁，操作方便简单，对操作者无编程语言上的要求。

具体使用的时候首先在MATLAB在界面上双机打开PDSI_Tool_Launcher.m，界面如图4所示。界面上一共拥有五个按钮，首先利用“输入AWC”和“输入降水和气温数据”将所需要的数据进行输入。根据站点顺序，分别编写temperature and precipitation.txt和AWC.txt，比如第1-50行是站点1近50年的降水等数据，对应的AWC.txt中的第一行是站点1的AWC数据。第51-100行是站点2近50年的降水等数据，对应的AWC.txt中的第二行是站点2的AWC数据。以此类推，将计算指标的原始数据放在一个文件中上传输入即可。这样可以一次计算多个站点的干旱指标，避免原先的计算方式中每计算一个站点需要运行一遍代码，实现多个站点同时计算，提高指标计算效率。

为PET的计算方法设定了两种，一个是Thornthhwaite算法，一个是Hamon算法，如果用户有需求可以通过增加一个.m文件，将算法编写下来，并在GUI界面上增加一个选项即可。具体的算法编写可以模仿Thornthwaite_PET.m和Hamon_PET.m。除此以外，自校正方式等代码也具有可更改性，方便使用者根据需求进行更正。

输入的数据通过选择的PET计算方法(Thornthhwaite算法或者Hamon算法)，调用对应的小程序(Thornthwaite_PET.m或Hamon_PET.m)计算出PET，即完成了图3的输入数据和计算PET。而后利用WaterBalance.m计算水量平衡系数，利用Zindex.m计算水分亏缺指数，根据Between0s.m的计算结果计算出现在是干期还是湿润期，分别调用Function_Ud.m和Function_Uw.m计算给定月份的结束概率PPe值，以判定是否回溯计算，根据BackTrack.m小程序进行回溯，最后通过Main.m给出三个备选项PX1、PX2、PX3，确定了PDSI和PHDI值将最终确定并在PDSI_Central.m函数中表化，输出结果完成计算。

本发明一种帕尔默干旱指数计算方法的技术效果如下。

1、与原有的计算方式相比，该计算可以同时进行多个站点的计算，大大提高了计算效率。除了在计算速率上保证之外，计算结果十分精确，通过与其他方法计算的指标进行对比，可以很好的反映干旱的时空变化过程。

2、计算过程简便，计算程序提供了一个界面，只需要点几个按钮即可完成多站点的干旱指标的计算过程，对代码编程能力要求低，就算是普通人也可以进行操作。

3、与原有计算方式来说，拥有开源的计算代码，可以根据研究区域的不同对指标进行相应改进，若有需求，可根据需要进行干旱指标的修正和代码的修改。

图15纵坐标为SPI(Standardized Precipitation Index，标准化降水指数)、SPEI(Standardized Precipitation-Evap-otranspiration Index，标准蒸发指数)、SCPDSI(Self-calibrating Palmer Drought Severity index，改进帕尔默干旱指数)和MPDSI(Modified Palmer Drought Severity Index，修正的帕尔默干旱指数)的指标值，横坐标为年份，图15显示的是中国1960年至2014年，共计55年的月尺度SPI、SPEI、SCPDSI和MPDSI干旱指标的时间序列曲线，其中0是干旱和湿润状态的分界线。总体而言，四种干旱指标的趋势大体一致，其中SPI和SPEI以及SCPDSI和MPDSI之间的一致性更强，这主要是由于两组指标在计算原理上差异较大造成，但总的来说干旱指标趋势大体一致。

为了更加直观的反映指标之间的相关性程度，利用2474个气象站点55年的数据进行探究，得出相关结果如图5所示。其中，图形矩阵的上三角为对应行和列所代表的干旱指标的相关系数，p＝0.001；图形矩阵的对角线为各指标值的密度分布图及边际地毯，对应的横纵坐标分别为上/下侧和右侧；图形矩阵的下三角为对应行和列所代表的干旱指标的分布位置，横纵坐标为左侧和上/下侧。

SPI和SPEI及SCPDSI和MPDSI的相关性分别为0.89和0.82，两组指标显著相关。从图形矩阵的下三角两组指标的散点分布情况可以看出，二者基本呈对角线分布。SPI与SCPDSI和MPDSI的相关性略差于SPEI，分别是0.33和0.46，但所有结果都通过了0.001的显著性检验。除此之外，根据密度分布图和边际地毯，发现指标存在明显的聚集现象，SPI和SPEI的指标值主要集中在-2到2之间，最大密度为0.4左右，而对于SCPDSI和MPDSI，指标主要集中在-4-4之间，最大密度为0.2左右。-2和-4恰好是对应SPI、SPEI及SCPDSI、MPDSI指标的极度干旱的分界线，这与指标原先的分级方式不谋而合，也就意味着，MPDSI按照原先的干旱指标分级标准是可行的，因此在MPDSI的指标等级划分上不做处理，保持与SCPDSI一致。

基于本发明一种帕尔默干旱指数计算方法采用MATLAB进行帕尔默干旱指数自动计算的步骤包括：

S1：准备计算干旱指标所需要的数据。

采集的数据应该被分成2个文本文件(.txt)。第一文件包含各个气象站点的AWC数据。第二文件应该包含纬度、年份、温度和降水等月值数据，用于计算PDSI。第一文件和第二文件都不应该有任何列标题。若有多个气象站点，则在两个文件中要按照站点的顺序对应排列。

S2：打开MATLAB，开启计算系统的操作界面，输入相关数据。

选择PET计算方法和矫正时间，而后帕尔默(Palmer)自动计算系统会自动利用气温数据计算蒸散发量PET和进行指标计算。

S3：等待运行完毕后，将所需的变量导出。

输出结果文件包含11个变量，分别是Latitude、Year、PET、Z-Index、PPe、PX1、PX2、PX3、X、PDSI、PHDI，选择所需的列将结果导出，常用的为PDSI指数。其中Latitude表示气象站点的纬度，Z-Index表示水分亏缺指标。

采用MATLAB进行帕尔默干旱指数自动计算的系统安装及运行基本要求如下：

硬件要求：

中央处理器CPU Intel Core i3 2.0G以上；

内存2G以上；

显示分辨率1024*768以上；

显示颜色24位彩色以上；

软件要求：

操作系统Microsoft Windows 7Profesional及以上；

Matlab 2012及以上，案例示范中选取的版本为Matlab R2016a。

软件安装完毕后，将压缩文件(帕尔默干旱指数(PDSI)自动化计算系统.zip)进行解压缩，解压后如下图6所示，包含13个代码文件和2个测试数据(temperature andprecipitation.txt和AWC.txt)。

运行Matlab并将Matlab的工作目录设置为相应的帕尔默干旱指数(PDSI)自动化计算系统解压后的存储路径，设置完成后即可在Matlab左侧当前文件夹中出现所有的代码文件和附加资料。如下图7所示：

在运行系统之前，必须下载适当的文件，并将数据按照一定的格式整理完毕。在本节中，列出了所需要的数据，并说明了正确的数据输入格式。

数据应该被分成2个文本文件(表1)。第一个文件包含各个气象站点的AWC，如图8所示。第二个文件应该包含纬度、年份、温度和降水，如图9所示，为月值数据，用于计算PDSI。文件都不应该有任何列标题。若有多个站点，则在文件中要按照站点的顺序排列，排列要一致。

表1数据类型及参数

上面是示例数据文件，需要注意的是，由于中国区域气象站点众多，为保证不同文件的站点排列顺序一致，不仅需要比较其纬度，还需要比较其经度数据，一般将经纬度相同的视为同一站点，此步骤为数据预处理，可自行选择软件完成，不多赘述。

图8反映的是纬度为53.466999°的站点的AWC值；图9对应的是该站点1960-2014年间的逐月的纬度、年份、温度和降水数据。注意，每个条目都有年份和纬度，整个文件按照站点顺序排列，在某一站点内部则按照年月排列。

本节提供了运行该系统的详细步骤及相关细节、帮助菜单上的详细信息以及对用户可选择的不同选项的解释。

Matlab可以对工作文件夹和文件位置非常严格。为了使系统顺利运行，请遵循下面的Matlab说明在运行系统之前，确保所有的脚本(包括.fig文件)都必须位于同一个文件夹中，并将该文件夹设为工作文件夹，以便系统能够正常运行。未压缩的文件夹可以用作工作文件夹，或者脚本可以移动到另一个合适的工作文件夹，只要确保使用的文件夹包含所有下载的脚本.m和.fig文件。

在运行该工具时，包含脚本的文件夹必须是出现在Matlab屏幕左侧的“当前文件夹”(图7)。

数据文件(.txt)可以存储在任何文件夹中，而不必存储在工作文件夹中。数据导入工具允许在导入数据时导航到不同的文件夹。但是，如果用户怕遇到问题，可以将数据移动到工作文件夹，以避免任何潜在的报错。

操作步骤：

(1)安装软件，下载压缩文件后解压，并将所有Matlab脚本和.fig文件放在同一个文件夹中。

(2)打开Matlab，将放置Matlab脚本的文件夹设置为当前文件夹。

(3)如果由于某种原因用户无法在Matlab中设置工作目录，则可以在代码所在文件夹通过右键单击PDSI_Tool_Launcher.m，选择“用Matlab打开”。

(4)Matlab工作目录设置完毕之后，从左侧的窗口中双击打PDSI_Tool_Launcher.m。

(5)点击运行按钮”运行”，系统操作界面会弹出(图4)。

(6)点击“输入AWC”按钮；

(7)导航到用户存储其AWC数据文件的文件夹。

(8)选择AWC数据文件并单击打开。

(9)预览屏幕应该打开，检查数据格式，如果数据正确，按“下一步”。

(10)在下一个屏幕上按“完成”，不要点击任何其他按钮。

(11)此时在命令行窗口，应该出现一个确认AWC数据已成功加载的信息，“AWC数据已成功输入”。

(12)点击“加载气温和降水数据”按钮；

(13)对温度和降水数据重复步骤(7)–(10)，预览信息如图11所示。

(14)正确加载的数据在完成时应该会出现提示信息“气温和降水数据已成功输入”；

(15)如果在加载数据时出错，则从步骤(4)重新开始。

(16)选择所需的PET计算方法，包括Thornthwaite算法和Hamon算法。默认是Thornthwaite算法。

(17)选择所需的校准周期。默认选项是NOAA(1931-1990)。根据计算需要选择全部数据。

(18)点击运行按钮，当系统完成PDSI的计算时，会弹出提示消息：

运行成功！输出文档-"Palmer.txt"-，如图12。

保存位置为C:\Users\Administrator\Desktop\专利软著\专利和软件著作权---yu\帕尔默干旱指数(PDSI)自动化计算系统。

该运行系统有一个内置的帮助菜单，要访问帮助菜单，请在GUI中单击“帮助”按钮。若要退出“帮助”菜单，请再次单击“帮助”按钮。

选项说明

PET的算法

该运行系统为用户提供了提供了两种PET计算方法:Thornthwaite和Hamon。在PDSI最初提出时，Palmer使用了Thornthwaite算法计算潜在蒸散发，随着科学不断探索，提出了更加全面和准确的Hamon算法和Penman算法，但同时它对数据要求较多，而Thornthwaite对数据要求不多。

校准周期的选择

校准周期是用于计算一个地点的气候适宜态降水及修正因子等。NOAA使用1931年1月至1990年12月期间作为其校准周期(见Karl,1986)，但当在没有长数据记录或使用更全面的时间间隔时，还提供使用完整记录作为校准周期的选项。

结果输出

一旦系统完成运行，所有运行结果都会保存在一个输出文件Palmer.txt，将出现在Matlab的当前文件夹中(如图13)，如果多次运行该工具，请确保移动或更改输出文件的名称，因为在多次运行期间原输出文件会被覆盖。

输出文件包含11个变量，分别是Latitude、Year、PET、Z-Index、PPe、PX1、PX2、PX3、X、PDSI、PHDI(图14)；输出变量的详细信息可以在Palmer(1965)和Alley(1984)中找到，PET的输出单位是英寸(inch)。

图2为本发明一种帕尔默干旱指数计算系统结构示意图，如图2所示，一种帕尔默干旱指数计算系统包括：

数据采集模块201，用于采集多个气象站点的气象数据，所述气象数据包括纬度、年份、温度、降水量和可利用水量；

气象站点确定模块202，用于根据多个所述纬度确定各气象站点；

时间确定模块203，用于根据所述年份确定各所述气象站点记录数据的开始时间和结束时间；

文件记录模块204，用于根据各气象站点和各所述气象站点记录数据的开始时间和结束时间，按照第一列为维度、第二列为年份的方式依次将各气象站点的纬度、年份、温度和降水量保存于第一文件；将各可利用水量记录为一列存储于第二文件，各所述可利用水量与所述第一文件中的维度和年份相对应；

帕尔默干旱指数计算模块205，用于根据所述第一文件和所述第二文件记录的数据，采用MATLAB计算各所述气象站点的帕尔默干旱指数。

所述帕尔默干旱指数计算模块205，具体包括：

水分平衡分量计算单元，用于根据各气象站点的所述气象数据计算各气象站点的水分平衡分量。

水分亏缺量计算单元，用于根据各气象站点的水分平衡分量计算各气象站点的水分亏缺量。

水分亏缺指标计算单元，用于根据各气象站点的水分亏缺量计算各气象站点的水分亏缺指标。

帕尔默干旱指数计算单元，用于根据各气象站点的水分亏缺指标计算各气象站点的帕尔默干旱指数。

所述水分平衡分量计算单元，具体包括：

潜在蒸散发计算子单元，用于采用Hamon算法或Thornthwaite法计算各气象站点的潜在蒸散发。

可选地，所述水分平衡分量计算单元，具体还包括：

水分平衡分量计算子单元，用于根据各气象站点的降水量和潜在蒸散发，采用双层土壤模型计算各气象站点的实际蒸散发量、实际土壤水分补充量、实际产流量、实际土壤失水量、潜在土壤水分补充量、潜在产流量和潜在土壤失水量。

所述帕尔默干旱指数计算单元，具体包括：

拟合方程获取子单元，用于根据各气象站点逐次干旱的历时时间和逐次干旱的累积水分亏缺指标进行拟合，获得拟合方程；

帕尔默干旱指数的计算公式确定子单元，用于根据所述拟合方程获得所述帕尔默干旱指数的计算公式。

本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。对于实施例公开的系统而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。

本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上所述，本说明书内容不应理解为对本发明的限制。