基于多元线性回归模型和Anusplin的秦巴山区≥10℃积温空间模拟比较

张善红, 白红英, 齐贵增, 梁 佳, 赵 婷

(1.西北大学 城市与环境学院, 西安 710127; 2.商洛学院 城乡规划与建筑工程学院, 陕西 商洛 726000;3.西北大学 陕西省地表系统与环境承载力重点实验室, 西安 710127; 4.陕西省气象服务中心, 西安 710014)

政府间气候变化委员会(IPCC)第五次评估报告指出，19世纪 80年代到2012年，地球平均温度升高了0.85℃(IPCC，2014年)。在过去的30 a里，每10 a的地表平均气温都高于人类有记录以来的任何10 a，且2000年以来的十几年气温是有史以来人类记录中最高的(IPCC，2013年)。以全球变暖为特征的气候变化及其有可能产生的负面影响正在成为国际社会普遍关注的问题[1]。特别是近些年来，随着自然因素和人类活动对气候变化的双重影响[2]。气候变化已经严重威胁到农业安全、能源安全、生态安全、水资源安全、公共卫生安全等各方面[3]。气候变化会改变农业气候资源的分布，特别是热量资源的分布。随着热量资源的变化，一个地区的农业种植制度、农业结构和农作物产量都会随着热量资源的变化而产生深刻的影响[4]。积温是衡量一个地区的热量资源极为重要的一个指标，也是衡量作物生长发育过程中所需热量条件的一种指标[5-7]。10℃是喜温植物适宜生长的起始温度，日平均气温稳定≥10℃的积温体现喜温作物所需热量资源的高低[8]，是进行农业区划的重要依据。准确了解地形起伏的山区≥10℃积温的空间分布，对山区农业规划和农业布局具有十分重要的意义。但在现实中各气象站点往往受地理条件、技术条件及经济条件等多种因子的限制，各地区的气象站点的空间分布并不均匀，一般以县为单位作为采集点。而气象站点的实时观测数据只能代表其所在位置的气象状况，无法表征整个区域的气象状况，因此寻找适合山区气象要素的空间插值方法，准确快速地获取山区高精度空间气象数据是近年来地理科学的重要任务之一[9]。

对于积温的研究，国内外各学者相继提出来一些研究模型和方法[10-13]，自20世纪90年代，随着GIS技术的高速发展，气温精细空间化研究在国外得到很大发展，国内的学者开始应用GIS技术研究积温的精细空间分布。刘新安等[14]采用空间插值与三维二次趋势面分析相结合的方法，对中国陆地≥0℃，5℃，10℃，15℃积温进行空间化研究(空间分辨率1 km×1 km)，且空间的精度基本达到实用要求；尹洪涛等[15]利用GIS和数学模拟技术，通过小网格推算模型模拟了辽西地区的热量资源，并通过小地形进行订正，生成了辽西地区100 m×100 m的热量资源空间分布图；戴声佩等[16]采用多元线性回归模型的积温空间模拟方法，模拟并研究了华南地区1960—2011年≥10℃积温的时空变化特征，认为多元线性回归模型适合于较大范围的积温数据空间模拟；高蓓等[17]基于多元线性回归模型的空间模拟认为秦岭南北近50年来≥10℃积温的空间变化整体呈增加趋势。也有一些学者通过多系统气候模式预测较大尺度地区21世纪植被的生长季和活动积温的变化[18]。目前，国内外学者使用Anusplin专业气候空间插值软件等新方法形成了一系列较为权威的气温空间化数据集[19]。此方法逐步成为空间插值的主要手段之一，该方法可考虑多个地理因子作为线性协变量，其插值精度优于传统的插值方法，被广泛应用于气候要素的空间插值研究中。

多元线性回归模型模拟积温方法，较传统的插值方法在复杂的空间插值中有着较高的准确度，因此被广泛应用积温的空间插值[16-22]。国内外学者在积温插值方法方面多考虑多元线性回归模型与Kriging，IDW，Spline这些传统的插值方法对比，证明多元线性回归模型较传统的插值方法精度要高。目前，国内外学者运用Anusplin方法对气温插值的研究较多，但运用该方法对积温插值的研究较少。本文选取地形起伏较大的秦巴山区作为研究区，利用研究区内和周边91个气象站点，采用Anusplin插值和多元线性回归模型两种方法，实现研究区内≥10℃积温的空间模拟，并对预留气象站点实测数据进行检验，分析Anusplin和多元线性回归模型方法效果的优劣，评估更适合我国秦巴山区的气象数据插值方法，为复杂山区积温模拟方法的选择提供科学依据，优化插值效果。

