空间分析能力是GIS区别于其他信息系统的主要特征^[1]，利用空间分析功能来进行辅助决策是GIS的核心功能。如今，许多商用软件已经拥有相对完善的空间分析功能，但是开源软件对于空间分析中的许多功能并不支持。因此，要基于开源软件进行系统开发，首先需要实现系统所需的空间分析功能。

服务范围分析是网络分析的组成部分之一，主要的功能是计算服务设施对其周围的作用范围和影响程度^[2]。在实际应用中，可以通过缓冲区分析来计算服务范围，但是通过欧氏距离来表达可达性只能作为概略的结果，往往要比实际的服务范围要大^[3]。真正的服务范围计算，其设施点只能够沿着路网拓扑图形中拓扑邻接和拓扑关联的网络节点、网络边进行由近及远的传播^[2]，符合要求的网络节点、网络边的集合就是设施点在特定阻值下的服务范围。本文在最短路径算法的基础上，改进搜索过程中初始点的选取顺序，根据实际需求增加障碍几何图形，对分析结果用最小多边形进行包围。利用开源函数库Geotools进行算法实现，并将其应用于系统开发中。

1 地理网络中心服务范围算法

依据服务范围的概念求算大范围网络密集的路网时存在严重的效率问题^[4]，本文借鉴Dijkstra算法的搜索过程^[5]计算服务范围，按照距离由小到大找出符合要求的各个节点，在寻找节点的同时，筛选出符合要求的网络边。按照上述思想，通过一次遍历就能求出服务范围，保证了不会产生漏点漏边的情况，并且本算法能够根据实际情况添加障碍几何图形，使分析过程更加符合实际情况。

1.1 求算网络节点集合、网络边集合

假设已经根据路网数据构建出网络数据集，如图 1所示，A、B、C、D、E、F为网络数据集中的节点，AB、AC、BC、BD、CD、CE、DE、DF、FE为网络数据集中的网络边。

为方便叙述，本文将某些相关概念用英文加以表示，设定地理网络集合表示为S=(NE，c)，其中N代表地理网络集合中节点的集合，E代表地理网络集合中网络边的集合，c代表服务范围分析中设定的设施点，中心点阻值为dist，e_ij表示节点n_i和节点n_j之间的网络边，其阻值为w_ij，r_ic表示节点n_i和中心设施点c之间的最短路径，则中心点服务范围表示为：

$ \begin{array}{l} S = \{ {n_i}|{r_{ic}} \le dist, {n_i} \in N\} \cup \\ \{ {e_{ij}}|{r_{ic}} + {w_{ij}} \le dist, {e_{ij}} \in E\} \end{array} $

(1)

已经遍历点集合为tempNodes，另外根据实际情况需要，加入障碍几何参数barrers，算法步骤如下：

1) 求出距离中心点最近的节点，作为起始点，并且将起始点的阻值设为0，加入tempNodes。

2) 将起始点加入服务区点集合N，获取在起始点周围的网络边集合edges。

3) 遍历网络边集合edges，每一条边e_ij的阻值为D=r_ic + w_ij，如果D < dist，并且E中没有该网络边，则将其加入E。在遍历网络边的同时，获得该网络边两端除起始点的另一个节点n_j，如果N中不包含该节点，则将n_j及其阻值D加入tempNodes，如果tempNodes中已经存在n_j，则将阻值D改为与之前相比较小的值。

4) 将起始点从tempNodes中删除，求算tempNodes中阻值最小的网络节点作为起始点，如果该点阻值大于dist或者tempNodes为空，则遍历结束，否则执行步骤2)。

由于加入了障碍几何图形barrers，所以在步骤3) 计算阻值时，如果网络边与障碍几何图形相交，则该网络边的权值w_ij为无穷大，就不会走这条路径，实现对障碍几何图形的避让。该算法步骤4) 在寻找新的起始点时不是直接从tempNodes中取，而是取了其中的最小值，由于网络边的阻值为r_ic + w_ij，其中r_ic为当前的最小值，w_ij是网络边固有的权值，所以保证了网络边在第一次被访问时就取得了最小值。比较以前的最短路径算法^[3]，该算法不用记录网络边的前置节点、后置节点，只用一次遍历就能取得网络节点和网络边集合，相比广度优先算法^[4]，需要大量的后续处理来解决“漏边”现象，该算法大大减少了算法的复杂度。

