王磊之，胡庆芳，王银堂，胡 艳，刘 勇

(1.南京水利科学研究院水文水资源与水利工程科学国家重点实验室，江苏南京 210029; 2.太湖流域管理局水文局(信息中心)，上海 220000)



基于广义可加模型的降水空间估算模型

王磊之1，胡庆芳1，王银堂1，胡 艳2，刘 勇1

(1.南京水利科学研究院水文水资源与水利工程科学国家重点实验室，江苏南京 210029; 2.太湖流域管理局水文局(信息中心)，上海 220000)

采用太湖流域浙西区33个站点1989—2013年的降水观测资料，基于广义可加模型建立了该区域年降水空间估算模型，得到了该区域分辨率为1 000 m×1 000 m的降水栅格数据。借助多种统计指标，分析了不同预测因子对降水空间估算结果的影响，特别比较了是否将高程作为预测因子情况下的年降水量估算差异。结果表明：加入高程因子作为广义可加模型的预测变量后，无论是模型的优良性还是降水空间估算精度均有所提高，能够更合理地描述浙西区降水空间分布受地形影响的特征。此外，从模型偏态效应图来看，浙西区降水量随着纬度的增加明显减小，随着高程的增加显著增大，而随经度的变化则不明显。

降水空间估计； 年降水量； 广义可加模型； 太湖流域

降水是重要的水文气象变量之一，精确获取降水空间分布对于气象、水文、生态及地质等领域均具有重要的理论和实际意义[1]。但是，受大气运动、海陆位置和下垫面状况(地形等)等多种因素的影响，降水的空间分布十分复杂，是最难估算的水文气象要素之一。

由于降水空间定量估算在众多科学领域和生产实践中的基础性作用，长期以来国内外针对这一问题进行了大量研究。到目前为止，降水空间估算方法有数十种之多。从降水空间估计所采用的信息源来看，传统的降水空间估算方法可分为两种。一种是仅采用地面雨量站网所提供的降水空间自相关性信息进行估算，这类方法包括Thiessen多边形、反距离权重、普通克里金等方法；另一种在降水空间自相关性信息的基础上，进一步利用降水与地形地貌因子之间的空间互相关性进行估算。这两个方面的文献均较多。如钱永兰等[2]在月和日时间尺度上，比较了薄板样条函数(Thin plate spline，TPS)、反距离权重(Inverse distance weighting，IDW)和普通克里金法(OK) 估算区域降水中的差异，指出 TPS 的精度要高于后两种算法；舒守娟等[3]建立了以经纬度、高程、坡度、坡向为预测变量的年、季降水回归模型，估算了青藏高原的降水空间分布；石朋等[4-6]均发现加入高程信息的协克里金法(Co-kriging,CK)对降水的估算精度要高于OK法。同时，除了传统的降水空间估算法，近年来，一些学者又从多源降水信息融合的角度来研究降水空间估算，如胡庆芳[7]根据卫星降水信息提供的背景场与地面降水观测进行融合，建立了3种降水融合模型；潘旸等[8]采用贝叶斯融合法实现了雷达估测降水、卫星反演降水与地面站点观测降水3种资料的融合。

在众多空间估算方法中，广义可加模型(Generized Addictive Model，GAM)是一种基于非参数回归和平滑技术发展起来的方法，其优点是：没有参数方法对先验分布的限制条件；能灵活地将参数和非参数函数加入到模型中[9]，便于揭示不同预测变量与响应变量的非线性影响。目前，该方法在降水等气象水文要素的空间估算中已得到一定应用，如J.Aalto等[10]采用GAM和外部漂移的克里金法(KED)分别对芬兰1981—2010年的月均气温和月均降水进行插值，发现GAM的月均气温插值精度比KED更高，但在月均降水插值精度上较KED略低；C.Huang等[11]采用GAM中的薄板平顺样条函数融合遥感降水数据和地面雨量站网数据，估算了新西兰的日降雨空间分布，其估算结果精度较高，但研究中未将高程作为预测变量。