1 研究区概况

秦巴山区位于我国中部(102°54′—112°40′E，30°50′—34°59′N)，总面积约为22.35万km2，地跨陕、甘、蜀、豫、鄂、渝5省1市[23](图1)，地势西高东低，海拔高度悬殊。地貌类型多以山地丘陵为主，间有盆地，形成“三山夹两川”格局。全区受季风气候影响，兼有暖温带和亚热带特征，水资源丰富，是汉江、嘉陵江、丹江等主要河流的发源地，年均降水量450～1 300 mm，年均温度12～16℃。植被类型多样且地带特征明显，秦岭主体为暖温带落叶阔叶林为优势的植被类型，秦岭以南的大巴山区为北亚热带常绿—落叶阔叶混交林，是我国人文、地理、气候、生物等南北过渡区，也是气候变化的敏感区和生态环境的脆弱区[24]。

图1 研究区气象站点分布

2 数据与方法

2.1 数据来源

本文选取研究区及周边91个气象站点1960—2019年逐日平均气温观测数据，来源于陕西省气象局、中国气象局和中国气象科学数据共享服务网。为保证研究数据的完整性和连续性，部分站点缺测数据，通过站点已有数据与邻近站点实测数据的一致性检验进行插补，结果见图2。由图2可以看出，缺测数据站点与临近站点的数据一致性良好，R2均大于0.97以上，说明采用临近站点的插补方法是可行的。数字高程(DEM)数据分辨率为30 m×30 m，来源于国家基础地理信息中心。

注：前边站点为所缺数据站点名称，后面站点为临近站点。

2.2 研究方法

2.2.1 ≥10℃积温的计算 本文的界限温度确定为10℃，采用5日滑动平均法确定日平均气温稳定≥10℃起止日期，其初日定义为5日滑动平均气温≥10℃的日期，终日定义为5日滑动平均气温<10℃的日期[25]。初日和终日之间的日平均气温总和即为日平均气温稳定≥10℃期间积温[16]。

2.2.2 基于多元回归模型的积温模拟 以各站点的经度、纬度、海拔这三要素为自变量，运用多元回归模型模拟秦巴山区积温的小网格空间模型，得到秦巴山区积温的空间分布。针对秦巴山区地形地貌复杂的特点，为提高积温的拟合精度，必须计算积温的残差部分。气候资源要素的残差部分，通常代表着该气候要素受微观地形因素影响所产生微小波动[26]。

本研究方法≥10℃积温的空间分布地表现形式为：

Y=F(λ,φ,h)+σ

(1)

式中：F(λ,φ,h)表示受经度、纬度、海拔影响的≥10℃积温；σ为残差项，即小地形因子对≥10℃积温的影响和随机误差。将建模站点≥10℃的积温及其经度、纬度、海拔代入上式，在SPSS 22.0进行建模，得出回归方程如下：

Y=22309.374-88.891λ-215.718φ-1.584h

(2)

该模型的置信水平α为0.001，复相关系数为0.971，R2为0.943，达到建模的精度要求。

2.2.3 Anusplin插值方法 Anusplin是基于一般薄盘和局部薄盘样条函数插值理论。局部薄盘光滑样条是对薄盘光滑样条原型的扩展[27-28]，它除通常的样条自变量外，允许引入线性协变量子模型，如海岸线距离、海拔等多个地理因子作为线性协变量[29]。局部薄盘光滑样条的理论统计模型为：

Zi=f(xi)+bTyi+ei(i=1,2,3,…，n)

(3)

式中:Zi为位于空间i点的因变量;f(xi)为关于xi的未知光滑函数；xi是独立变量;yi为p维独立协变量;b为独立协变量yi的p维系数;ei为具有期望值为0的自变量随机误差。当公式(3)中缺少第2项,即协变量(p=0)时,模型简化为普通薄盘光滑样条;当缺少第1项独立自变量时,模型变为多元线性回归(Anusplin中不允许这种情况出现[30-31])。Anusplin在独立变量、协变量和样条次数方面有多种组合,共有18种模型。经多次试验,本文在对≥10℃积温插值的最优空间插值模型以经度、纬度为自变量,海拔为协变量的三变量局部薄盘光滑样条函数,样条次数设为2,能够保证精度达到最高。

