云南省16个州市降水量特征分布研究——基于克里格插值法

Study on the Characteristic Distribution of Precipitation in 16 Cities in Yunnan Province—Based on Kriging Interpolation Method

摘要: 云南省是我国山地灾害最为严重的区域之一，掌握该地区16个州市降水的变化规律是山地灾害成灾演化的关键控制因素。本文使用实验变差函数克里格插值法，拟合相适应的球状模型，并计算出每个各点周围降水数据的权重。将得到的权重与已知的降水数据相乘得到相应格点的降水量，并使用Matlab软件画出降水分布图，以此模拟出云南省16个州市的降水量的降水图。研究发现：月降水量多雨区主要分布在省内腹地，少雨区主要分布在省内西北部，从南向北逐渐减少。降水量的空间分布预测对于云南省区域水资源环境及气候变化研究展具有重要的意义。

Abstract: Yunnan Province is one of the regions with the most serious mountain disasters. The key control factor for the evolution of mountain disasters is to master the variation rule of precipitation in 16 cities and towns in this region. In this paper, the experimental variation function Kriging interpo-lation method is used to fit the suitable spherical model and calculate the weight of precipitation data around each point. The obtained weight is multiplied by the known precipitation data to the precipitation of the corresponding grid points, and the Matlab software is used to draw the pre-cipitation distribution map, so as to simulate the precipitation map of 16 cities in Yunnan Province. The results show that the areas with heavy monthly rainfall are mainly distributed in the hinter-land of the province, while the areas with little rainfall are mainly distributed in the northwest of the province, decreasing gradually from south to north. The prediction of the spatial distribution of precipitation is of great significance to the study of regional water resources environment and climate change in Yunnan Province.

1. 引言

全球气候变化导致全球范围内极端天气气候事件频繁发生，且发展趋势不断增加，IPCC第六次报告指出全球变暖将会导致强降水事件发生的频、强度和量级增加 [1] (黄雪松等，2008)；中国极端气候、天气事件发生的概率明显增多，西南区域性极端降水事件发生频率显著上升，南方地区单日极端降水事件显著增加 [2] (韦志刚等，2021)。极端降水事件引发的气象灾害(洪涝、山体滑坡、泥石流等)频发，对人民生命财产、社会经济和生态环境的可持续发展造成的巨大危害备受关注。近年来，国内外关于极端小时降水事件、极端小时降水阈值的确定及其时空分布特征的研究越来越多。研究表明，由于受下垫面、海陆位置、季风、太阳辐射、山谷热力环流等因素影响，极端小时降水呈现出明显的季节性、区域性、不对称性和局地性，地形地貌条件比较复杂的区域极端小时降水总量和日变化具有明显的区域性，不对称性更为突出 [3] [4] (潘欣等，2019；宇如聪等，2014)。目前有关极端降水事件的研究使用的资料多为24 h累积降水量，由于降水事件在时空分布上的不对称性和不均匀性，24 h累积降水有可能高估或低估由持续性弱降水和短时强降水事件引发的极端降水强度(江远安等，2017) [5] ，而逐小时降水或更高分辨率的降水数据可避免24 h累积降水带来的不确定性，在资料筛选中减少误差，能更准确地表征实际降水强度并详细地阐释降水过程特征(杨霞等，2020) [6] 。

在研究降水规律方面，Fan等(2013)基于太湖流域中96个雨量站的年最大降水，利用L矩阵方法，发现极端降水逐渐集中在城市区域 [7] ；董满宇等(2017)等基于1960~2017年逐日降水数据，运用R/S分析等方法揭示了年降水总量中，中雨量占比例最大，为32.05% [8] ；包云轩等(2018)采用EOF正交经验法、线性倾向率法研究环太湖地区近40年来的气候变化，发现区域降水表现为北部增加，南部减少 [9] ；杨凤等(2019)利用MK趋势检验分析流域及其水利分区的降水变化趋势，发现其降水指标均呈上升趋势，且东部地区上升趋势更明显 [10] 。