为深化对广义可加模型在降水空间估算中的适用性和可靠性的认识，本文以太湖流域浙西区为研究区域，采用区域及周边33个雨量站连续25年(1989—2013年)的降水观测数据，开展广义可加模型估算年降水量空间分布的研究。针对前人在GAM降水空间估算中未引入高程作为预测变量，在此研究中引入高程因子，并采用多种定量精度指标，基于连续25年的模型拟合结果，探讨GAM考虑不同预测因子对年降水空间估算结果的影响，特别是针对是否将高程作为预测变量情况下降水估算结果的差异加以分析，其研究结果可为地形条件复杂的山丘区降水空间估算、水文分析计算及水资源评价提供科学依据。

1 研究方法

1.1 模型基本原理

GAM是广义线性模型(GLM)的半参数扩展，由T.J.Hastie和R.J.Tibshirani在1990年提出[12]。与广义线性模型相似，GAM包括1个随机成分，1个可加成分，以及1个连接函数(用于连接这两个成分)。该模型假设函数是相加的，函数的组成成分是非参数的光滑函数：

g{μ(X)}=α+f1(X1)+f2(X2)+…+fp(Xp)+ε

(1)

式中：α+f1(X1)+…+fp(Xp)为可加成分;g(μ)为连接函数，μ=E(Y/X1，X2，X3，…，Xp);α为截距;ε为残差变量;fi(Xi)(i=1,2,…)是针对每个解释变量的任意单变量光滑函数，且各项没有具体的参数形式，响应变量的分布可以是指数分布族中的任意一员，如二项分布、Poission分布、Gamma分布等，即：

(2)

式中：θ为自然参数;φ为尺度参数。

与广义线性模型相比，广义可加模型更加强调对数据进行分参数的探索，因而大大降低了线性设定带来的模型设定风险[13]。由于对样本的要求较少，使得广义可加模型的应用较为广泛。

1.2 模型建立方法

以雨量站的经度、纬度和高程为预测变量，以年降水量为响应变量，建立GAM模型，并按照有无高程因子作为预测变量，模型1采用各站点的经向、纬向坐标作为预测变量，模型2在经纬向坐标的基础上增加高程作为预测变量。两个模型的数学形式为：

P1=α0+f1(x)+f2(y)+ε

(3)

P2=α0+f1(x)+f2(y)+f3(h)+ε

(4)

式中：P1和P2为各站年降水量;α0为截距；f1,f2,f3分别为经向、纬向坐标和高程对应的非参数光滑函数；ε为残差变量。

按照式(3)和(4)，对每一年分别建立GAM估算模型，建模过程如下：分别采用各站点年降水数据及经向、纬向坐标和高程值拟合式(3)和(4)，估算出样条函数f1，f2，f3的系数和β0，ε，该建模过程通过R语言环境下“mgcv”包[14]中的“gam”函数实现；当参数估算完成后，即可计算各空间位置对应的降水量，从而得到降水空间估算结果。计算过程通过“mgcv”包中的“predict”函数完成。

本次建模中，“mgcv”包为预测变量提供了不同的平滑函数[15]以及它们的连接方式，平滑函数包括平滑样条、张量积平滑和张量积相互作用函数，本次拟合采用最为常用的平滑样条函数。平滑样条函数是分段多项式拟合函数，常用极小化惩罚平方和来求解：

(5)

平滑函数节点间的连接方式包括薄板光滑样条、三次回归样条、P-样条和自适应样条。因此，GAM比通常采用的薄板光滑样条函数更具有普适性。根据交叉验证的精度统计结果，选用拟合精度最高的3次回归样条来表示各预测变量中的平滑项f1，f2，f3，其表达式为三次多项式：

(6)

1.3 模型性能评价

采用AIC准则[16]以及修正后的决定系数(AdjustedR2，修正R2)来比较不同模型拟合的优良性。AIC准则同时考虑了模型的统计拟合度以及参数数量的影响。MAIC值越小，模型性能越好；对决定系数进行修正是为了体现GAM中预测变量数量不同对估算结果的影响，修正后决定系数较大的模型应该优先选择。

(7)

(8)

(9)

根据留一法交叉验证[17]得到的4项精度指标反映所建模型的计算精度及预测能力，分别为各站点实测值与拟合值的平均绝对误差、平均相对误差、均方根误差以及线性相关系数。

