全文概述
在空间数据的回归分析中，同一因素在不同地区的影响程度往往会有差异，这种现象称为空间非平稳性。地理加权回归( , GWR)是顾及空间非平稳性的空间统计方法之一[1]，其在普通线性回归( , OLR)模型的基础上加入样本点的地理位置特征，利用局部加权最小二乘方法进行逐点参数估计，通过获取各个回归点随地理位置变化的参数估计值，探测地理空间过程中的空间非平稳性，目前已被应用于社会经济、城市地理、气象、农业、公共卫生等诸多学科领域。文献[2]利用GWR模型检验PM2.5、气溶胶光学深度、气象参数与土地利用之间的关系，得到很好的拟合结果。文献[3]提出地理加权回归克里格法( ，GWRK)，用于检查美国宾夕法尼亚州环境变量与土壤有机碳存量之间的关系，结果表明，与普通回归克里格法相比，GWRK提高了估计精度；文献[4]将GWR用于检验城市化指标与农业景观格局之间的空间变化关系，发现GWR的解释度优于OLR。文献[5]探究了人口、道路、企业等因素对生态用地演变的影响，结果表明，GWR的解释效果优于OLR。文献[6]使用GWR模型可用于解释儿童人口密度、气候因素与手足口病发病率之间的空间异质性关系。文献[7]的研究表明非欧距离度量校准的GWR能够很好地改善房价模型拟合结果。
GWR是一种全样本回归模型，计算时对各个回归点分别求地理权重，过程需实现多层迭代循环计算，整个流程耗费大量的计算时间。为了提高GWR的计算效率，文献[8]研究了大规模GWR在Spark框架上面的实现，发现可利用集群计算机进行并行式的GWR运算，但集群编程开发与测试对于普通开发者而言比较困难。普通商用GPU的多核心高并发的特性促成了GPGPU( GPU)的出现，可将复杂的地理分析计算并行化之后，放在GPU上面进行加速运算，文献[9—11]在此基础上提出了HPC-G和方案，成果显示GPGPU在计算效率提升方面是传统算法的千倍以上。目前，利用GPU实现GWR模型由传统串行模式转变为并行模式的过程尚需深入研究，此外，传统GWR算法用到的权重矩阵的存储方式是全矩阵存储，存储空间是样本数据数量的2次幂，权重矩阵的存储方式也需要进行升级改进。本文提出基于CUDA的快速并行地理加权回归(fast , FPGWR)算法，通过并行计算提高GWR的处理速度和处理能力。
1 经典GWR模型
进行地理现象分析时，由于全局回归未考虑数据的空间异质性，因此结果只能反映出数据在研究区域内的平均水平。学者利用局部多项式光滑思想，提出了地理加权回归模型[1,12-14]。
1.1 GWR模型概念
经典的GWR模型的公式如下
(1)
式中，xim为第i(i=1, 2, …,n)个样本点的第m个自变量；yi为第i个样本点的因变量；(ui,vi)为第i个样本点在研究区域内的坐标(一般使用经纬度)；βm(ui,vi)为第i个样本点的第m个回归系数，是以ui和vi为自变量的函数。
为了简便使用，将式(1)简写为
(2)
将GWR模型相关的变量用矩阵形式表示。其中，样本因变量Y矩阵记为
(3)
样本自变量X矩阵记为
(4)
在研究区域内特定位置i的样本点相对于各样本点j的空间权重wij共有n项，在GWR模型中，第i个点的权重矩阵Wi取对角方阵，记为
(5)
第i个样本点的回归系数估计值
记为
(6)
式中，Xi代表X矩阵的第i行向量。帽子矩阵S记为
(7)
1.2 最优带宽选择准则
由样本值与估计值可计算残差e
(8)
残差平方和记作RSS( sum of )
(9)
根据AIC最小准则寻找最优带宽bw[15-16]，公式为
(10)
利用AIC准则寻找最优带宽bw的原则是，找到AIC最低得分时对应的带宽bw值，搜索方法通常有循环遍历法和黄金分割法。
2 GWR算法原子化
GWR模型的回归过程主要涉及两个方面：第一，依照最小化AIC准则搜寻权函数的最优带宽；第二，利用最优带宽进行结果回归和模型诊断。GWR算法的传统实现方式多基于串行编程计算模式，数据样本量与算法的时间复杂度呈现为幂次或指数关系。本文采用GWR算法原子化的方式将传统串行算法设计模式改为并行计算模式。
2.1 特殊的中间矩阵
普通线性回归模型的矩阵形式如下
(11)
采用普通最小二乘估计可以找到回归系数的估计值
(12)
可得模型的回归结果
(13)
利用普通最小二乘法(OLS)进行无偏估计的回归模型有共同的中间矩阵M，记为
(14)
在逐点回归过程中，计算矩阵M可以作为原子操作。下面详细分解中间矩阵M的计算过程。对于矩阵XT
(15)
矩阵XT与对角方阵Wi相乘具有特殊性，结果矩阵A如下
(16)
令
(17)
式中，矩阵B为p+1阶方阵。在实际应用中，p+1通常在10以内，因此对矩阵B进行求逆计算所花费的时间成本可以忽略。
2.2 算法原子化的实现
区别于传统的全矩阵乘法，算法原子化之后，数据各样本点参与回归计算的过程被固化独立出来，算法计算过程可重复是算法并行化的前提条件。原子过程要解决带宽校准与模型诊断的问题矩阵相乘的前提条件，其中模型诊断需要用到带宽校准的结果。首先将带宽校准过程中产生的AIC得分和估计结果
按照单例模式存储起来；然后，比较AIC得分存储较低AIC得分值及其相对应的
；最后结合Y进行结果的统计分析。在给定测试带宽bw的条件下，算法的原子过程步骤见表 1。
表 1
原子过程算法
原子过程：最小AIC得分策略优选带宽
1. 给定测试带宽(bw)和原子过程索引(z)
2.遍历a=1, 2, …,p+1, 计算
3.遍历b=1, 2, …,p+1, 计算
4.遍历j=1, 2, …,n, 计算
5.Bab+=xja×wij×xjb
6.结束
7.结束
8.结束
9.计算B-1
10.遍历a=1, 2, …,p+1, 计算
11.遍历b=1, 2, …,p+1, 计算
12.+=xzb×Bba-1
13.结束
14.Sz+=×xza×wzz
15.遍历i=1, 2, …,n, 计算
16.+=xia×wzi
17.结束
18.
+=×
19.结束
20.输出Sz,
3 FPGWR
FPGWR算法的实现步骤如下：首先，利用Host与GPU 通信，通知GPU准备开始工作，同时设定公用静态变量并传输至GPU常量内存中；其次，将存储在Host内存中的样本数据输入到GPU全局内存中，本文选择全局内存进行数据存储是因为样本数据量十分庞大，带宽较高的共享内存与本地内存都无法满足空间需求；再次，基于CUDA标准实现原子化核函数的编程；然后，核函数指令以Warp为单位被分配至SM中执行，GPU启动其线程流水线以应对海量线程的执行压力；最后，GPU将计算完毕的模型回归结果与诊断信息反馈至Host，之后将设备资源销毁释放。FPGWR算法的具体工作流程如图 1所示。
图 1基于CUDA的FPGWR算法流程
图选项
4 FPGWR与GWR的对比分析
4.1 内存
FPGWR算法改进了GWR运算过程中的矩阵与数据存储模式，采取了按需存储和矩阵向量化两方面措施。传统GWR方法获得的权重矩阵Wi被存储在n×n的对角方阵之中，空间复杂度为O(n2)。改进算法FPGWR中，原子化过程只会按照需求存储必需的数据，同时权重矩阵与其他矩阵也被分解向量化，空间复杂度被降低至O((p+1)n)(p+1通常小于10)。对比FPGWR与传统GWR的空间复杂度(见表 2)发现：样本量较小时(小于100 000)，两者的内存需求有差异但都处于GPU设备可以承受的范围。在样本量为10 000 000的条件下，传统GWR需要GPU设备分配大约364 TB的显存空间，现在的单GPU设备是不可能符合要求的。相对而言，FPGWR仅仅需要380 MB的显存空间，最大限度降低了存储空间。
表 2
FPGWR与传统GWR内存使用对比
样本点数量
FPGWR
传统GWR
100
3.9 KB
39 KB
1000
39 KB
3.8 MB
10 000
390 KB
380 MB
100 000
3.8 MB
38 GB
1 000 000
38 MB
3.8 TB
10 000 000
380 MB
364 TB
注:所有的数值均采用32位浮点数存储；自变量数量p=9。
表选项
4.2 计算时间
传统GWR的全样本回归过程需要经历两个过程：一是每个样本帽子矩阵Si的独立计算，二是全样本帽子矩阵S的迭代计算。权重矩阵Wi为特殊的对角方阵，在矩阵乘法计算过程中不会增加时间复杂度；矩阵A的时间复杂度为O(1)，矩阵B的时间复杂度为O((p+1)2n)，由于p+1远远小于n，矩阵B-1的计算时间可忽略不计；在固定权重带宽的情况下，帽子矩阵Si的时间复杂度为O((p+1)2n)，全部样本的帽子矩阵S需将Si迭代n次，时间复杂度为O((p+1)2n2)。改进算法FPGWR将各点的回归过程看作一个个独立的线程任务，大幅降了低帽子矩阵S的时间复杂度：矩阵A的时间复杂度为O((p+1)n)，矩阵B的时间复杂度为O((p+1)2n)，由于p+1远远小于n，矩阵B-1的计算时间可忽略不计；帽子矩阵Si以并列相加的方式将矩阵A与矩阵B结合在一起，帽子矩阵Si的时间复杂度为O((p+1)(p+2)n)。
5 试验分析
为了测试FPGWR算法的有效性，采用真实数据对FPGWR和传统GWR的计算时间、内存占用及结果准确性进行对比分析。
5.1 数据来源
夜间灯光遥感数据所反映的夜间灯光亮度可以在一定程度上反映社会经济发展情况。选取贵州省为研究区域，利用FPGWR研究城市兴趣点(point of , POI)及道路对夜间灯光亮度的影响。夜间灯光遥感数据获取自NPP-VIIRS卫星，本文选取了2016年的夜间灯光年合成数据，空间分辨率为500 m，如图 2所示；POI数据获取自百度地图，包括住宅、美食、金融、购物、休闲娱乐等五类；道路数据获取自Open Map。
图 2夜间灯光遥感数据
图选项
5.2 数据预处理
首先，分别计算每类POI点及道路的密度；然后，将研究区域划分为500 m分辨率的网格，获取每个网格中心点上的夜间灯光亮度、POI点密度和道路密度，并将夜间灯光亮度作为因变量，将每类POI点的密度及道路密度作为自变量，得到92万余条数据。由于贵州山区比例高，93%的区域夜间灯光亮度为0，样本的严重不平衡会影响模型的稳定性，因此，本文按照2:1的比例选取夜间灯光亮度≥0的数据，最终得到99 121条数据用于真实数据试验。
5.3 结果与分析
在完全相同的计算环境下，对样本量为99 122条的真实数据分别利用FPGWR算法与传统GWR算法进行计算。结果表明，两种算法的计算精度基本一致，表 3展示了真实数据下两种算法的内存使用、运行时间及FPGWR相对于传统GWR的加速比。
表 3
真实数据下FPGWR与传统GWR算法的内存使用与运行时间对比
FPGWR
传统GWR
加速比(GWR/FPGWR)
内存使用
运算时间/ms
内存使用
运算时间/ms
2.3 MB
22 134.71
36.6 GB
n/a
∞
注:所有的数值均采用32位浮点数存储；自变量数量p=6。
表选项
从表中可以看出，传统GWR需要GPU设备分配大约36.6 GB的显存空间，而FPGWR仅仅需要2.3 MB的显存空间，大大降低了存储空间；FPGWR的运行时间为22 134.71 ms，传统GWR无法完成回归任务。这是因为传统GWR在其运算过程中需要存储完全的权重矩阵，超出了计算机的存储能力。FPGWR优化权重矩阵的存储方式降低了空间复杂度，引入并行计算模式降低了时间复杂度，充分提高了模型运算效率。
6 结语
传统GWR模型采用串行计算模式矩阵相乘的前提条件，在进行大数据量计算时存在内存占用大、计算效率低的问题，难以处理大数据量的地理空间数据。本文提出了FPGWR算法，基于CUDA实现GWR的并行，优化权重矩阵存储方式，将传统的全样本点大循环原子化分解为可并行子过程，显著降低GWR的运行时间和内存占用。为证明FPGWR的实用性，本文使用一套真实数据集进行了试验。结果证明，FPGWR可以提升GWR的运行速度和可处理样本量，并且数据量大时加速效果更为显著。FPGWR的提出使GWR可应用于更大范围研究区域、更高分辨率数据并引入更丰富的数据源。接下来的研究可聚焦于研究FPGWR在不同领域和不同类型数据的应用。
作者简介
作者简介：刘振涛(1977-), 男, 博士生, 高级工程师, 主要研究方向为计算机软件与理论。E-mail:
通信作者：杨毅。E-mail:
权威 | 专业|学术| 前沿
微信、抖音小视频投稿邮箱 |