1.2 构造服务区多边形

如图 2所示为某市的交通网络图，中心原点为某快递集散中心，分析该中心5 000 m范围之内的服务范围，按照上述算法得出的分析结果如图 2(a)中加粗黑线所示，对比图 2(b)ArcMap的分析结果显示并不直观。若按照最小凸包多边形来包围网络点集合和网络边^{[6, 7]}，如图 2(c)所示，则分析结果中又会包含不符合条件的网络节点和网络边，所以并不是合适的选择。因此，本文在原有凸包多边形算法的基础上进行改进，得到包围网络节点集合和网络边集合的最小多边形。

利用§1.1计算出的网络节点集合N、网络边集合E分析的主要步骤如下：① 遍历节点集合N，求出Y坐标最小的点fNode；② 遍历网络边集合E，求出以fNode为起始点、极坐标最大的边fEdge，以fEdge作为起始边，并且放入areaEdges集合中；③ 按照逆时针方向，依次求出E中与起始边相交的夹角最大的边作为起始边，并且放入areaEdges集合中，将新找出的边作为起始边，循环直至返回起始节点fNode；④ 将areaEdges集合中的网络边按照存储时的顺序构成多边形。其中关键的步骤就是找到起始网络边，通过步骤①、② 保证了起始网络边的准确性和唯一性。

如图 2(d)为最终的分析结果，和ArcMap的分析结果基本重合，并且与ArcMap效率相同，能够满足现实的分析需求，多边形结果相较图 2(a)更加直观。

2 系统实现

随着开源GIS的发展，一些开源软件在功能和性能上已经和商业GIS软件相差无几，同时开源GIS安全性高，用户定制能力强，自由度大，开发灵活，开源免费，这些特点使其有着广阔的发展前景^[8]。本文利用开源软件设计实现了基于B/S架构的系统，系统架构如图 3所示。

数据存储层采用开源数据库PostGIS对路网和底图等空间数据进行存储和管理；业务逻辑层中采用开源地图服务器geoserver从数据库中取出数据进行渲染发布地图服务；表现层利用开源JS库openlayers3进行地图和分析结果的展示。

服务范围分析位于业务逻辑层中的空间分析服务器中，空间分析服务器采用Java语言在Geotools基础之上进行二次开发，实现数据读取和拓扑构建。Geotools中的几何运算主要依赖于JTS(Java Topology Suit)^[9]，JTS实现了基于二维空间数据模型分析算法，并且利用这些模型分析算法提供了7种空间几何运算。本文主要用到其中的Interaction来判断障碍几何图形与网络边是否相交^[10]。

前端的分析功能界面如图 4(a)所示，分析时，首先选择服务设施点，可以从设施类型中选择在地图上自动定位，也可以在地图上点选中心点的坐标；选择分析方式，既可以是距离成本，也可以是时间成本；然后输入中心点阻值的大小；最后可以根据实际情况添加障碍几何图形(点、线、面)。点击确定进行分析，空间分析服务器通过前端传来的参数和从数据库中读取的数据进行分析，将分析结果返回前端由openlayers3进行解析，并将分析结果展现在地图上。

本文选取某市的路网数据，选定设施类型为医院(图 4(b)中的蓝色点)，阻值的类型为距离，自定义成本为5 000 m，并且按照实际分析需求在部分不能通行的区域添加障碍面和障碍线，分析结果如图 4(b)所示。

3 结束语

本文以计算中心服务范围为前提，在基于Dijkstra思想的算法基础上进行改进，使算法更加简便高效，根据实际情况加入有障碍几何图形的服务范围分析，同时用一种更为准确、直观的方式显示分析结果。在开发方式上，主要采用开源中间件Geotools进行Java开发，依托各种开源软件进行系统实现。相比商用软件开发，采用开源技术能够避免对二次开发组件和运行许可的依赖^[10]。经验证，分析结果在功能和性能上都与商用软件基本无异，能够满足实际的分析需求。