2 研究区域与数据

浙西区是太湖流域8个一级水利分区之一，位于太湖流域西南(如图1(a)所示)，面积为5 931 km2，占太湖流域面积的16.8%。浙西区的地形条件比较复杂，地面高程在2～1 578 m之间，山地和丘陵面积为4 608 km2，约占整个区域面积的78%。浙西区各站多年平均降水量为1 022～1 879 mm，是太湖流域降雨量最高的地区，加之该区位于太湖流域上游，使得该区成为太湖水量的主要来源地区之一。据统计，该地区入湖水量约占太湖总入湖水量的50%[18]。受地形影响，浙西区降雨空间分布差异较大。

目前，浙西区已经建立了比较密集的雨量站网，且具备较长系列的降雨观测资料，可供建立降水空间估算模型。根据水利部水文局刊印的水文年鉴及太湖流域管理局提供的资料，整理了1989—2013年(共25年)浙西区33个雨量站的年降水量，雨量站空间分布如图1(b)所示。同时，获取了这33个雨量站点的Lambert投影坐标，并采用DEM提取了这些站点的高程。

3 结果分析与讨论

3.1 年降水估算模型精度

按照是否将高程作为预测因子，针对各年降水量建立了两个GAM模型。其中，模型1的预测变量为经纬向坐标，模型2为经纬向坐标和高程。表1给出了1989—2013年两个模型精度指标的统计参数。由该表可知：与模型1相比，模型2各年相关系数的均值、最大值、最小值均明显提高，各年平均绝对误差，平均相对误差和均方根误差的均值、最小值、最大值均明显下降。同时，模型2各项精度指标的标准差均低于模型1。因此，相比于模型1，模型2的估算精度更高、更稳定。这说明，对于浙西区而言，高程对降水量的空间分布具有重要影响，在降水空间估算过程中对高程因子的作用应当予以考虑。

表1 两个模型拟合结果定量精度指标统计(1989—2013年)

Tab.1 Statistical results of quantitative precision indexes of two GAM models (1989—2013)

精度指标模型1模型2均值最小值最大值标准差均值最小值最大值标准差平均绝对误差(mm)122.487.6176.928.6117.975.9149.121.6平均相对误差(%)10.06.013.00.018.05.011.00.02均方根误差157.3114.7225.438.7154.098.9213.533.1相关系数0.710.520.880.210.820.620.930.19

根据各年份模型精度指标，选择精度较高的2007和2009年，绘制各站点年降水量GAM拟合值和实测值散点图(如图2)。1989—2013年两个模型的MAIC值和修正后的R2值见图3。

图2 两个模型实测-拟合值对比(2009年)Fig.2 Comparison between measured and fitted annual rainfall in two GAM models in 2009

图3 两个模型MAIC值及修正R2值比较Fig.3 Comparison between MAIC and adjusted R2 values in two GAM models

由图2可知，相比于模型1，加入了高程因子作为预测变量的模型2的点据更为集中在y=x附近，其实测-拟合回归线与y=x基本重合，拟合程度更好。由图3可知：在1989—2013年，除个别年份外，模型2的MAIC值均低于模型1，且修正R2值均高于模型1。这一结果进一步说明，相比于模型1，模型2能更好地拟合出浙西区的降水空间分布，在建模过程中应当考虑高程。

3.2 年降水量影响因子

GAM中各预测因子对降水空间分布的独立影响可以用模型偏效应图[9]来表示。图4～5分别为模型1和模型2的偏效应图(以2013年为例，其他年份类似)。其中，x轴表示各预测因子，y轴代表该预测因子的光滑函数，阴影部分表示95%的贝叶斯置信区间。由图可知：对于模型1，年降水量随纬度的增加明显降低，随经度增加有一定变化，但其规律较为杂乱；对于模型2，年降水量随纬度的增加明显减少，随高程的增加明显增加，但随经度的增加没有明显变化，其趋势线比较平缓。可见，在预测变量中直接加入了高程因子后，经向坐标对年降水量的影响已基本被排除。因此，综合图4和5，浙西区年降水量的主要影响因子是纬向坐标和高程，经向坐标的实际影响较小。