2. 研究方法与数据来源

2.1. 研究方法

2.1.1. 实验变差函数

变差函数是描述随机场和随机过程空间相关性的统计量，被定义为空间内两空间点之差的方差。在实际应用中，把区域化变量 $Z\left(x\right)$ 在空间点x， $x+h$ (h代表滞后距)两点处的值之差的方差之半定义为 $Z\left(x\right)$ 在x方向上的变差函数，记为 $\gamma \left(x,h\right)$ ，即

$\gamma \left(x,h\right)=\frac{1}{2}Var\left[Z\left(x\right)-Z\left(x+h\right)\right]$ (1)

在满足二阶平稳假设与本征假设的情况下， $\gamma \left(x,h\right)$ 与x取值无关，只依赖于h。于是得到：

${\gamma }^{*}\left(h\right)=\frac{1}{2N\left(h\right)}{\displaystyle {\sum }_{i=1}^{N\left(h\right)}{\left[Z\left(x_i\right)-Z\left(x_i+h\right)\right]}^2}$ (2)

其次，将滞后距h间隔取为0.2，结合数据稳定性共取26个数值，同时求得变差函数值 ${\gamma }^{\text{*}}=\left(h_i\right)$ 。

2.1.2. 克里格插值法

克里格法也称为克里金(Kriging)，是地址统计学的主要内容之一，从统计意义上说，是从变量相关性和变异性出发，在有限区域内对区域化变量的取值进行无偏、最优估计的一种方法；从插值角度讲是对空间分布的数据求线性最优、无偏内插估计一种方法(白瑞和，2010) [11] 。克里格法的适用条件是区域化变量存在空间相关性。

设 $Z\left(x\right)$ 是二阶平稳的点承载物区域变化量，在这里对以 $x_0$ 为中心的盘区平均值

$Z_V\left(X_0\right)=\frac{1}{V}{\displaystyle {\int }_{V}Z\left(X\right)\text{d}x}$ 进行估计测算。在克里格插值法中，用待估点变量周围影响范围内n个已知变量

值的线性组合来求其估计值。

定义公式为

$Z\left(x_0\right)={\displaystyle {\sum }_{i=1}^n{\lambda }_{i}Z\left(x_i\right)}$ (3)

其中， $Z\left(x_i\right),i=1,2,\cdots ,n$ 表示已知的观测值，它们分别代表区域内已知的观测值，待估点变量 $x_0$ 是一个未采样点，需要用克里格法进行最优、线性、无偏估计。在这里为了避免系统出现误差，使 $X\left(x_0\right)$ 的估计无偏，必须要选取适合的 ${\lambda }_i$ ，使得 ${\displaystyle {\sum }_{n=1}^n{\lambda }_i}=1$ ，根据未采样点需满足线性无偏、最有估计条件，用拉格朗日乘数法得到普通克里格方程组即上述等式(3)可变成为：

