基于IGGⅢ的地理加权回归模型研究

于志英，张福浩，仇阿根，赵阳阳

(中国测绘科学研究院，北京 100830)

地理加权回归(geographically weighted regression,GWR)是地学领域常用的空间分析方法，被广泛用于空间非平稳性探测。其基本思想是将采样点的空间位置嵌入回归参数中，利用局部加权最小二乘方法逐点进行参数估计。通过构建回归模型，探究空间模式背后的影响因素或预测空间现象的发展情况。

诸多学者对地理加权回归模型进行了深入研究。文献[1]提出了混合地理加权回归模型，将回归模型中随空间位置变化的参数称为局域参数，不受空间位置影响的参数称为全域参数，并采用迭代的方法对参数进行近似估计。文献[2]在此基础上，对混合地理加权回归模型进行推导，得到了全域参数和局域参数的精确表达。文献[3]将时间要素融入地理加权回归模型，提出了时空地理加权回归。文献[4]对时空地理加权回归模型估计方法、核函数选择、因子选择、多重共线性检验和参数估计过程进行了详细介绍。文献[5]提出了基于半监督的地理加权回归方法，利用有标记样本训练无标记样本，选择置信度高的结果扩充有标记样本，解决了地理加权回归样本量较少情况下模型精度不高的问题。文献[6]针对现有模型无法充分拟合复杂非线性关系的问题，提出了地理时空神经网络加权回归模型，有效拟合时间邻近和空间邻近的非线性融合作用。另有学者从中间过程出发，对地理加权回归模型进行改进。文献[7]将非线性主成分用于回归模型变量选择，既避免了变量间的共线性，又保留了原始影响因素的主要信息。文献[8]用绝对值交叉验证(absolute value cross-validation，ACV)替代交叉验证(cross-validation,CV)[9]进行了最优带宽选取，避免了CV法中二次误差标准放大离群值的影响。

通过以上研究现状发现，各方法从不同角度对地理加权回归模型进行改进，但缺乏异常值检测与处理过程。因此，本文提出基于IGGⅢ的地理加权回归模型，将IGGⅢ权函数用于地理加权回归参数估计过程，对邻域观测值进行降权或剔除。首先对方法原理进行详细介绍，并介绍方法流程，最后将模拟数据和真实数据与GWR、ACV-GWR进行对比试验，利用均方误差(mean square error，MSE)、平均绝对误差(mean absolute error，MAE)和R2作为指标进行评价。

1 基于IGGⅢ的地理加权回归原理

1.1 方法原理

地理加权回归模型为

(1)

式中，y为因变量；x为自变量；(u,v)为采样点位置；βk(ui,vi)为第i个采样点的第k个参数，β的取值与采样点位置有关；ε为随机误差，符合正态分布。

未知参数β的估计采用最小二乘方法实现。当数据中存在离群值时，离群值参与回归点的参数估计过程，其二次误差会主导残差平方和的值，影响参数估计结果。为削弱邻域点中的离群值对回归点参数估计的影响，本文在地理加权回归模型中引入基于权函数的粗差处理方法。常用的权函数有Huber法[10]、Hample法、Turkey法、Danish法、IGG方案[11]、IGGⅢ方案等。考虑Turkey权函数为有界连续函数；Danish法实质上为淘汰法，没有抗差上的论证[11]；IGG法为有淘汰区的M估计，权因子变化平缓[12]，性能优于Huber法、Hample法[12-13]；IGG为跳跃函数，IGGⅢ权函数为连续函数。因此，本文选择IGGⅢ方案[14]中的权函数用于加权最小二乘参数估计过程。该函数采用三段法进行权重定义，对正常段的观测采用最小二乘估计，对可用观测采用权因子降权，权因子取0～1之间的变值，对达到淘汰界的离群值进行剔除。权函数如下所示

(2)

将其绘制成图直观展示，如图1所示。

将空间距离权重和观测点可靠性权重同时纳入地理加权回归参数估计模型，模型表达为

(3)

式中，wij为空间距离权重，在地理加权回归中常采用Gauss和Bi-square两种核函数进行计算。由于不同核函数对模型参数估计的影响相差不大，最优带宽的选取对参数估计结果影响较大[15]，因此，选择最优带宽值以确定合适的空间权重矩阵对模型解算至关重要。常用的带宽选取方法有Akaike信息量准则(Akaike information criterion，AIC)、贝叶斯信息准则(Bayesian information criterion,BIC)及交叉验证等方法，AIC准则和BIC准则通过极大似然估计计算，CV法采用二次误差标准进行计算。

1.2 方法流程

方法流程如图2所示，流程说明：

(1) 构建地理加权回归模型。利用相关性分析和共线性分析选择变量，并构建地理加权回归模型。

(2) 计算最优带宽。利用AIC准则、BIC准则、CV或ACV进行最优带宽选取。

(3) 计算空间权重矩阵。根据最优带宽，利用Gauss核函数、Bi-square核函数或自适应核函数计算空间权重矩阵。