2.2.4 模型验证方法 为检验两种方法对气象数据插值的优势，通过81个建模站点进行积温的空间模拟，获取预留未参加建模的10个气象站点(蓝田、勉县、石泉、安康、山阳、宕昌、西和、西峡、郧县、开县)的模拟值，对10个预留气象站点的模拟值与实测值在SPSS中进行回归分析，进行误差检验，运用表征估计值可能的误差范围的平均绝对误差(MAE)、反映利用数据估值的灵敏度和极值效应的均方根误差(RMSE)、反映插值相对精确性的平均相对误差(MRE)3个参数[28,32]评估基于多元回归模型的空间积温模拟和Anusplin插值方法模拟秦巴山区≥10℃积温插值的效果。具体公式为(4—6)。

(4)

(5)

(6)

式中：Zi为第i个站点的实际观测值；Z′i为第i个站点的估计值；Z′i为第i个站点的估计值；n为参与检验的站点的个数。参数值越小，说明插值方法的准确度越高。

3 结果与分析

3.1 基于多元回归模型的≥10℃积温空间模拟

(1)参数获取：利用研究区DEM数据在ArcGIS 10.2软件上获取模拟所需的研究区经度、纬度、海拔栅格数据。

(2)获取预测栅格图：利用获取的3个模拟参数和回归模型公式(2)求得研究区1960—2019年日平均气温稳定≥10℃积温期间积温的预测栅格图。

(3)模拟≥10℃积温空间分布图：利用回归模型模拟的建模站点的残差进行反距离权重(IDW)插值，获得模拟残差空间分布图；根据公式(1)模拟计算，获得研究区1960—2019年≥10℃积温空间分布图(图3A)。

(4)插值验证：10个未参与建模站点的实测数据与模拟结果的数据进行回归分析验证，经检验，在p=0.001水平上，R2为0.959 7，实测值与模拟值的线性一致性较强，满足模拟结果的精度要求，说明基于多元回归模型方法模拟秦巴山区≥10℃的空间积温可行。

3.2 基于Anusplin插值≥10℃积温空间模拟过程

3.2.1 数据处理 应用Anusplin软件，首先需要将≥10℃积温数据通过SPSS软件处理成程序要求的标准格式，包括站点代码、经度、纬度、高程信息和积温等数据集，在变量视图中修改数据名称、类型、宽度等，数据格式详见表1，以研究区内5个气象站点的数据为例。Anusplin软件所需要的DEM数据需要转换成 ASCII形式，利用ArcGIS软件对DEM数据投影变换为WGS_1984_UTM_Zone_47N投影，然后利用转换工具将DEM转为ASCII格式文件。

表1 1960-2019年均≥10℃积温数据格式(AAT10℃.dat)

3.2.2 编写运行脚本，反复检验，确定最优模型 Anusplin插值软件调用 splina.exe和lapgrd.exe程序，首先将插值数据和DEM数据转成ASCII格式文件，然后编写运行需要的两个cmd文件脚本(lapgrd和splina)，并运行脚本形成空间插值结果(.grd)，然后对空间插值结果进行投影。

Anusplin在独立变量、协变量和样条次数方面有多种组合，共有18种模型。本文对日均温≥10℃积温插值的最优空间插值模型以经度、纬度为自变量，海拔为协变量的三变量局部薄盘光滑样条函数，样条次数设置为2，保证精度达到最高。

3.2.3 获得研究区≥10℃积温空间分布图 最优模型的空间插值投影后，以研究区为掩膜进行裁剪，获得秦巴山区≥10℃积温空间分布图，见图3。

3.2.4 插值验证 10个未参与建模站点的实测数据与模拟结果的数据进行回归分析验证，经检验，在p=0.001水平上，R2为0.983 2，实测值与模拟值的线性一致性较强，满足模拟结果的精度要求，说明基于Anusplin方法模拟秦巴山区≥10℃的空间积温可行。

3.3 两种插值结果比较

