1 引 言
空间依赖性是地理学第一定律的主要推论，描述了地理现象在空间上不是独立或随机分布的，而是表现出空间自相关性。空间依赖性的普遍存在是进行空间分析与统计的基础，要进行空间依赖性的分析，首先要构建空间权重矩阵，通过距离测度定量表达地理要素的相似程度。不同的距离测度方法能够构建不同的空间权重矩阵，对空间权重矩阵的相关研究，集中于对距离测度方法的分析[1, 2, 3, 4, 5]。二进制邻接矩阵、Queen邻接矩阵、N近邻矩阵、阈值矩阵、距离权重矩阵、Cliff-Ord权重矩阵等是出现较早的空间权重矩阵形式，主要通过邻接与非邻接对空间邻近关系进行描述，使用0-1距离或欧氏距离(Euclidean distance)等定量表达空间邻近程度。高斯距离衰减(Gaussian distance decline)、Tri-cube距离衰减(Tri-cube distance decline)等是较新的空间权重矩阵构建方法，主要通过高斯核函数等方式取代欧氏距离表达空间邻近关系[6, 7, 8, 9, 10, 11]。
空间权重矩阵被广泛应用于空间插值、空间回归、空间点模式识别等空间数据分析领域，空间权重矩阵的定义方式对于空间数据分析结果的准确性有着显著影响，例如在空间插值中利用已知点与未知点的空间权重，进行未知点属性值的预测；在空间回归分析中，构建空间权重矩阵进行回归变量的分析[12, 13, 14]。但现有空间权重矩阵构建方法，大都建立在空间依赖性的基础上，认为在空间上较为邻近的地理要素具有较高的相似度，这种构建思路忽略了地理要素属性值在空间上分布的不均匀性。在空间聚类中，已有研究指出了类似问题，提出了综合考虑地理要素空间与非空间属性的双重距离聚类方法[15]，在遥感影像分类中，也有研究指出应考虑影像单元属性值变化的不均匀性[16, 17, 18]。实际上，地理要素的空间特征除具有空间依赖性外，还具有空间非平稳性。这种局部空间特征的存在，使得地理要素在子空间中具有不尽相同的特征[19, 20, 21]。以欧氏距离或高斯核函数等方法对整个空间地理要素的相似程度进行统一描述并不恰当，应使用局部空间分析方法(local spatial analysis)对空间要素的非平稳性进行分析，在空间权重矩阵的构建中综合考虑空间依赖性与空间非平稳性。本文从这一思路出发，提出了累积相似度表面的概念，以地理要素间的累积相似度作为距离测度，构建空间权重矩阵。 2 方 法 2.1 地理要素相似度
定义地理要素空间T={Di,i=1,2,…,k}，其中地理要素D=(x,y,z)，x、y为地理要素的空间位置，z为属性值，定义p维属性向量z=(z1 z2 … zp)。
定义地理要素的相似度Sij=F(ASij,SSij)，其中ASij为地理要素的属性相似度(attribute similarity)，SSij为地理要素的空间相似度(spatial similarity)。地理要素的相似度由ASij与SSij共同描述，函数F刻画了两者相互作用的方式。
定义地理要素Di与Dj间的属性相似度为ASij=F(zi,zj)，函数F定义了属性相似度的计算方法。函数F可以有多种定义方式，这里F使用高斯核函数(Gaussian kernel)或称径向核函数(radial basis function kernel)
式中，‖zi－zj‖为要素属性值间的欧氏距离，距离值越大，要素间的属性相似度越低，σ是核函数的带宽，用于描述属性相似度减小的速度，当zi=zj时，ASij=1，当‖zi－zj‖趋向于无穷大时，ASij=0。
定义地理要素Di与Dj间的空间相似度为SSij=F(xi,yi,xj,yj)，F使用高斯核函数
如果地理要素的属性在空间上变化均匀，那么属性相似度ASij与空间相似度SSij的变化趋势相似，即满足空间依赖性，空间相似度高的地理要素趋向于具有较高的属性相似度，此时，可通过地理要素的空间相似度SSij度量属性相似度ASij，进而通过已知属性值进行未知属性值的预测，此时相似度函数退化为Sij=SSij。利用欧氏距离、高斯核函数距离或其他现有距离构建空间权重矩阵，实际上就是利用空间相似度SS描述相似度S与属性相似度ASij，距离越大，空间相似度、属性相似度与相似度越小。
由于空间非平稳性的存在，地理要素的属性值在空间上的变化并不均匀，在各局部空间呈现不同的变化特征，利用均匀变化的空间相似度SS对相似度S进行度量已不再准确，在地理要素空间T具有较高空间非平稳性时，这种度量有时甚至会产生较大误差。因此，需要在综合考虑空间相似性与空间非平稳性的基础上，重新定义Sij=F(ASij,SSij)中的函数F，反映空间相似度与属性相似度实际相互作用方式与地理要素空间特征变化趋势。 2.2 累积相似度表面
定义地理要素空间T={Di,i=1,2,…,k}为二维连续表面，存在多种划分连续表面，构建地理要素集的方法，本文采用正四边形进行地理要素空间的划分。同时，为了便于定义地理要素的相邻关系，本文以单个正四边形代表一个地理要素，不考虑多个正四边形组合成地理要素的情况。
在地理要素空间T中，定义地理要素Di的属性相似度表面ASi
式中，ASij为地理要素Di与Dj的属性相似度，对于i=j，ASij=1。
定义地理要素Di的累积相似度表面(accumulated similarity surface，ASS)si={sij,j=1,2,…,K}，其中，sij表示以Di为起点，到达终点Dj的相似度累积最大值。
引入曲线演化理论(curve evolution theory)与快速行进方法(fast marching method)进行累积相似度表面的构建[22]。对于T中的两个地理要素D1、D2，通过AS1构建势能图P：Ω→R+，P在AS1较小的位置具有较大的值，定义P=1/AS1。能量函数表达为Ε:AD1,D2→R+
式中，AD1,D2是D1、D2间路径的集合；γ是两要素间的一条路径；l是长度参数。使能量函数E(γ)具有全局最小值的γ称作D1与D2间的最短路径，记为CD1,D2，D1的最小能量图(minimum energy image)U1:Ω→R+。地理要素Dm在最小能量图上U1的值等于Dm与D1间最短路径CD1,Dm对应的能量函数值
U1满足Eikonal方程
快速行进方法利用逆向差分求解Eikonal方程稳定解
式中，hx、hy为二维连续地理要素空间中x、y方向的要素单元大小。
快速行进方法定义已到达点(alive点)、活动点(active点)、远离点(faraway点)。将起始地理要素D1标记为已到达点，其余点均标记为远离点。开始行进后，计算起始点四邻域点的最小能量并将其均标记为活动点，放入活动点集合。从活动点集合中选择具有最小能量的单元点作为当前行进到的点，改其标记为已到达点，根据该点的最小能量更新其非已到达四邻域点的最小能量(图 1)。循环从活动点集合中选择具有最小能量的单元点，进行单元点最小能量值更新，直至地理要素均被标记为已到达点，最终得到D1的最小能量图U1。更新当前单元点非已到达四邻域点最小能量的方法如式(8)所示
|
| 图 1 快速行进过程示意图 Fig. 1 Fast marching method |
生成D1最小能量图的过程，可以视为水波从D1出发，逐步向外扩散的过程，势能图P可以视为扩散位置的阻力，属性相似度AS1m可视为水波在扩散至Dm位置时的速度1/P。与D1属性相似度AS1较大的Dm，由于扩散速度较快，倾向于具有较小的最小能量U1，但由于扩散是从D1出发由近向远扩散，如果Dm与D1间空间相似度SS1m较小，即Dm与D1空间距离较远，最终得到的最小能量U1仍有可能较大。可见，在扩散进行的过程中，求取的最小能量U1，包含了地理要素间空间依赖性与空间非平稳性等信息。
对于归一化后的最小能量U1，累积相似度表面s1=1－U1，s1m表示从D1出发，沿着属性相似度表面AS1中累积相似度最大(行进速度最大)的方向，行进至D1时所得的相似度累积最大值。 2.3 空间权重矩阵构建
定义地理要素的相似度Sij=F(ASij,SSij)=sij，用累积相似度表面i、j位置的累积相似度sij描述地理要素Di与Dj的相似程度，F为根据ASij与SSij求解sij的函数，即求解最小能量图Um(m=1,2,…,k)的方法，在求解过程中，通过从起始点扩散的先后顺序，加入了空间相似度ASij的影响，通过属性相似度的累积，加入了地理要素空间依赖性与空间非平稳性的影响，体现了地理要素的局部空间特征，可以将累积相似度表面作为一种局部空间分析方法。
地理要素的累积相似度sij反映了地理要素间的相关程度，用其表达空间权重矩阵中的权重值wij，构建基于累积相似度表面sm(m=1,2,…,k)的空间权重矩阵
通过累积相似度表面构建的空间权重矩阵，自身就是对地理要素空间变化趋势较为准确的描述，综合反映了地理要素的空间相似度与属性相似度，与全局空间自相关、局部空间自相关等空间特征描述方法类似，可作为分析地理要素空间特征的重要手段。同时，通过该方法构建的空间权重矩阵，应用于空间插值、空间回归、空间点模式识别等空间数据分析领域时，可以提高空间数据分析结果的准确性。例如，地理加权回归是用于分析空间非平稳性的回归方法，但其使用的空间权重矩阵构建方法仍是基于欧氏距离，利用本文提出的方法构建空间权重矩阵，可以改进地理加权回归方法。再者，通过累积相似度表面构建的空间权重矩阵，可以作为相似性度量，应用于遥感影像分割、分类与特征提取。 3 试 验 3.1 趋势面模拟数据试验
趋势面是用于模拟地理要素在空间上的分布及变化趋势的数学曲面，它抽象并过滤了一些随机因素的影响，使地理要素的空间分布规律明显化。试验建立一次至三次的多项式趋势面，用于分析不同空间矩阵构建方法对插值结果的影响(式(10)~式(12))。
设定趋势面的系数均为1，构建101像素×101像素大小，以(50，50)为中心的1~3次多项式趋势面(图 2(a)~(c))。生成A点(50，50)、C点(25，25)的累积相似度表面(图 2(d)~(f)、(g)~(i))。
|
| 图 2 趋势面模拟数据累积相似度表面 Fig. 2 The accumulated similarity surface of trend surface |
在趋势面模拟数据中，点A的属性值最大，从A点向外扩散，单元点属性值单调递减，与A点距离逐渐增加，属性相似度与空间相似度的变化趋势一致，单元点与A点相似度应逐渐降低，累积相似度表面很好地反映了这一变化趋势：A点与其自身具有最高的累积相似度，越向外扩散，累积相似度越低。从点C出发向外扩散，单元点属性值的变化不再单调，属性相似度与空间相似度的变化趋势不一致，利用点C的累积相似表面同样可以很好地反映这一变化趋势：从点C出发，扩散沿着与点C属性相似度最大的方向快速行进，使得累积相似度较高的区域呈现为环状特征，趋势面中心与四角区域，由于距离C点较远，属性值与C点相差较大，具有较低的累积相似度。
从点A到点B，单元点属性值单调递减(属性相似度降低)，与A点间距增大(空间相似度降低)，与起点A的相似度越低，欧氏距离可以很好地描述AB间单元点的变化趋势(图 3(a))；从点C至点D，单元点属性值变化不再单调，而是先上升再下降，从点C到点E(50，25)，单元点的属性值增大，与起点C的相似度降低，变化趋势反映正确，但从点E到点D，单元点属性值减小(属性相似度增大)，与起点C的欧氏距离增大(空间相似度降低)，空间相似度与属性相似度趋势不一致，而欧氏距离反映从E点到D点，单元点与起点C的相似度单调递减，因此，使用欧氏距离描述CD间单元点的变化趋势并不准确(图 3(b))。
|
| 图 3 相似度-属性值变化曲线 Fig. 3 Attribute value correlated with similarity |
相比于欧氏距离，累积相似度表面可以较好地反映趋势面数据的空间变化趋势：从点A到点B，单元点属性值单调递减，AB间累积相似度降低(图 3(c))；从点C到点E，单元点属性值增大，与起点C的累积相似度降低，到达E点后，与点C欧氏距离仍旧增大(空间相似度降低)，但与点C的属性值差异减小(属性相似度增大)，初始阶段，空间相似度的降低占主导地位，单元点的累积相似度仍旧下降，但随着单元点与点C属性值越来越相似，属性相似度开始占主导地位，对于一次与二次趋势面，单元点与点C的累积相似度开始上升，对于三次趋势面，单元点与点C的累积相似度降低速率变慢(图 3(d))。考虑到常用于地理要素空间拟合的多项式趋势面不会高于三次，累积相似度表面对趋势面数据空间变化趋势的描述比利用欧氏距离描述更为准确。 3.2 数字高程模型数据试验
试验选择ASTER GDEM数据，随机采样获取32个已知点数据(图 4)。
|
| 图 4 ASTER GDEM试验数据与随机采样点 Fig. 4 ASTER GDEM experimental data and random sample points |
利用欧氏距离与基于累积相似度表面的方法构建空间权重矩阵，进行反距离权重插值如下
式中，zi为待插值点的属性值；wij表示zi与zj间的权重，若zj同样为待插值点，则wij=0，d(zi,zj)表示zi与已知点zj间的距离；p表示距离权重，试验取d(zi,zj)=‖zi－zj‖，距离权重p=2。对基于累积相似度表面的空间权重矩阵，定义wij=sij，在计算属性相似度AS时，核函数带宽取10。反距离权重搜索策略为选择与单元点空间权重最大的10个已知点。通过欧氏距离与累积相似度表面构建空间权重矩阵，插值结果如图 5所示。
|
| 图 5 反距离权重空间插值结果 Fig. 5 IDW interpolation results |
试验数据变化趋势复杂，包含3个方向特征显著的高值区域(图 4中区域A、B、C)，1个呈沟壑状的低值区域(图 4中区域D)。利用欧氏距离与基于累积相似度表面方法均可以插值得到3个高值区域，但由于区域B高程为710 m和790 m的采样点相距较远，利用欧氏距离得到的区域B插值结果表现为以两个采样点为中心的圆环状区域，基于累积相似度表面的插值结果可以较为准确地反映区域B高程的空间分布特征。对于低值区域D，利用欧氏距离得到的插值结果同样表现为以插值点为中心的圆环状区域，无法准确反映该区域高程的空间分布特征。基于累积相似度表面的插值结果则反映出了区域D沟壑状的高程值分布特征(图 5)。由此可以看出，相比于基于欧氏距离构建的空间权重矩阵，基于累积相似度表面构建权重矩阵可以很好地反映数字高程模型数据的空间变化趋势，尤其是能较为准确地反映数据空间分布的细节特征。
空间插值的目的在于对待插值点属性值进行预测，试验将图 7中的空间插值结果与试验数据相比较，分析两种方法对单元点属性值预测的准确性。构建残差平方和评价指标
式中，zi为单元点i的插值结果；z0i为单元点i的实际值。计算得到利用欧氏距离插值的指标J1=3.18×1010，基于累积相似度表面插值的指标J2=4.24×109，J1约为J2的7.5倍。基于累积相似度表面的插值相比于利用欧氏距离的插值，具有较小的残差平方和，插值结果更为准确。 4 结论与讨论
本文将地理要素的相似度定义为属性相似度与空间相似度，提出累积相似度表面的概念，引入曲线演化理论和快速行进方法生成累积相似度表面，构建空间权重矩阵。基于累积相似度表面构建的空间权重矩阵，以属性相似度体现了地理要素属性值的相似程度，从起点出发，沿着相似度最高的方向快速行进，行进过程中得到的累积相似度包含了地理要素的局部空间特征。通过趋势面模拟数据与数字高程模型数据的试验分析可知，相比于利用欧氏距离等方法，通过累积相似度表面构建的空间权重矩阵综合考虑空间依赖性与空间非平稳性，体现了地理要素的局部空间特征，能够更为准确地描述地理要素空间特征变化趋势。
目前，尚未出现定量描述空间非平稳性的指标或方法，更多的是通过回归、聚类等方法证明这种性质的存在。本文指出累积相似度表面中包含了空间非平稳性等空间变化趋势信息，能更好地描述地理要素的相似程度，但尚未提出如何进一步分析统计累积相似度表面，获取对空间特征变化趋势的定量描述，如何由累积相似度表面及其建立的空间权重矩阵，定量描述空间非平稳性，需要进一步的分析探讨。
通过累积相似度表面生成空间权重矩阵时，需要计算每个地理要素的累积相似度表面，快速行进方法的计算复杂度为O(NlogN)，则遍历所有地理要素生成累积相似度表面的计算复杂度为MO(NlogN)，M为地理要素的数量，而利用欧氏距离计算空间权重矩阵，只需计算每个地理要素与起始要素的欧氏距离，计算复杂度为O(N)。因此，基于累积相似度表面的方法，在算法运行效率上，不如基于欧氏距离的方法，不利于处理包含大量地理要素的空间数据。后续研究可以从以下两方面进行算法效率的提升：① 各个地理要素的累积相似度表面计算过程独立，可以通过算法并行化实现累积相似度表面的快速计算；② 空间权重矩阵应用于空间插值或空间回归分析时，对于一个地理要素，往往只需要获取与其空间权重最大的m个地理要素，或空间权重大于某一阈值的地理要素，因此可采用相似度搜索的相关算法，生成部分地理要素的累积相似度表面，即可完成空间权重矩阵的构建。
| [1] | ALDSTADT J, GETIS A. Using AMOEBA to Create a Spatial Weights Matrix and Identify Spatial Clusters[J]. Geographical Analysis, 2006, 38(4):327-343. |
| [2] | GETIS A. Spatial Weights Matrices[J]. Geographical Analysis, 2009, 41(4):404-410. |
| [3] | GETIS A, ORD J K. The Analysis of Spatial Association by Use of Distance Statistics[J]. Geographical Analysis, 1992, 24(3):189-206. |
| [4] | CLIFF A D, ORD J K. What Were We Thinking?[J]. Geographical Analysis, 2009, 41(4):351-363. |
| [5] | WANG Hongliang, HU Weiping,WU Chi. Analyzing the Effect of Spatial Weighted Martix on Spatial Autocorrelation—Taking Hunan’s Incoming GAP between Urban and Rural Areas as a Case[J]. Journal of South China Normal University:Natural Science Edition, 2010(1):110-115. (王红亮,胡伟平,吴驰.空间权重矩阵对空间自相关的影响分析——以湖南省城乡收入差距为例[J].华南师范大学学报:自然科学版, 2010 (1):110-115.) |
| [6] | FOTHERINGHAM A S, BRUNSDON C. Some Thoughts on Inference in the Analysis of Spatial Data[J]. International Journal of Geographical Information Science, 2004, 18(5):447-457. |
| [7] | PELLEGRINI P A, FOTHERINGHAM A S. Modelling Spatial Choice: a Review and Synthesis in a Migration Context[J]. Progress in Human Geography, 2002, 26(4):487-510. |
| [8] | FOTHERINGHAM A S, CHARLTON M, BRUNSDON C. The Geography of Parameter Space: an Investigation of Spatial Non-stationarity[J]. International Journal of Geographical Information Systems, 1996, 10(5):605-627. |
| [9] | DEMSAR U, FOTHERINGHAM A S, CHARLTON M. Combining Geovisual Analytics with Spatial Statistics: the Example of Geographically Weighted Regression[J]. Cartographic Journal, 2008, 45(3):182-192. |
| [10] | DEMSAR U, FOTHERINGHAM A S, CHARLTON M. Exploring the Spatio-temporal Dynamics of Geographical Processes with Geographically Weighted Regression and Geovisual Analytics[J]. Information Visualization, 2008, 7(3-4):181-197. |
| [11] | GETIS A, ALDSTADT J. Constructing the Spatial Weights Matrix Using a Local Statistic[J]. Geographical Analysis, 2004, 36(2):90-104. |
| [12] | STEIN A, BASTIAANSSEN W G M, DE BRUIN S, et al. Integrating Spatial Statistics and Remote Sensing[J]. International Journal of Remote Sensing, 1998, 19(9):1793-1814. |
| [13] | HARRIS P, FOTHERINGHAM A S, JUGGINS S. Robust Geographically Weighted Regression: a Technique for Quantifying Spatial Relationships between Freshwater Acidification Critical Loads and Catchment Attributes[J]. Annals of the Association of American Geographers, 2010, 100(2):286-306. |
| [14] | GENG Xiepeng, DU Xiaochu, HU Peng. Spatial Clustering Method Based on Raster Distance Transform for Extended Objects[J]. Acta Geodaetica et Cartographica Sinica, 2009, 38(2):162-167. (耿协鹏,杜晓初,胡鹏.基于栅格距离变换的扩展对象空间聚类方法[J].测绘学报, 2009, 38(2):162-167.) |
| [15] | LI Guangqiang, DENG Min, CHENG Tao, et al. A Dual Distance Based Spatial Clustering Method[J]. Acta Geodaetica et Cartographica Sinica, 2008, 37(4):482-488. (李光强,邓敏,程涛,等.一种基于双重距离的空间聚类方法[J].测绘学报, 2008, 37(4):482-488.) |
| [16] | ATKINSON P M. Spatially Weighted Supervised Classification for Remote Sensing[J]. International Journal of Applied Earth Observation and Geoinformation, 2004, 5(4):277-291. |
| [17] | FOODY G M. Geographical Weighting as a Further Refinement to Regression Modelling: an Example Focused on the NDVI-rainfall Relationship[J]. Remote Sensing of Environment, 2003, 88(3):283-293. |
| [18] | ATKINSON P M, LEWIS P. Geostatistical Classification for Remote Sensing: an Introduction[J]. Computers and Geosciences, 2000, 26(4):361-371. |
| [19] | FOTHERINGHAM A S, BRUNSDON C. Local Forms of Spatial Analysis[J]. Geographical Analysis, 1999, 31(4):340-358. |
| [20] | HARRIS P, FOTHERINGHAM A S, CRESPO R, et al. The Use of Geographically Weighted Regression for Spatial Prediction: an Evaluation of Models Using Simulated Data Sets[J]. Mathematical Geosciences, 2010, 42(6):657-680. |
| [21] | NAKAYA T, FOTHERINGHAM A S, BRUNSDON C, et al. Geographically Weighted Poisson Regression for Disease Association Mapping[J]. Statistics in Medicine, 2005, 24(17):2695-2717. |
| [22] | DESCHAMPS T, COHEN L D. Fast Extraction of Minimal Paths in 3D Images and Applications to Virtual Endoscopy[J]. Medical Image Analysis, 2001, 5(4):281-299. |