(4) 估计回归模型参数。构建地理加权回归参数估计模型，根据最小化损失函数原则进行参数求解。

(5) 计算残差。依据步骤(4)计算所得的参数值结果，计算因变量y的估计值。因变量观测值与估计值的差值即为残差。

(6) 计算评价指标。计算标准化残差及MSE、MAE、R2等评价指标。

(7) 判断是否存在离群值。通过标准化残差判断模型中是否存在离群值。若存在，利用IGGⅢ计算权因子，构建基于IGGⅢ的地理加权回归模型，迭代进行模型求解，直到模型中不存在离群值时，结束迭代；若不存在，输出结果。

2 试 验

2.1 模拟数据试验

2.1.1 数据生成

本文根据地理加权回归模型特性设计模拟数据。其中，自变量和因变量满足线性回归关系，系数与采样点空间位置有关，具体公式见表1。

表1 模拟数据生成公式

表1中，x1、x2为自变量，服从(0,1)均匀分布；u、v为位置变量，服从[0,20]均匀分布；ε为随机误差。此外，向模拟数据中添加高斯白噪声。

2.1.2 对比试验设置

本文将GWR、ACV-GWR和IGGⅢ-GWR进行对比试验，采用Gauss核函数计算空间权重矩阵，GWR和IGGⅢ-GWR采用CV确定最优带宽，ACV-GWR采用ACV确定最优带宽。ACV计算方法为

(4)

2.1.3 试验结果分析

利用模拟数据对以上3种方法分别试验40次，表2列举部分试验结果。表3展示IGGⅢ-GWR较ACV-GWR、GWR各指标性能平均提升情况。

表2 部分试验结果

表3 各指标平均提升情况

从MSE、MAE、R2性能提升百分比来看，IGGⅢ-GWR比GWR性能分别提升51.14%、23.77%、28.4%，比ACV-GWR分别提升49.96%、22.57%、27.1%。

2.2 真实数据试验

2.2.1 试验数据

本文选用2016年1月至2018年3月北京地区空气质量及其影响因素作为试验数据进行分析验证。计算CO、NO2、O3、PM10、PM2.5、SO2与空气质量指数的Pearson相关系数(见表4)，CO、NO2、PM2.5、SO2与空气质量指数显著正相关，O3与空气质量指数呈较强的负相关，PM10与空气质量指数呈较强的正相关。

表4 各影响因素与AQI间的Pearson相关系数

对包含CO、NO2、O3、PM10、PM2.5、SO2在内的6种污染物进行多重共线性分析，以方差膨胀因子小于2且条件索引小于10为限定条件进行因子选取，最终选择O3、PM2.5、SO2作为影响指标构建回归模型。北京地区空气质量监测站点分布如图3所示。

2.2.2 试验结果分析

分别采用GWR、ACV-GWR、IGGⅢ-GWR进行试验，以MSE、MAE和R2作为指标对试验结果进行评价,见表5。

表5 真实数据试验结果

从MSE、MAE、R2指标性能上看，IGGⅢ-GWR较GWR分别提升12.65%、7.44%、0.37%，较ACV-GWR分别提升11.85%、6.96%、0.34%。

随机选取任意月份绘制空气质量分布和标准化残差分布，比较各模型结果，本文以2017年1月为例进行结果展示，如图4—图7所示。

图4—图7展示了空气质量观测值和不同模型估计结果。从图4—图7可以看出，2017年1月北京市空气质量指数介于88～202之间，GWR估算结果介于29～150之间，ACV-GWR估算结果介于39～182之间，IGGⅢ-GWR估计结果介于87～208之间。从估算结果上看，IGGⅢ-GWR模型估算得到的空气质量指数更符合真实情况。

从空气状况空间分布来看，2017年1月北京市南部地区空气质量相对较差，GWR估计结果显示中部地区空气质量相对较差，ACV-GWR估计结果规律性不明显，IGGⅢ-GWR估计得到的空气质量空间分布情况与观测值更吻合。

绘制不同回归模型计算所得的标准化残差分布，如图8—图10所示。

从图8—图10可以看出，GWR标准化残差计算结果介于0～2.05之间，ACV-GWR标准化从残差计算结果介于0～2.04之间，IGGⅢ-GWR标准化残差计算结果介于0～1.98之间。从整体上讲，IGGⅢ-GWR估计效果更好一些。从局部来看，IGGⅢ-GWR对北京中南部地区空气质量估计效果优于其他两种方法。

3 结 语

本文提出了基于IGGⅢ的地理加权回归模型，将IGGⅢ方案应用于地理加权回归，降低了离群值对参数估计的影响，提高了地理加权回归模型对离群值的抵抗能力。通过模拟数据和真实数据与GWR、ACV-GWR进行对比试验，以MSE、MAE、R2作为评价指标进行验证。试验结果表明，IGGⅢ-GWR可用于空间非平稳性表达与未知量预测，当数据中存在离群值时，基于IGGⅢ的地理加权回归模型拟合效果更好。