运用协同克里金方法对空气相对湿度进行插值

　　东北地区是中国最大的商品粮基地和农业生产最具发展潜力的地区之一，同时也是中国重要的工业和能源基地。 东北地区位于北半球的中高纬度，是我国纬度最高的地区，是世界著名的温带季风气候区，是典型的气候脆弱区和受气候变暖影响最为敏感的地区。 近年来针对东北地区气候变化已展开不少研究。 已有研究较多地关注气温、降水等方面，李莎等将时空 Kriging 方法应用于东北地区气温空间插值研究中；贺伟等将 Morlet 小波分析应用于东北地区气温和降水变化趋势研究中；孙力等针对东北地区夏季降水变化的时空分布及变化规律进行研究。 这些研究对于正确认识东北地区气候变化规律具有重要意义。 但是以往的研究大多都是基于气温或降水的单一气象要素，而对东北地区空气相对湿度的研究却较少。 在气候资源中，气温、降水和湿度与农作物及生物的生长发育有密切关系。 气温和降水固然是气候资源中最主要的两个因子，但也不能忽视空气湿度这个气候因子。 空气湿度的大小是形成天气演变主要因素之一， 与生产、生活有密切关系，如农作物病虫害的发生直接取决于湿度的大小；空气湿度过大或过小对人体健康也有很大影响等。

　　本研究对东北三省（黑龙江、吉林、辽宁）的空气相对湿度进行空间插值，利用有限的空气湿度站点数据来估算未知点的空气湿度。 目前已发展了较多空间内插方法，如泰森多边形、克里金（Kriging）内插法、反距离加权平均法、趋势面分析法、多项式回归法等。 这些方法只是局限于观测站点的湿度值，没有考虑到降水量、气温、地形地貌、大气环流等因素对空气湿度的影响。 而空气湿度与降水量具有直接的关系，我国一般通过降水量来划分湿润区和干旱区。800 mm 的年降水量为湿润区与半湿润区界线，400 mm 的年降水量为半湿润和半干旱区界线，200 mm 的年降水量为半干旱与干旱区界线。

　　因此，在进行空气湿度插值时，考虑降水量的影响可以提高插值精度。 为此，本研究针对空气相对湿度，选取不同的样本数，结合降水量的观测值，采用协同克里金方法对空气相对湿度进行插值，并对插值结果的精度与普通克里金插值方法得到的结果进行对比分析与探讨，同时对不同样本数的插值结果进行对比。

　　1、研究区域概况与数据处理

　　1．1研究区域

　　研究的是中国东北三省 87 个气象站点， 站点地处 38．90°-52．97°N，119．70°-132．97°E 之间。 地面观测台站所观测的资料来自 1970－2009 年东北三省日空气相对湿度和日降水量。 由于气象台站观测资料不可避免存在缺测，本研究剔除了缺测较多的台站，共选择 77 个资料较好的观测台站。 空气相对湿度选用不同的样本数，选用两种不同数量的观测值分别是 77 个和 46 个， 而降水量共选用 77 个台站的观测值。 具体情况如图 1 所示。

　　1．2数据处理

　　结合降水量对东北三省空气相对湿度进行插值，以得到更加丰富和准确的空气湿度信息。 利用数据相对较全的 1970－2009 年的空气相对湿度观测数据， 也就是对观测站的空气湿度的 40 年平均值分别进行普通克里金插值和协同克里金插值，并分别选取不同的样本数量，对不同样本数量的两种不同空间插值方法进行比较。

　　40 年的平均空气相对湿度分别选取 46 个样本数和 77 个样本数时的统计结果（表 1）。 从表 1 中可以看到，通过 Kolmogorov－Smirnov 法对两种不同样本数的平均空气相对湿度分别进行正态检验的p－value 分别为 0．236 0 和 0．287 9， 均大于 0．05，因此认为符合正态分布。

　　

　　当选择多个变量信息进行插值时，第一信息有限时，合理有效地选择第二信息辅助插值是提高预测精度的关键。 辅助因子的选择需要符合较易获取、与第一信息关联度强等特点。 根据研究目的，将空气湿度作为主变量， 降水量作为辅助变量，将空气湿度和降水量这两种变量进行 Person 相关系数分析。 从表 1 可以看到，当空气湿度选取 77 个观测值时，这两种变量的相关系数为 0．548，同时达到了 0．01 显著水平；当空气湿度选取 46 个观测值时，这两种变量的相关系数为 0．596， 并且达到了 0．01的显著水平， 因而空气湿度与降水量密切相关，可以将降水量作为提高空气湿度插值精度的辅助变量。

　　2、克里金插值及降水量信息利用

　　空间插值是通过已知样本值估计未知点的属性值或由不规则观测点数据生成连续表面或等值线。 它的本质是根据给定的一组或多组离散点的属性值与空间位置拟合一个函数，该函数能充分反映样本点属性值与空间位置间的数学关系，从而推断研究区域范围内未知空间点的属性值。 空间插值结果的精度受到插值模型、模型算法、尺度范围和样本代表性等多种因素的影响。