图4 模型1偏效应(2013年)Fig.4 Partial effect graph for GAM based on annual rainfall interpolation model 1 in 2013

图5 模型2偏效应(2013年)Fig.5 Partial effect graph for GAM based on annual rainfall interpolation model 2 in 2013

3.3 年降水空间估算结果

根据率定好的GAM降水估算模型，得到了1989—2013年浙西区1 000 m×1 000 m分辨率的年降水量空间估算结果。根据各年份两个模型的MAIC值和修正R2值，同时考虑到模型的统计精度，选择2007年、2009年展示了浙西区年降水估算结果(见图6)。由图6可知，两个模型的降水空间估算结果均表现为南多北少的总体格局。但模型1只能描述年降水量在空间上的趋势性分布，而模型2加入高程作为预测变量后可以较好地再现年降水量受地形影响的特征，展现年降水量在空间上的局地性、细节性变化规律。两个模型估算结果之差在浙西区的分布主要有以下2个特征：在浙西区南部，模型2的估算结果基本要高于模型1，在浙西区北部则正好相反；两个模型估算结果之差的高值区与浙西区海拔高值区有较好的对应性(见图1(b))。

为验证估算结果的准确性，查阅了《太湖流域及东南诸河水资源公报》[19-20]以及《太湖流域片水情年报》[21-22]，2007年、2009年浙西区年降水量分别为1 319.2和1 459.1 mm，同时，年降水整体表现为由南向北递减，高值区均位于该区西南，高值区的年降水量超过2 000 mm，而2007年模型1、模型2相应值分别为1 370.9 和1 281.6 mm，2009年模型1、模型2 相应值分别为1 560.1和1 498.3 mm；对比两类模型的插值结果可知：模型2的降水空间估算结果与水资源公报中的降水空间分布更为吻合。可以得出：基于模型2的插值结果更能准确反映浙西区降水的空间分布。

图6 基于2个模型的浙西区降水空间估算结果及其计算结果之差Fig.6 Differences between estimated and calculated values based on annual precipitation spatial estimating results in two GAM models

4 结 语

采用太湖流域浙西区33个雨量站1989—2013年的降水观测资料，建立了基于GAM的年降水量空间估算模型，借助多种精度指标，评估了模型的可靠性，同时分析了经度、纬度、高程3种预测因子对降水空间分布估算的影响。结果表明，加入了高程因子作为GAM的预测变量后，浙西区年降水量空间估算精度更高、更稳定，且能够更加准确反映年降水量空间分布受地形影响的细节特征。因此，在浙西区降水空间估算过程中应当对高程的作用合理描述。同时，在空间上浙西区年降水量表现为随纬度增加而减少的规律，但随经度的变化规律不明显。

对GAM在降水空间分布估算中的应用作了初步研究，相比于国内外其他研究成果，重在分析地形因子对GAM降水空间估算模型精度的提升作用，为地形复杂的山丘区降水空间估算及水资源评价提供了一定参考。在后续研究中，将进一步开展基于GAM的月降水量、日降水量空间分布估算研究；此外，可考虑在预测因子中进一步增加坡度和坡向等地形因子，以提高模型的预测合理性。

[1]胡庆芳,杨大文,王银堂,等.赣江流域高分辨率卫星降水数据的精度特征与时空变化规律[J].中国科学(技术科学),2013,43(4)：447- 459.(HU Qing-fang,YANG Da-wen,WANG Yin-tang,et al.Accuracy characteristics and changing regularities of high-resolution satellite precipitation data in Ganjiang Basin[J].Sci China Tech Sci,2013,43(4)：447- 459.(in Chinese))

[2]钱永兰,吕厚荃,张艳红.基于ANUSPLIN软件的逐日气象要素插值方法应用与评估[J].气象与环境学报,2010,26(2)：7- 15.(QIAN Yong-lan,LV Hou-quan,ZHANG Yan-hong.Application and assessment of spatial interpolation method on daily meteorological elements based on ANUSPLIN software[J].Journal of Meteorology and Environment,2010,26(2)：7- 15.(in Chinese))