3.3.1 两种插值方法误差对比分析 利用预留的10个站点的实测数据对两种方法的模拟结果进行精度验证，见图3中精度检验。根据验证点两种方法的模拟值和观测值的散点分布可以看出，多元线性回归模型的R2为0.959 7(p=0.000)，AnusplinR2为0.983 2(p=0.000)，表明模拟结果符合要求，说明基于多元线性回归模型的积温空间插值方法和Anusplin插值方法可行，适合于秦巴山区的积温数据空间化模拟。但AnusplinR2略大于多元线性回归模型R2，表明基于Anusplin的插值方法模拟积温与实测积温的线性一致性更强。

图3 秦巴山区≥10℃积温空间模拟结果分布

对两种方法的模拟结果精度进行对比(表2)，结果表明，研究区域Anusplin≥10℃积温插值的RMSE，MRE和MAE均小于多元线性回归模型，表明Anusplin模拟的精度高，误差更小，更适于山区积温空间数据的模拟。

表2 两种模拟方法精度对比

3.3.2 插值结果的空间分析 由图3可以看出，基于Anusplin方法获得的秦巴山区1960—2019年日均≥10℃积温的空间分布可以真实地反映积温的空间分布状况，即≥10℃的积温分布整体呈现出“南多北少”“东多西少”的分布特征。南北方向表现出自北向南随纬度变化的地带性分布变化的阶梯状分布，即随着纬度的减小而增加，但在秦岭北坡地势平坦海拔较低的地方出现了>4 500℃的温度带；大巴山地区由于海拔的缘故，致使>4 500℃的温度带中间出现其他的温度带。随经度的变化幅度较大，≥10℃的积温分布随着海拔的增加而降低。其最高值出现在重庆云阳，其值在6 000℃以上，最小值分布在研究区西北部以及太白山等高海拔处。

为了解Anusplin插值方法模拟的积温与海拔高程的关系，以秦巴山区100 m高程为间隔，统计秦巴山区≥10℃积温均值，得到整个秦巴山区≥10℃积温随海拔变化情况，见图4。1960—2019年≥10℃积温均值与海拔高度的线性方程为y=-1.445x+5462，即海拔每升高100 m，≥10℃积温降低144.5℃，秦巴山区1960—2019年≥10℃积温随高程的增加呈直线下降趋势。

结合图3和图4可以得出，Anusplin插值方法获得的≥10℃积温插值表面细节突出，可以反映局部地形特征，使≥10℃的积温随海拔梯度变化直观，体现了海拔、经纬度对积温的影响，能真实地反映积温的空间分布状况。

图4 秦巴山区ATT10值与高程的关系

4 结 论

(1)基于Anusplin插值方法的秦巴山区≥10℃积温数据空间模拟较多元线性回归模型插值方法的精度更高。其中R2为0.983 2大于多元线性回归模型(0.959 7)插值方法的模拟结果。并且该方法的模拟结果随经纬度、海拔变化趋势明显，精细度也更加明显，体现了经纬度、海拔对积温的空间影响。

有学者[20]认为基于多元线性回归模型的≥10℃积温空间模拟方法较 Kriging，IDW，Spline等插值方法的精度较优。而本文通过多元线性回归模型与Anusplin插值结果对比，表明基于Anusplin的空间模拟方法能更好地获取山区高精度的空间积温模拟，为复杂山区高精度的积温空间数据的模拟提供了一个新的方法。

(2)通过分析1960—2019年秦巴山区≥10℃积温空间分布情况，可以看出研究区≥10℃积温最高值出现在重庆云阳，最低值出现在西北、太白山等海拔比较高的地区；并且≥10℃积温随纬度的增大而减小，随经度的变化规律性较差，但随海拔的升高呈直线下降趋势。

(3)在进行≥10℃积温空间模拟时，因为气象站一般设置在相对平坦开阔地区，插值模型引进坡度、坡向等因子没有通过检验，因此仅考虑了海拔、经度和纬度的影响。该插值模型没有考虑太阳辐射、坡度、坡向、地形以及天气过程等因子对积温的影响，但这些因素的存在对模拟结果存在一些误差。在以后的研究中，应该改进模型，考虑多个影响因子，也要考虑卫星遥感影像在积温研究中的应用，以提高模拟的精度，优化插值效果。