　　2．1区域化变量

　　当一个变量呈空间分布时，就称之为区域化变量。 这种变量反映了空间某种属性的分布特征。 地质、海洋、土壤、气象、水文、生态等领域都具有某种空间属性。 区域化变量具有双重性，在观测前区域化变量 Z（x）是一个随机场，观测后是一个确定的空间结构函数值。

　　区域化变量具有两个最显著的特征，即随机性和结构性。 一是区域化变量 Z（x）是一个随机函数，多数时候它具有局部的、随机的、异常的特征；其次是区域化变量具有一般的或平均的结构性质，即变量在点 X 与偏离空间距离为 h 的点 X＋h 处的随机量 Z（x）与 Z（x＋h）具有某种程度的自相关性，而且这种自相关性依赖于两点间的距离 h 与变量其他结构特征。 在某种意义上说这就是区域化变量的结构性特征。

　　2．2普通克里金方法

　　克里金法是地统计学的主要内容之一，从统计意义上说，是从变量相关性和变异性出发，在有限区域内对区域化变量的取值进行无偏、最优估计的一种方法；从插值角度讲是对空间分布的数据求线性最优、无偏内插估计的一种方法。 首先假设区域化变量 Z（x）满足二阶平稳假设和本征假设，其数学期望为 m，协方差函数 c（h）及变异函数 γ（h）存在。 即

　　假设在待估计点 x0的临域内共有 n 个实测点，即 x1，x2， … ，xn。 那 么 ， 普通克里 金法 （OrdinaryKriging）的插值公式：

　　

　　式中，Z*（x0）为待估计空气湿度栅格值 ；n 为用于空气湿度插值站点数；Z（xi）为第 i 个站点平均湿度值；λi为赋予第 i 个站点平均湿度的一组权重系数。 引入拉格朗日系数 u 进行推导可得：

　　

　　因此，权重系数（λi，i＝1，…，n）和拉格朗日乘数u 可由以上式（5）求得，从而由式（4）可得研究区域内任意点的插值估计。

　　2．3引入降水量信息的协同克里金方法

　　协同克里金法（CoKriging）是普通克里金法的扩展，把区域化变量的最佳估值方法从单一属性发展到两个或两个以上的协同区域化属性，用一个或多个辅助变量对所感兴趣的变量进行插值估算，这些辅助变量与主要变量都有相关关系，并且假设变量之间的相关关系能用于提高主要预测值的精度。 当研究区内的辅助信息（如降水量）较容易获取且变化平稳时，可将这类信息作为辅助影响因素引入协同克里金方法， 用协同克里金法对东北三省年均湿度进行插值时， 将年均降水量作为一个协变量有利于区域估值的结果， 在计算中要用到两属性各自的变异函数和协变异函数。 在二阶平稳假设下，其期望为：

　　

　　式中，Z2,CK*（x0）为 x0处空气湿度估算值；Z2（x2j）为各点空气湿度；λ2j为赋予各个站点空气湿度的一组权重系数；Z1（x1i）为各点降水量；λ1i为赋予各个站点降水量的一组权重系数；N1、N2分别为降水量和空气湿度用于湿度插值站点数，其中 N1＞N2。 引入两个拉格朗日系数 u1和 u2进行推导可得：

　

　　式中 γ11和 γ22分别是 Z1和 Z2的变异函数模型，γ12和 γ21是这两个变量的变异函数模型， 其中，γ12（h）＝γ21（h）求解线性方程组式（9）即可获得权重系数（λ1i，i＝1，2，…，N1；λ2j，j＝1，2，…，N2）以及两个拉格朗日乘数 u1和 u2可由以上两式求得，从而由（8）式可得研究区域内任意点的插值估计。

　　2．4估值精度评价方法

　　采用交叉验证（Cross－validation）比较不同模型和寻找最佳的插值结果。 在交叉验证过程中， 首先将观测值 Z（xi）暂时去除， 然后通过其他观察值Z（x1），…，Z（xn）和变异函数模型来预测 Z（xi）值， 最后将预测值放回到原始数据中去， 重复以上过程直到对所有观测点进行估值， 所以在所有的观测点处既有实际观测值又有估计值，用统计方法一一比较实际值和预测值。 均方根误差（RMSE）可以用来评价预测值与观测值的接近程度。 利用协同克里金（COK）相对于普通克里金（OK）的均方根误差（RMSE）减少的百分数（RRMSE）表示预测精度的提高程度。

　

　　式中，Z（xi）和 Z*（xi）分别为在 xi处的测量值和预测值；n 为样点数。

　　3、结果与分析

　　从表 2 可以看出，当空气湿度观测值个数为 46时， 协同克里金的 RMSE 值比采用了全部观测值（77 个） 的协同克里金的 RMSE 值降低了 20．58％；比普通克里金的 RMSE 值降低了 7．66％。试验证明，引入降水量作为协变量的协同克里金与普通克里金相比可以提高插值精度。 同时，也可以发现，当空气湿度样本数为 46 个时，比样本数为 77 个的插值结果要好。 因此，在样本采样存在困难时，可以采用协同克里金插值，减少采样点个数，从而节省人力、物力和财力。