[3]舒守娟,喻自凤,王元,等.西藏地区复杂地形下的降水空间分布估算模型[J].地球物理学报,2005,48(3)：535- 542.(SHU Shou-juan,YU Zi-feng,WANG Yuan,et al.A statistic model for the spatial distribution of precipitation estimation over the Tibetan complex terrain[J].Chinese Journal of Geophysics,2005,48(3)：535- 542.(in Chinese))

[4]石朋,芮孝芳.降雨空间插值方法的比较与改进[J].河海大学学报(自然科学版),2006,33(4)：361- 365.(SHI Peng,RUI Xiao-fang.Comparison and improvement of spatial rainfall interpolation methods[J].Journal of Hohai University(Natural Sciences),2006,33(4)：361- 365.(in Chinese))

[5]GOOVAERTS P.Geostatistical approaches for incorporating elevation into the spatial interpolation of rainfall[J].Journal of Hydrology,2000,228(1)：113- 129.

[6]徐天献,王玉宽,傅斌.四川省降水空间分布的插值分析[J].人民长江,2010,41(10)：9- 12.(XU Tian-xian,WANG Yu-kuan,FU Bin.Interpolation analysis of precipitation spatial distribution in Sichuan Province[J].Yangtze River,2010,41(10)：9- 12.(in Chinese))

[7]胡庆芳.基于多源信息的降水空间估计及其水文应用研究[D].北京：清华大学,2013.(HU Qing-fang.Rainfall spatial estimation using multi-source information and its hydrological application[D].Beijing:Tsinghua University,2013.(in Chinese))

[8]潘旸,沈艳,宇婧婧,等.基于贝叶斯融合方法的高分辨率地面-卫星-雷达三源降水融合试验[J].气象学报,2015,73(1)：177- 186.(PAN Yang,SHEN Yan,YU Jing-jing,et al.An experiment of high-resolution gauge-radar-satellite combined precipitation retrieval based on the Bayesian merging method[J].Acta Meteorologica Sinica,2015,73(1)：177- 186.(in Chinese))

[9]余黎,雷相东,王雅志,等.基于广义可加模型的气候对单木胸径生长的影响研究[J].北京林业大学学报,2014,36(5)：22- 32.(YU Li,LEI Xiang-dong,WANG Ya-zhi,et al.Impact of climate on individual tree radial growth based on generalized addictive model[J].Journal of Beijing Forestry University,2014,36(5)：22- 32.(in Chinese))

[10]AALTO J,PIRINEN P,HEIKKINEN J,et al.Spatial interpolation of monthly climate data for Finland：comparing the performance of kriging and generalized additive models[J].Theoretical and Applied Climatology,2013,112(1-2)：99- 111.

[11]HUANG C,ZHENG X,TAIT A,et al.On using smoothing spline and residual correction to fuse rain gauge observations and remote sensing data[J].Journal of Hydrology,2014,508：410- 417.

[12]HASTIE T J,TIBSHIRANI R J.Generalized additive models[M].CRC Press,1990.

[13]冯国双,陈景武.广义可加模型及其 SAS 程序实现[J].中国卫生统计,2007,24(1)：82- 84.(FENG Guo-shuang,CHEN Jing-wu.Generalized addictive model and its SAS program implementation[J].Chinese Journal of Health Statistics,2007,24(1)：82- 84.(in Chinese))

[14]WOOD S,WOOD M S.Package ‘mgcv’[Z].R package version,2015：1.7- 29.

[15]WOOD S N.mgcv：GAMs and generalized ridge regression for R[J].R News,2001,1(2)：20- 25.

[16]杨贵军,孟杰,王双喜.基于赤池信息准则的分类回归决策树剪枝算法[J].计算机应用,2014,34(增2)：147- 150.(YANG Gui-jun,MENG Jie,WANG Shuang-xi.Pruning algorithm of classification and regression decision tree based on Akaike information criteria[J].Journal of Computer Applications,2014,34(Suppl2)：147- 150.(in Chinese))

[17]ALLEN D M.The relationship between variable selection and data agumentation and a method for prediction[J].Technometrics,1974,16(1)：125- 127.

