为加强工程建设的成本效益管理，提高山区农村居民生活和生产水平，《农村饮水安全工程建设管理办法》(发改农经〔2013〕2673号)指出^[1]：农村饮水安全工程建设应当合理布局，优先采取跨村、跨乡镇联片集中供水等方式，实现供水到户，确保工程投资效益.然而，多数山区农村根据行政辖区来划分供水范围，即使在同一行政区内供水范围也只是根据各村户生产生活用水量进行划分，然后采用网络优化的方法实现管网的优化布置，没有协调好供水保障率和供水范围重叠之间的关系，供水管网的空间分布不能满足成本效益管理要求.因此，利用山头地势布设水池的既定供水点的服务范围划分问题，是山区供水工程设计必须解决的基础问题.

区划问题一般出现在水能资源规划^[2-4]、生态河岸带功能区划^[5-8]、滑坡地质灾害危险性区划及洪灾风险区划^[9-11]中，且大多采用GIS技术^{[2, 9-11]}构建区划模型.加权Voronoi图也常被用于城市空间影响范围研究^[12-13]、地市级旅游区划研究^[14]以及城市湿地功能分区模型研究^[15].可见，Voronoi图在区划模型中得到广泛应用.针对山区复杂地形供水问题，文献[16]对地表进行Delaunay三角网格划分，并将三角网格转换成带边权的网络模型，运用改进的Dijkstra算法计算两点间最短管线线程，为优化供水管网提供依据.然而，供水管网不仅有一级管网的优化目标，而且更需要对末端水池的二级管网优化进行研究.本文采用Voronoi图构造末端水池供水区划模型，在Matlab编程语言环境下，依托贵州省纳雍县山区农村饮水安全工程对既定水池供水区划问题进行案例研究，为科学解决末端水池二级供水管网的精细优化问题提供新的思路.

1 基于Voronoi图的既定水池区划模型 1.1 增量法原理

增量法是构造Voronoi图最普遍采用的算法.采用增量法实现既定水池区划优化的基本思想是^[17]：逐步加入水池点，在已构造的Voronoi图的基础上利用局部特性，采取局部修改方法，生成新的Voronoi图，即对于平面离散水池点集P={p₁, p₂, …, p_n}，在生成水池点集P_k={p₁, p₂, …, p_k}(k＜n)的Voronoi图的基础上，加入p_k+1，通过局部修改来构造P_k+1={p₁, p₂, …, p_k, p_k+1}的Voronoi图.如此不断加入新的水池点，最后可得到VD(P).

实现增量法思路的主要步骤是定位和局部修改，以下详细分析.

1.2 新加水池点定位

查询新加入水池点p_k+1在VD(P_k)中的位置，即其所在的Voronoi区域VR(P_i)(1≤i≤k).可通过将p_k+1与p_k中每个水池点的距离进行比较，找到离其最近的水池点p_i，即可确定p_k+1.

1.3 对Voronoi图局部修改

水池点p_k+1加入后，在其周围形成一个由离p_k+1最近的点组成的Voronoi区域VR(P_k+1), 见图 1.p_k+1落在VR(P_i)中，用p_k+1和p_i的中垂线把VR(P_i)分成两部分，每部分的点分别离p_k+1和p_i最近.它和VR(P_i)交于t和q两点.用p_k+1和p_j的中垂线把VR(P_j)分成两部分，同理，每部分的点分别离p_k+1和p_j最近.该线交VR(P_j)于q和r两点.用r代替q循环上述工作，直到t点被垂直平分形成一个多边形为止(图 1虚线多边形).最后删除该多边形中间的边和顶点，便形成p_k+1的Voronoi多边形VR(P_k+1)，见图 2.

其中VD(P_k)表示k个水池点的Voronoi区域的并集, VR(P_i)表示水池点p_i的Voronoi区域.

2 Voronoi图的增量算法构造

采用增量法原理，利用Matlab编程语言作为辅助工具，构造平面离散水池点的Voronoi图.首先在Matlab界面中输入平面离散水池点集P={p₁, p₂, …, p_n}，进行一系列运算，输出P的Voronoi图VD(P).具体逻辑关系如下：

1) 如果n≤1，返回.

2) 计算VR(p₁)，VR(p₂).//通过计算p₁，p₂的中垂线得到.

