利用GAS软件和Matlab编程实现空间统计分析.doc

13. 利用GAS软件和Matlab编程实现空间统计分析空间统计分析，目前比较流行的计算软件主要有Asellin等人编写的Space Stat、Geodata等，随着ARCGIS的发展，一些空间统计模块中也嵌套了类似的计算功能，如在ARC/INFO的工作空间WORKSPACE环境，其中的GRID模块中，可以调用命令MORAN来实现全局空间自相关指数，当然这主要是针对栅格数据的计算，对于一些其他格式的数据，可以通过数据的交互格式来转换。另一方面，弄清楚空间统计分析的基本原理之后，也可以自己编程来实现其计算过程。下面主要将如何在软件Geodata Analysis Software直接实现全局统计、局部统计以及在Matlab中编程实现空间统计的过程。例：以书本上的应用实例来说明，图1给出了上海市部分街道（50个样区）的邻居关系，各街道编号、名称、所属区，以及2000年的人口密度数据如表1所示（原始数据为shape格式的矢量数据）。那么，人口分布是否存在空间上的相关关系呢？为了回答这一问题，就需要进行空间统计分析。图1上海市部分街道（样区）的邻居关系表1 上海市部分街道（样区）的人口密度及其ESDA计算结果街道编号街道名称所属区人口密度(单位：人/km)GiZ(Gi)1大桥街道杨浦区26 089.380.072 9-0.606 42四川北路街道虹口区43 155.370.144 50.885 13江浦路街道杨浦区37 526.670.104 10.043 64新港路街道虹口区57 142.640.096 40.345金杨新村街道浦东新区18 434.640.036 6-1.036 76平凉路街道杨浦区33 715.290.092 5-0.575 67宝山路街道闸北区60 764.30.119 40.869 78宜川路街道普陀区45 052.610.053 6-0.201 69嘉兴路街道虹口区45 666.770.127 60.532 110芷江西路街道闸北区51 118.210.110 90.674 111提篮桥街道虹口区35 986.240.1530.181 812天目西路街道闸北区27 606.910.1710.507 713洋泾街道浦东新区17 597.010.083-1.077 314乍浦路街道虹口区51 422.870.157 30.260 815北站街道闸北区63 370.290.209 10.781 216长寿路街道普陀区25 291.350.100 4-0.034 517梅园新村街道浦东新区18 287.020.143-0.345 118外滩街道黄浦区29 876.990.232 80.800 319江宁路街道静安区43 080.80.122 2-0.004 420南京东路街道黄浦区55 084.40.053 3-0.651 221钦洋镇浦东新区4 813.090.020 7-0.642 422石门二路街道静安区57 311.380.118 0.333 23人民广场街道黄浦区32 799.650.237 1.256 624曹家渡街道静安区42 568.480.095 2-0.142 425金陵东路街道黄浦区66 937.740.109 1-0.257 426潍坊新村街道浦东新区32 958.370.020 7-0.642 427小东门街道黄浦区41 752.640.163 50.372 128南京西路街道静安区35 618.770.144 50.030 429豫园街道黄浦区63 508.910.152 1.039 530静安寺街道静安区27 676.190.117 1-0.10231老西门街道黄浦区75 864.040.241 81.339 432淮海中路街道卢湾区71 677.720.180 31.112 833花木镇浦东新区6 453.420.034 6-1.403 534江苏路街道长宁区35 509.790.074 1-0.173 135瑞金二路街道卢湾区38 087.970.154 30.205 236湖南路街道徐汇区29 796.740.107 4-0.289 837塘桥街道浦东新区18 085.240.063 6-0.800 938董家渡街道黄浦区36 815.940.146 9 0.073 139打浦桥街道卢湾区44 442.230.195 3 0.546 440天平路街道徐汇区29 699.730.107-0.297 941半淞园路街道黄浦区41 259.690.149 20.114 242五里桥街道卢湾区27 068.640.100 6-0.421 143南码头路街道浦东新区21 973.130.088 6-0.977 144斜土路街道徐汇区29 112.080.092 7-0.571 645徐家汇街道徐汇区24 479.340.058 3-0.537 846枫林路街道徐汇区41 794.510.074 2-0.5847周家渡街道浦东新区25 452.990.078 9-0.836 948六里镇浦东新区12 026.160.036-1.049 749上钢新村街道浦东新区16 569.860.059 5-0.886 50龙华镇徐汇区10 603.980.067-0.729 41利用Geodata Analysis Software实现将GIS数据（譬如*.shp格式），导入到Geodata 中，该软件系统界面如下图（图2）所示：图2 软件分析界面菜单命令及主要工具按钮如下： close all windows，即关闭所有已经打开的窗口；add centroids，即添加中心点；如果打开数据表，其形式如下：在做空间统计分析之前，首先需要计算空间邻居关系的权重矩阵。为此，可以在上述图形界面引导下，在“Tools”工具栏中，利用“Weights” 工具栏中的“Create”产生一个权重矩阵文件，譬如，对于以上人口分布的各个街道数据，可以产生一个权重矩阵文件density_weights.GWT。然后，利用权重矩阵，就可以做一系列的空间统计分析，譬如：（1）计算全局Morans I并做出Moran散点图，如下：（2）做出LISA显著水平图，如下：（3）还可以进一步通过回归分析，研究变量之间的相互关系。，在图形界面引导下，按步骤操作，就可以实现该软件提供的一些空间统计分析功能。 2.局部G统计计算的Matlab程序除了Geodata Analysis Software软件外，还可以利用Matlab编程进行空间统计分析。譬如，局部G统计的Matlab计算程序如下：function GI,ZG=Getis(Wij,Xij)% 计算G统计值：G值是不包含自己的G统计值% 输入：Wij空间权值矩阵% 输入：Xij研究区域的空间属性数据% 输出：GI 空间局部自相关指标-G统计值% 输出：ZG 对于GIJ的检验Z值的计算结果% GI=Getis(Wij,Xij,N)M,N=size(Wij);GI=zeros(N,1);GIn=0;Xn=0;for m=1:N for n=1:N if n=m GIn=GIn; Xn=Xn; else GIn=GIn+Wij(m,n)*Xij(n); Xn=Xn+Xij(n); end end GI(m)=GIn/Xn; GIn=0; Xn=0;endsave GI GI -ascii;%计算Z值EG=zeros(N,1);VarG=zeros(N,1);ZG=zeros(N,1);WI=0;for m=1:N %计算E（Gi(d)） for n=1:N if n=m WI=WI; else WI=WI+Wij(m,n); end end EG(m)=WI/(N-1); WI=0;endWI=0;Xj=0;Xj2=0;%Yi1=zero(N,1); %计算Yi1%Yi2=zero(N,1); %计算Yi2Wi=zeros(N,1); %计算Wifor m=1:N for n=1:N if n=m WI=WI; Xj=Xj; Xj2=Xj2; else WI=WI+Wij(m,n); Xj=Xj+(Xij(n)*Xij(n); Xj2=Xj2+(Xij(n)*Xij(n)/(N-1); end end Wi(m)=WI;