[18]李恒鹏,刘晓玫,黄文钰.太湖流域浙西区不同土地类型的面源污染产出[J].地理学报,2004,59(3)：401- 408.(LI Heng-peng,LIU Xiao-mei,HUANG Wen-yu.The non-point output of different landuse types in Zhexi hydraulic region of Taihu basin[J].Acta Geographica Sinica,2004,59(3)：401- 408.(in Chinese))

[19]水利部太湖流域管理局.2007年度太湖流域及东南诸河水资源公报[R].上海：水利部太湖流域管理局,2008.(Taihu Basin Authority of Ministry of Water Resources.Water resources bulletin of Taihu Lake basin and southeast China basins in 2007[R].Shanghai：Taihu Basin Authority of Ministry of Water Resources,2008.(in Chinese))

[20]水利部太湖流域管理局.2009年度太湖流域及东南诸河水资源公报[R].上海：水利部太湖流域管理局,2010.(Taihu Basin Authority of Ministry of Water Resources.Water resources bulletin of Taihu Lake basin and southeast China basins in 2009[R].Shanghai：Taihu Basin Authority of Ministry of Water Resources,2010.(in Chinese))

[21]水利部太湖流域管理局.2007年太湖流域片水情年报[R].上海：水利部太湖流域管理局,2009.(Taihu Basin Authority of Ministry of Water Resources.Hydrological information annual report of Taihu Lake basin in 2007[R].Shanghai：Taihu Basin Authority of Ministry of Water Resources,2009.(in Chinese))

[22]水利部太湖流域管理局.2009年太湖流域片水情年报[R].上海：水利部太湖流域管理局,2010.(Taihu Basin Authority of Ministry of Water Resources.Hydrological information annual report of Taihu Lake basin in 2009[R].Shanghai：Taihu Basin Authority of Ministry of Water Resources,2010.(in Chinese))

Precipitation spatial estimation based on generalized additive model and its preliminary application

WANG Lei-zhi1，HU Qing-fang1，WANG Yin-tang1，HU Yan2，LIU Yong1

(1.StateKeyLaboratoryofHydrology-WaterResourcesandHydraulicEngineering,NanjingHydraulicResearchInstitute,Nanjing210029,China; 2.HydrologyBureau(InformationCenter),TaihuBasinAuthority,Shanghai220000,China)

With the aid of the observation data obtained from 33 precipitation stations during 1989—2013 in the western region of Zhejiang Province within the Taihu Lake watershed,a predictive model for the annual precipitation spatial estimation has been developed based on the generalized additive models (GAM),by use of them an annual precipitation data set with 1 000 m×1 000 m spatial resolution in the study area was gained.By several statistical indexes making analysis of influences of different predictive variables on the precipitation spatial estimation results,the accuracy and reliability of GAM based on the interpolation method were evaluated,especially comparing the differences between two kinds of the models in estimating the annual precipitation space whether or not it selects an elevation as the predictive variable.The analysed and estimated results show that,after taking the elevation as the predicative variable,whether the GAM performance or the precipitation spatial estimation accuracy were improved apparently,which can reasonably reflect the precipitation spatial distribution influenced by the topographic factors in the study area.In addition,it is found from the GAM partial effect graphs that the annual precipitation amount obviously decreased with the increase of the latitude and remarkably increased with the increase of the altitude in the study area,however,there were no obvious changes in the annual precipitation amount with the increase of the longitude.

precipitation spatial estimation; annual precipitation amount; generalized additive model; the Taihu Lake watershed

10.16198/j.cnki.1009-640X.2016.06.007

王磊之，胡庆芳，王银堂，等.基于广义可加模型的降水空间估算模型[J].水利水运工程学报,2016(6):44-51.(WANG Lei-zhi,HU Qing-fang,WANG Yin-tang,et al.Precipitation spatial estimation based on generalized additive model and its preliminary application[J].Hydro-Science and Engineering,2016(6):44-51.)

2015-09-24

国家自然科学基金资助项目(51109136)；水利部科技推广计划项目(TG1528)；水利部公益性行业科研专项项目(201301075，201501014)

王磊之(1991—)，男，江苏句容人，博士研究生，主要从事水文水资源方面的研究。E-mail:wanglz@nhri.cn 通信作者：胡庆芳(E-mail:hqf_work@163.com)

P332.1

A

1009-640X(2016)06-0044-08