$\{\begin{array}{l}{\displaystyle {\sum }_{i=1}^n{\lambda }_i\gamma \left(x_i,x_j\right)}+u=\gamma \left(x_i,x_0\right)\\ {\displaystyle {\sum }_{i=1}^n{\lambda }_i}=1\end{array}$ (4)

其中，u是极小化处理时所求出的拉格朗日乘数， $\gamma \left(x_i,x_j\right)$ 是随机变量Z在采样点 $x_i$ 和 $x_j$ 之间的半方差值，是求得的最优拟合实验变差函数值。

2.2. 数据来源

本文以云南省2012年7月126个气象站点降水量为样本。该组数据的均值是13.6 mm，中位数为14 mm，从均值降水量均值和中位数差值来看，2012年7月云南省降水量波动较小，分布较为均匀。2012年7月云南省最大降水量为22.4 mm，最小降水量为8.2 mm，极差为14.2 mm，说明2012年7月云南省126个气象站点降水量的极值间有着较大差距。表1是该组样本数据的描述。

Table 1. Basic statistical information of precipitation (Unit: mm)

表1. 降水量的基本统计信息(单位：mm)

运用观测到的数据画出2012年7月云南省的降水量分布图，如图1所示。图1中x轴表示经度，y轴表示纬度。降水量的分布在图中表示了整个云南省16个州市的地理位置分布，也即云南省的地图。从图中可以看出，降水量集中分布在云南省的中东部，西部降水量较少；且北部降水量略少于南部，降水量分布不均。其中，每个格点的数据由周围已知的降水量数据计算得到，且为了结果的稳健性，本文采取格点周围不同象限的降水量来测度格点的降水量。

Figure 1. Precipitation distribution map of Yunnan Province

图1. 云南省降水量分布图

3. 结果分析

本文主要选取滞后距h = 0.2、0.4、0.6等26个距离，用excel处理数据，根据等式(2)可得到实验变差函数值 ${\gamma }^{\text{*}}\left(h\right)$ ，表2可知：

Table 2. Experimental variogram values

表2. 实验变差函数值

为了估计区域化变量的未知值，再用球状模型将实验变差函数拟合成理论变差函数，球状模型函数表达式为

$\gamma \left(h\right)=\{\begin{array}{l}0,h=0\\ C_0+C\left(\frac{3}{2}\frac{h}{a}-\frac{1}{2}\frac{h^3}{a^3}\right),0<h\le a\\ C_0+C,h>a\end{array}$ (5)

式中，C₀为块金常数， $C_0+C$ 为基台值，C为拱高，a为变程， $\gamma \left(h\right)$ 为拟合降水量，h为滞后距。

经过计算，其中 $a=3.89$ ， $C_0=1.41$ ， $C=8.16$ 。

本文拟合的球状模型为：

$\gamma \left(h\right)=\{\begin{array}{l}1.41,h=0\\ 1.41+8.16\left(\frac{3}{10}h-\frac{1}{250}{h}^3\right),0<h\le 3.89\\ 8.16,h>3.89\end{array}$ (6)

依据表2的实验变差函数值得到的散点图与上述理论球状模型函数图拟合，得到拟合图如图2所示：

Figure 2. Spherical model fitting

图2. 球状模型拟合

云南省的经度大约在97˚~107˚，纬度大约在20˚~30˚，取经度和纬度1度为单位画线，以线的交点为本论文的一个待估点，在每一个待估点附近的不同方向上选择距离待估点最近的三个样本点进行预测。设距离待估点S₀最近的三个样本点分别为S₁，S₂，S₃，其对应的三个样本点的降水量分别为y₁，y₂，y₃，则待估点S₀的估计降水量为 $y_0^{\text{*}}={\displaystyle {\sum }_{i=1}^3{\lambda }_i{y}_i}$ 。则列出普通克里格的矩阵形式为：

$\left[K\right]\left[\lambda \right]=\left[M\right]$

其中：

$\left[\lambda \right]=\left[\begin{array}{c}\begin{array}{c}{\lambda }_1\\ {\lambda }_2\\ {\lambda }_3\end{array}\\ {\lambda }_4\end{array}\right]$ ， $\left[K\right]=\left[\begin{array}{cccc}{C}_{11}& C_{12}& C_{13}& C_{14}\\ C_{21}& C_{22}& C_{23}& C_{24}\\ C_{31}& C_{32}& C_{33}& C_{34}\\ 1& 1& 1& 0\end{array}\right]$ ， $\left[M\right]=\left[\begin{array}{c}\begin{array}{c}{C}_{01}\\ C_{02}\\ C_{03}\end{array}\\ C_{04}\end{array}\right]$ (7)

将云南省2012年7月降水量观测点降水量图分成1˚ × 1˚的经纬度网格图。依次利用球状模型表达式及普通克立格方程组及方差的矩阵模拟出每个交叉点出的降水量，最后得出云南省2012年7月份降水量如表3所示。

Table 3. Spherical simulated precipitation value table

表3. 球状模拟降水量值表

利用Matlab软件，依据球状模型模拟的降水量数据得出云南省降水量空间分布图，如图3。由图3可知云南省降水量空间分布图，可以知道云南省从南到北降水量由充沛转为稀少，省内降水量差异明显。针对这种差异，可以因地制宜地安排适宜的农业生产和加强水利等基础建设。如在降水量丰富的滇中、滇西南区域种植喜水的经济作物，合理安排农业生产；在降水量稀少的滇西北区域种植耐旱的经济作物，同时注意兴修水利设施，保障当地居民基本的用水需求。近70年来，由于全球气候变暖云南省气温总体也呈现显著上升趋势，致使大部分地区相对湿度有所下降，降雨量略有减少，但是部分地区趋势明显，同时大部分地区日照时数呈增加趋势。众所众知，气候变化对参考作物需水量影响区域差异明显(乔伟等，2021) [12] 。云南省省内主要表现为滇中、滇东、滇东南、滇东北、滇南的需水量呈减小趋势，滇西、滇西北、滇西南的需水量呈现增大趋势。结合云南省历年统计年鉴可以发现，2006~2016年，省内农作物播种面积逐年增加，种植结构变化较大。其中，谷类种植面积占比降低5%，蔬菜种植面积占比增加了5%，同时基础粮食种植面积与经济作物种植面积比呈现减小的趋势。同时，随着经济发展以及对于水利基础设施的兴建，云南省节水灌溉面积显著增加，增加幅度达到70%，农田亩均灌溉用水量减少26.7%。

Figure 3. Spherical model simulation diagram

图3. 球状模型模拟图

总的来说，在种植面积、种植结构和节水灌溉等人类活动的综合作用下，云南省农业用水量基本保持稳定。但是，未来形势依然严峻，因为在气候变化、种植结构、灌溉管理等变化环境下农业需水规律发生了深刻变化，从而影响农业用水格局。云南省大部分地区由于气候变化、作物种植面积增加，导致农业需水量呈增加趋势，农业用水供需矛盾加剧，且随着未来环境变化不断恶化，而通过农田水利工程建设提高农业用水效率的潜力有限。因此，今后应加强变化环境对农业用水及粮食生产安全的影响研究，为制定保障云南省农业用水安全的适应性对策提供理论依据。

4. 结论与启示

气象学家指出：近年来中国气温变化的总体特征以增温为主，增温增湿的地区面积最为广大，季节变化构型也以增温增湿为主，秋季略有异常，与此同时，冰川融化增强，降水量与河道径流增加。本文发现，造成云南省降水分布的因素有地形因素：① 主要与云南西北特殊地形有关，滇西北紧挨青藏高原，位于其东南侧，处于高原与中低纬季风地区。气候特点复杂多样，首先在晴空条件下，无雪盖的高山白天太阳直接辐射强度和夜间有效辐射强度随高度增加而增大。因坡向不同，阳坡和阴坡得到的太阳辐射不同，并因此影响气温和气流的分布。② 在一定高度以上的山地，由于气流中水汽含量减少，降水量又随高度增加而减少。云南西北一角的降水量在近年气候变化的大背景下多与西部高原地区的气候变化保持一致。

此外，本文认为气候因素也是影响云南省降水的重要原因：① 云南省所处的地区是热带季风气候，云南中部和西南地区的月降水量主要受热带季风气候影响。云南省位于中国西南边疆，地属亚热带低纬高原季风型气候，山地立体气候显著，降水适中干湿分明，但存在降水时空分布不平衡，降水变异增强和极端降水频发等现象，主要受青藏高原和东亚、南亚两支季风影响显著。② 近年来随着经济和科技的发展，全球和地区人类经济发展活动引起的全球气候变化也是影响该区域气候变化的重要因素之一，以西双版纳、普洱、临沧为代表的云南南部和西南部的北热带、南亚热带雨林属于高温、静风地带。滇西南地区较大的月降水量主要与上述几个原因有关。