空间权重矩阵的那些事(八)-球面距离权重矩阵

前段时间有人向我咨询了根据经纬度计算球面距离的方法，希望我出一篇文章，所以就有了这篇文章。必要文件可通过后台回复「地理经纬度」获取。

我首先想到的是matlab的「jplv7」工具箱里的「distance」函数，下面是具体的程序内容：

function D = distance(xc,yc)
% PURPOSE: Computes the list of pairwise distances for a given set of locations (loc).
% ----------------------------------------------------------
% Usage: D = distance(xc,yc)
% where: xc,yc are vectors of latt-long coordinates for each location
% ----------------------------------------------------------
% Returns: D = (n x n)-matrix of pairwise distances

% Written by: TONY E. SMITH, 2/10/98

n = length(xc) ;  %number of locations
% Start procedure.
X = xc ; %column vector
Y = yc ; %column vector
U = ones(n,1) ; %column vector
XX = X*U‘ - U*X‘ ;
YY = Y*U‘ - U*Y‘ ;
D = (XX.^2 + YY.^2) ;

根据程序内容，距离计算公式是：

$$
D=\,\,\left( x1-x2 \right) ^2+\left( y1-y2 \right) ^2
$$

其中，(x1, y1) 是点 A 的坐标，(x2, y2) 是点 B 的坐标。这明显不是球面距离。实际上，如果经纬度数据如果是未投影的数据（即空间坐标，单位为 °），经过上面的函数计算得到的距离仍然不是地理矩阵，因为其单位不是 mi 或者 km 等距离单位。使用该函数时需要特别注意！总之，只有平面坐标通过上面的函数计算的结果还是真实的地理距离。

通过查找资料，发现matlab本身就存在可以计算球面距离的函数，而且函数名也是「distance」。这样，如果在matlab中已经导入了「jplv7」包的情况下再想使用matlab自带的distance函数，就可能存在冲突。解决方法是，调正函数的优先级。具体操作是：主页 > 设置路径，选中jplv7的所有路径，点击「移至底端」、「保存」。

空间权重矩阵的那些事(八)-球面距离权重矩阵

这样再使用「distance」函数时，系统就会自动调用matlab本身的函数，而不是「jplv7」下的「distance」函数。

具体的计算原理：

空间权重矩阵的那些事(八)-球面距离权重矩阵

图中，空间坐标系 A(x1, y1), B(x2, y2)，我们需要计算的球面AOB距离为弧长。可以通过公式：半径 × ∠AOB来完成。地球赤道半径 6378.140 km，地球极地半径 6356.755 km 。那么，现在面临的问题就是∠AOB的计算了。

matlab本身自带的「distance」函数就可以解决这个问题。这里，取赤道半径6378.1km，下面展示计算过程：

>> D = distance(23.0, 101.1, 23.06, 113.34)

D =

   11.2613

>> 11.2613/180*pi*6378.1

ans =

  1.2536e+003

或者

>> D = distance(23.0, 101.1, 23.06, 113.34, 6378.1)

D =

   1.2536e+03

可以看到。两种方式计算的结果是一致的。

下面，我们将演示怎样用matlab自带的这个函数进行空间权重矩阵的构建。首先，使用ArcGis导出地理坐标系下的经纬度数据（单位为 °），排好序后以列向量的格式导入matlab。具体过程不再演示，相关经纬度数据可通过后台回复「地理经纬度」获取，下面直接贴出代码内容：

% 球面距离空间权重矩阵
E = zeros(33, 33);
for i = 1:32;
    for j = 1 + i:33;
        E(i, j) = E(i, j) + distance(X(i), Y(i), X(j), Y(j), 6378.1);
    end;
end;
wd= E + E‘;

下面给出了，基于赤道半径的球面距离和欧几里得距离法下的空间权重矩阵（这两份数据也已上传到百度云，回复可获得）的部分结果：

空间权重矩阵的那些事(八)-球面距离权重矩阵

球面距离空间权重矩阵

空间权重矩阵的那些事(八)-球面距离权重矩阵

欧几里得距离空间权重矩阵

可见，两者存在一定差距。毕竟，地球并不是标准的球体，以赤道半径作为地球半径本身存在一定误差。所以，在没有特殊要求的情况下，还是建议使用投影后的平面坐标来计算地理距离。

空间权重矩阵的那些事(八)-球面距离权重矩阵

更多精彩内容，扫码关注微信公众号！

空间权重矩阵的那些事(八)-球面距离权重矩阵

Canethui