　　图 2a 是选取 77 个站点， 经普通克里金内插的东北三省平均湿度分布图， 图 2b 是选取 77 个站点，由协同克里金插值后的东北三省平均湿度分布图。图 3a 是选取 46 个站点，经普通克里金内插的东北三省平均湿度分布图， 图 3b 是选取 46 个站点，由协同克里金插值后的东北三省平均湿度分布图。

　

　　由图中可以看出，本试验区内最大空气湿度值集中分布在东北平原，并沿东、西、北逐渐降低。 最小值则集中在西边的大兴安岭一带。

　　降水量作为协变量来进行 CoKriging 的插值，从表 2 可以看到，与 Ordinary Kriging 相比可以提高插值精度。 当空气湿度的样本数为 77 个观测值时，CoKriging 方法比 Ordinary Kriging 的插值方法的均方根误差 RMSE 值降低了 2．02％； 当减少主变量的样本数， 空气湿度的样本数从中选取较完整的 46个站点时，CoKriging 方法比 Ordinary Kriging 的插值方法的均方根误差 RMSE 值降低了 7．66％。 从中可以发现，引入了高程作为协变量的协同克里金插值方法效果更佳。 另一方面，减少了主变量的（只选取 46 个站点的观测值）的 CoKriging 方法比 77 个站点的 CoKriging 方法均方根误差 RMSE 降低了20．58％，插值结果更好。 因此 ，试验结果表明，CoKriging 不但可以提高插值精度 ， 而且可以在一定程度上减少实地观测样本数，从而可以省时省力节省资金。

　　

　　4、小结与讨论

　　以东北三省东经 119．7°～132．97°，北纬 38．9°~52．97°区域为试验区。 利用该地区 1970－2009 年全年气象站的空气相对湿度观测数据的 40 年平均值， 分别采用普通克里金和协同克里金空间内插方法， 对两种不同方法的实验结果进行对比分析。

　　而且选用不同的样本数量来估计该区的空气湿度，对不同样本数量的空气湿度的协同克里金的插值结果进行对比分析。 试验结果表明，CoKriging 不但可以提高插值精度，而且可以在一定程度上减少实地观测样本数，从而可以节省人力财力和物力。

　　由于空气湿度不仅与降水量有关， 还与气温、日照时间、经度纬度、海拔等都有关系，要想提高空气相对湿度的插值精度，还要进一步考虑这些因素的影响。 同时，由于气象台站分布稀疏，且空间分布不均匀，而且不可避免的存在缺测，导致局部地区误差较大，对空间插值而言，要想提高精度，能够反映数据空间变化趋势和周期的合理的采样设计是必要的前提。 在后续研究中，对缺测数据，可以先进行填补再进行插值，或者采用专门针对这种复杂数据进行处理的统计模型来进行研究。 而且东北地区冬天干燥夏天湿润，在下一步研究中还可以考虑时间的周期性， 在空间插值过程中将时间结合起来，以解决空气湿度的时空分布问题。

　　参考文献：
　　［1］ 姜晓燕，刘树华，马明敏，等．东北地区近百年降水时间序列变化规律的小波分析［J］．地理研究，2009，28．（2）：354－362．
　　［2］ 吴正方，靳英华，刘吉平，等．东北地区植被分布全球气候变化区域响应［J］．地理科学，2003，23（5）：564－570．
　　［3］ 李 莎 ， 舒 红 ， 徐正全 ． 利 用时空 Kriging 进行气温插值研究［J］．武汉大学学报（信息科学版），2012，37（2）：237－241．
　　［4］ 贺 伟，布仁仓，熊在平，等．1961－2005 年东北地区气温和降水变化趋势［J］．生态学报，2013，33（2）：519－531．
　　［5］ 孙 力，安 刚，丁 立，等．中国东北地区夏季降水异常的气候分析［J］．气象学报，2000，58（1）：70－82．
　　［6］ 闫庆武，卞正富，王 红．利用泰森多边形和格网平滑的人口密度空间化研究———以徐州市为例［J］.武汉大学学报（信息科学版），2011，36（8）：987－990．
　　［7］ OLIVER M A，WEBSTER R． Kriging： A method of interpola-tion for geographical information systems ［J］. International Jour-nal of Geographical Information Systems，1990，4（3）：313－332．
　　［8］ LOUIS D E． Pollution models and inverse distance weighting：Some critical remarks ［J］. Computers ＆ Geosciences，2013，52：459－469．
　　［9］ 胡绪福，唐 玮，王 博，等．趋势面分析法应用效果研究———以YL 油田某区块为例 ［J］. 长 江大学学报 （ 自然科学版 ），2010，7（3）：536－538．
　　［10］ 刘小飞，张寄阳，刘祖贵，等．喷灌大田土壤水分不同空间插值方法对比分析［J］. 灌溉排水学报，2008，27（4）：116－118．
　　［11］ 薛 毅，陈立萍．R 统计建模与 R 软件［M］. 北京：清华大学出版社，2007．