3) for(k=3；k≤n；k=k+1){

4) 计算水池点p_k所属的Voronoi区域，设为VR(p_i).

5) 计算水池点p_k的Voronoi区域VR(p_k)}.

3 任意用水点的优化查询

利用点定位问题的分层方法，用O(n²)的时间和空间复杂度作预处理，花费O(logn)时间便可找出任意用水点所属的Voronoi区域^[16].

预处理过程：经过每一个顶点作一条水平线，对这些水平线的高度进行排序，记为数组A.Voronoi图平面被分成n+1个水平长条，称为层.整个Voronoi图被分成若干个梯形或三角形，且这些梯形或三角形要有指针指向原来所在的多边形.在每一层中，按位置在X₁方向上对以上梯形或三角形进行排序，把第i层排序结果记为数组B_i，具体见图 3.

给定任意用水点，先用二分法搜索数组A所在位置i层，再搜索B_i所在的梯形或三角形.数组A和B_i中的元素个数最多为O(n)，所以任意用水点的定位查询的时间复杂度为O(logn).

4 工程应用 4.1 工程背景

纳雍县地处中国贵州省西北部、毕节市南部.该县东西长56 km，南北宽48 km，总面积2 448 km².全县辖27个乡镇：雍熙街道办事处、居仁街道办事处、文昌街道办事处、王家寨镇、中岭镇、阳长镇、龙场镇、维新镇、乐治镇、百兴镇、张家湾镇，勺窝乡、新房乡、厍东关乡、董地乡、化作乡、寨乐乡、老凹坝乡、沙包乡、水东乡、曙光乡、姑开乡、锅圈岩乡、羊场乡、昆寨乡、左鶂嘎乡、猪场乡，具体见图 4.

4.2 根据山头地形确定水池位置

为满足纳雍县所有乡镇供水需求，分别在王家寨镇、中岭镇、阳长镇、维新镇、龙场镇、乐治镇、百兴镇、张家湾镇、羊场乡、锅圈岩乡、化作乡、猪场乡、沙包乡、水东乡设置供水水池.根据地形图，确定各乡镇山头高地作为水池修建地点，以重力进行输水.14个水池的位置坐标见表 1.

4.3 既定位置水池的区划Voronoi图

根据表 1中的数据将所定的14个水池移入直角坐标系X₁~X₂中，见图 5.

采用Matlab编程语言，首先输入点的坐标X=[11, 30, 12, 37, 21, 36, 54, 65, 59, 75, 36, 28, 54, 39]；Y=[76, 72, 58, 60, 44, 50, 50, 47, 39, 32, 28, 20, 17, 7]，然后定出点的位置，最后调用函数polytope得既定水池位置的重力式输水区划Voronoi图.

由图 6可知，每个Voronoi图多边形所包含的村庄到所在多边形内的水池取水路径最短，位于Voronoi图多边形边上的村庄到其两边水池取水路径一样长.

4.4 纳雍县供水范围区划图

结合图 4-6，得纳雍县供水范围区划图，见图 7.

4.5 任意用水点优化查询

为体现任意用水点的优化查询，以勺窝乡为例，给定用水点，采用任意用水点优化查询原理，将图 6分为18层，先用二分法搜索勺窝乡所在位置13层，再搜索B₁₃所在的梯形，最终确定勺窝乡居民在中岭镇所设的水池中取水路径最短，见图 8.

5 结论

针对山区供水布设既定水池位置的区划问题，提出一种基于Voronoi图的空间区域划分模型及算法，以纳雍县供水为研究实例，得出如下结论：

1) 在纳雍县的14个乡镇布置水池，利用Matlab编程得到既定水池区划Voronoi图，水池供水范围实现优化划分，山区农村供水得到保障且不会重叠覆盖.

2) 利用点定位问题的分层方法，用O(n²)的时间和空间复杂度作预处理，花费O(logn)时间便可定位任意用水点所属的Voronoi区域，大大提高任意用水点的优化查询.

3) 由于本文水池布设点较少，所以采用增量法的原理构造既定水池Voronoi图，算法简单且易于理解与编程.但增量算法总的时间复杂度为O(n²)，对于布设点较多的图形中，可用增量法先构造Voronoi图的对偶图——Delaunay三角化^[18]，再计算Voronoi图会简化时间的复杂度.