区域化是在给定一组约束的情况下将平面面单位分割成若干空间相邻区域的过程。区域化的研究可以看做基于费舍尔的均匀区域理论^[1]，均匀区域内部在一组属性上如经济、土地显示出高度的相似性。区域化问题广义而言可以划分为分区问题，而分区问题在土地利用分区优化上有广泛的研究应用^[2-6]。区域化问题也被广泛用于包括学区、保健区等各类研究区的划分^[7-8]。

区域化是一种特殊类型的聚类问题，受空间邻接约束^[9]。在实践中，研究人员一般具有关于区域的局部约束的领域知识，但经常较难提前确定区域数。Duque等基于p区域问题^[10]提出了最大p区域问题，在满足约束的情况下能够产生最大数量的区域^[11]。此后研究者基于p区域问题和最大p区域问题作出一系列扩展，包括考虑区域功能的p功能区模型^[12]，区域形状紧凑性区域模型，以及道路网络约束的最大p区域模型^[14]等。

区域化问题通常可以形成为整数规划问题^[10]。但混合整数规划问题计算密集^[11]，不适合中大型的分区问题。因此研究者一般使用启发式算法来对分区模型进行求解。常用算法包括模拟退火^[15]、自适应贪婪随机^[16]、禁忌搜索^[17]等。禁忌搜索算法由于其尽量避免局部最优解的能力被作为处理分区问题较好的解决方案^[11]，被研究者较多地融合在区域化问题的求解步骤中^[11]。

在实际的区域划分过程中，研究者往往希望划分的区域能够满足一定边界约束，以满足由于功能、权属问题所产生的区域不可跨越问题。直观地，这些划分的单元需尽量包含于所处边界，若有单元跨越多个边界的情况，则需要将这种不确定性的情况结合区域化问题的原有优化目标共同进行优化求解。另外，如果区域不能满足阈值约束而产生跨越边界的特殊情况，则将整个边界划分到形成的区域当中以保持划分的一致性。目前，国内外尚无针对最大p区域问题的边界约束相关研究，本文提出一种边界约束最大p区域问题，以方便研究者和实践者使用。针对该非确定性多项式(nondeterministic polynomially，NP)难题问题，本文基于禁忌搜索设计并实现了一种启发式算法，并使用模拟数据集和武汉市数据进行实验分析，以验证方法的有效性。

1 边界约束最大p区域问题 1.1 基础概念

边界约束最大p区域问题针对二维平面对象进行研究。基础概念包括空间单元、区域、单元属性^[11]。本文在模型定义中加入了边界，并扩展了邻接矩阵用于启发式计算。

边界由一组多边形Q={q₁, q₂…q_m}构成，Q完整覆盖整个区域S，并保证多边形之间没有间隙。与边界q_k相交的空间单元有两类，一类是完整覆盖的，即这些空间单元只与q_k一个边界相交；另一类是与空间单元内部相交的，但空间单元整体可与多个边界相交。

空间邻接矩阵用于预先存储空间单元之间，以及空间单元与内部、边界之间的邻接关系。其中空间单元矩阵W为一个n×n的二值对称矩阵，本文以常见的皇后(Queen)邻接矩阵形式进行计算^[18]；空间单元与边界的邻接关系矩阵E为一个n×m的二值矩阵，其中e_ij代表空间单元g_i内部是否与边界q_i相交；边界之间的邻接关系矩阵F为一个m×m的二值矩阵，用于启发式算法中边界之间关系的判断。

区域化的分区结果P是一组区域的集合。区域个数p=|P|，边界约束最大p区域问题需要在计算得出满足阈值约束的前提下，最大化区域个数p并最小化顾及空间单元从属边界的不确定度加权差异性的空间划分P_opt。

1.2 模型定义

本文模型主要在最大p区域问题中加入了边界约束，为此扩展了最大p区域问题的模型参数、决策变量、目标函数和约束条件。

1.2.1 模型参数

1) G为空间单元集合, G= {1…n};

2) i、j为单元索引, i, j∈G;

3) k为潜在区域索引, k= {1…n};

4) c为空间单元赋值顺序, c= {0…h}, h= {n-1 };

5) ${w_{ij}} = \left\{ \begin{array}{l} 1, 空间单元i、j邻接, i, j \in G;i \ne j;\\ 0, 不邻接 \end{array} \right. $

6) $ N_{i}=\left\{j | w_{i j}=1\right\}$, 空间单元i邻居;

7) d_ij为空间单元i、j差异性, i，j∈G, i < j;

8) t_i为用于阈值计算的单元属性, i∈G; T为区域内部单元属性t求和所需达到的最小阈值;

9) $o=1+\log \left(\sum\limits_{i} \sum\limits_{j | j>i} d_{i j}\right) $, 代表$ {\sum\limits_i {\sum\limits_{j|j > i} {{d_{ij}}} } }$向下取整的位数, i, j∈I;

10) Q为边界集合, Q= {1…m};

11) q为边界索引, q∈G;

12) $ {e_{iq}} = \left\{ \begin{array}{l} 1, 空间单元i与边界q相交\\ 0, 不相交i \in G, q \in Q \end{array} \right.$

1.2.2 决策变量

1) ${s_{ij}} = \left\{ \begin{array}{l} 1, 单元i、j同属于一个区域k, i < j;\\ 0, 其他 \end{array} \right. $

2) $x_i^{kc} = \left\{ \begin{array}{l} 1, 单元i是属于k的c阶单元;\\ 0, 其他 \end{array} \right. $

3) $ y_k^q = \left\{ \begin{array}{l} 1, 区域k属于边界q\\ 0, 其他 \end{array} \right.。$

1.2.3 目标函数

最小化：

$ \begin{aligned} z &=\left(-\sum\limits_{k=1}^{n} \sum\limits_{i=1}^{n} x_{i}^{k 0}\right) \times 10^o+\\ & \sum\limits_{i} \sum\limits_{j | j>i} d_{i j} \cdot s_{i j} \cdot f(i, j) \end{aligned} $

(1)

其中，

$ f(i, j)=\sum\limits_{k=1}^{n} I(i, j, k) \cdot\left(1+\lambda \cdot u_{i k} \cdot u_{j k}\right) $

(2)

$ I(i, j, k)=\sum\limits_{c=0}^{h} x_{i}^{k c} \cdot \sum\limits_{c=0}^{h} x_{j}^{k c} $

(3)

$ u_{i k}=1-\sum\limits_{q=1}^{m} y_{k}^{q} \cdot e_{i q} \cdot \frac{a(i, q)}{a(i)} $

(4)

$ \lambda \in[0, +\infty) ; u_{i k}, u_{j k} \in[0, 1) $

(5)

1.2.4 约束条件

$ \sum\limits_{i=1}^{n} x_{i}^{k 0} \leqslant 1, \forall k=1 \cdots n $

(6)

$ \sum\limits_{k=1}^{n} \sum\limits_{c=0}^{h} x_{i}^{k}=1, \forall i=1 \cdots n $

(7)

$ \begin{array}{l} x_i^{kc} \le \sum\limits_{j \in {N_i}} {x_j^{k(c - 1)}, } \forall i = 1 \cdots n;\\ \forall k = 1 \cdots n;\forall c = 1 \cdots h \end{array} $

(8)

$ \sum\limits_{i=1}^{n} \sum\limits_{c=0}^{h} x_{i}^{k} t_{i} \geqslant T \cdot \sum\limits_{i=1}^{n} x_{i}^{k 0}, \forall k=1 \cdots n $

(9)

$ \begin{aligned} s_{i j} & \geqslant \sum\limits_{c=0}^{h} x_{i}^{k}+\sum\limits_{c=0}^{h} x_{j}^{k c}-1 \\ \forall i, j &=1 \cdots n | i <j ; \forall k=1 \cdots n \end{aligned} $

(10)

$ x_{i}^{k_{c}} \in\{0, 1\}, \forall i=1 \cdots n ; \forall k=1 \cdots n ; \forall c=0 \cdots h $

(11)

$ s_{i j} \in\{0, 1\}, \forall i, j=1 \cdots n | i <j $

(12)

$ \sum\limits_{q=1}^{m} e_{i q} \cdot y_{k}^{q} \geqslant \sum\limits_{c=0}^{h} x_{i}^{k c}, \forall i=1 \cdots n ; \forall k=1 \cdots n $

(13)

$ \begin{array}{l} \left( {\sum\limits_{q = 1}^m {y_k^q} - 1} \right) \cdot \left( {\sum\limits_{c = 0}^h {x_i^{kc}} - {e_{iq}} \cdot y_k^q} \right) \ge 0, \\ \forall i = 1 \cdots n;\forall k = 1 \cdots n;\forall q = 1 \cdots m \end{array} $

(14)

$ y_{k}^{q} \in\{0, 1\}, \forall k=1 \cdots n ; \forall q=1 \cdots m $

(15)

$ \sum\limits_{q=1}^{m} y_{k}^{q} \geqslant 1, \forall k=1 \cdots n $

(16)

其中，模型参数中相比最大p区域问题，本文增加了参数10)~12)，用以表示边界及其与空间单元邻接关系。决策变量中添加变量3)用以表示区域k从属的边界q。一般情况下，区域k只从属于一个边界q。在出现阈值条件不满足的情况下，区域k可能从属多个边界。

边界约束最大p区域问题的目标函数中第一部分$ \left(-\sum\limits_{k=1}^{n} \sum\limits_{i=1}^{n} x_{i}^{k 0}\right) \times 10^{o}$来源于原始最大p区域问题，其保证了区域个数最大化优先考虑。而本文在目标函数第二部分则提出融入空间单元从属边界的不确定度函数f(i，j)对单元差异性进行加权，其物理意义在于对于空间单元与多个边界相交时, 使目标函数能够综合考虑单元与不同区域和边界的从属关系。f(i，j)代表了单元i和单元j从属某同一区域k的整体不确定度，其定义和最大p区域问题中的d_ij·s_ij对应。为符合混合整数规划的严格形式，因此f(i，j)函数中需基于决策变量来定义。f(i，j)函数输入为两个单元i、j，在其组成中，I(i，j，k)为指示函数，代表单元i、j是否同属某一区域k。λ参数为不确定度权值，用以控制不确定度加权的影响程度，其物理含义对应于研究者和实践者实际应用需要。举例而言，当i、j不同属任何区域时，f (i，j) ≡0；当i、j同属一个区域k时，若λ=0，则f(i，j) ≡1，否则f (i，j) ≥1，此时目标函数会驱使单元从属于不确定度较低的边界。例如，λ越大，则f(i，j)越大，当λ=0时，则f(i，j)=1，即不对目标函数产生任何影响。u_ik和u_jk代表单元i、j在区域k内的不确定度，以u_ik为例，形式如式(4)所示，其考虑一个区域包括多个边界的情况。本文使用相交面积比例这一相对简化的度量来代表不确定度，其基本假设是单元内部的均匀性，a(i，q)和a(i)分别代表单元i和边界q的相交面积和单元i的面积。若单元i与区域从属边界的相交面积越大，则u_ik越小；若相交面积为0，则u_ik=0。受后续约束限制，单元一定与某一区域从属的边界有交集，因此u_ik小于1，而不会等于1。

约束条件式(6)~(12)代表最大p区域问题中的原有约束条件，其中式(6)~(8)采用顺序赋值方法^[19]保证区域中的空间单元邻接性，c代表一个单元划分给区域的顺序。式(6)代表一个区域k只有一个根单元(即c=0)，式(7)代表一个单元仅能被赋予给某一个区域k并且只有一个顺序值c，式(8)保证了区域单元之间的空间邻接性，即除根单元外，一个单元只有在其任意邻接单元被划归给某区域k后才能被划归给区域k。式(9)保证区域满足阈值限制，式(10)代表s_ij所反映的两个单元同属一个区域的限制。式(11)和(12)保证决策变量x_i^kc和s_ij的值域。

本文在最大p区域问题基础上进一步提出条件式(13)~(16)，代表边界约束条件。其中，式(13)保证若区域k属于边界q，则区域k内的空间单元也必在q之内。式(14)保证若区域k从属于多个边界，则这些边界内的单元必都属于区域k。式(15)和式(16)保证决策变量y_k^q的值域并确保每个区域最少从属于一个边界。

2 启发式算法

区域化问题启发式算法一般包括两个主要步骤：初始解生成以及局部搜索。

2.1 初始解生成

初始解生成主要包括以下6个步骤。

1) 阈值过滤。遍历所有单元，若单元变量t超过阈值T，则将该单元单独构成一个区域R加入初始解P。

2) 区域综合。遍历所有边界，若有边界内所有空间单元集合G_q对于阈值变量t的和小于T，则将该边界内的所有单元合成为一个特殊单元g_q，称为边界单元，加入边界单元集合O。从整个空间单元集合G中删除G_q。

3) 随机区域生长。从剩余没有赋值的空间单元S中，随机选择单元作为种子单元生长区域R。该贪婪选择过程将待选单元与区域内已有单元分别求取加权差异值d_ij·f (i，j)并加总，最终选择最小化该值的邻接单元，并在选择过程中始终保证R中的空间单元包含在同一边界内。区域生长过程中一直需要维持并更新其相邻单元列表W，以加速单元选择的过程。当达到阈值条件时，停止R生长，将其添加入P。若生长过程中W为空(即没有找到邻接的单元可供选择)，则删除R。最终剩余未分配的空间单元为飞地单元。

4) 飞地赋值。遍历所有飞地单元，计算其与P中所有邻接区域的加权边界单元交界不确定度和差异性值，并加入最小化该加权值的区域。

5) 边界单元合并。遍历步骤2)中形成的边界单元集合：将边界单元与周边区域或边界单元进行贪婪形式的合并形成区域。优先选择其他边界单元，若没有，则选择在第3)步中划分区域较少的边界；若区域相等，则选择加权边界单元交界不确定度和差异性较小的边界。添加后仍达不到阈值则继续选择，直到满足阈值条件为止。最终形成的区域加入P。

6) 区域数过滤。在反复上述步骤1)~5)生成初始解集合M后，得到最大的区域个数值p，并获取集合I的子集M′={P|numRegion (P) = p, P∈M}。

2.2 局部搜索

局部搜索的整体框架与求解最大p区域问题及其他扩展形式类似^[14]。其对§2.1获取的初始解集合，其中每个初始解进行局部搜索优化，最终解为优化解集合中目标函数最小的解。每一次移动对应两个区域中一对空间单元的交换。区别于最大p区域问题，对于仅与一个边界相交的空间单元(即边界内部的单元)，交换仅在各边界内的区域集合里进行。而与多个边界相交的空间单元则可在相交的多个边界内的区域集合里进行交换。由独立边界单元形成的区域不加入优化过程。其余每个边界q内部维持一个单元和区域的列表G_q，G_q包含所有不违背空间连接性的潜在移动。局部搜索将随机从G_q中选择一个移动, 并由此将原解P产生新解P′。采用禁忌搜索允许P′的目标函数值大于P，即通过允许局部次优解的形成来尽量避免局部最优解，以此更为全面地探索解空间。

3 实验分析

本节使用模拟数据和实际数据对方法有效性进行验证。算法使用Java实现。程序运行配置为：中央处理器为i5英特尔CPU(central processing unit)，内存为16GB DDR3。本文侧重对边界约束最大p区域问题的性质和划分结果进行分析，因此在局部搜索阶段统一使用禁忌长度为85，并设定初始解生成数量等于5 000。

3.1 模拟实验

本文采取最大p区域模型所采用的模拟数据进行有效性验证，空间单元数目为1 024。其中，差异属性为一个属性，由参数0.9的空间自回归过程模拟生成，阈值属性由一个[10, 15]区间的均匀分布产生。本文进一步设计不规则的模拟边界对算法进行验证(m=12)。图 1(a)显示了模拟边界叠合空间单元效果。下面分别以阈值参数T和不确定度权值参数λ对模型结果进行分析。图 1(b)显示了最大p区域问题的求解结果(阈值T=100)。

图 2显示了不同阈值T下的边界约束最大p区域问题的分区结果(取λ=0)。图中不同颜色代表不同归属边界, 黑色加重边框的区域代表其有多个从属边界。相比原始最大p区域问题，边界约束的模型中区域生长过程明显受边界制约，在阈值增至100时，则发生了边界合并的情况。这是由于被合并的边界内部单元阈值属性加总达不到100，只有合并相邻的边界才能使得单元阈值属性加总超过100，产生满足阈值约束的区域。

图 3显示了不同不确定度权值λ下的边界约束最大p区域问题求解比较(取T=100)。通过与边界关系可见，λ值越大，与多边界交叉的空间单元越倾向于属于与其交叉面积大边界包含的区域，这一倾向过程同时受到目标函数中差异值最小化的制约。当λ值越大，制约效果则相对越小。

3.2 武汉市数据实验

实际数据采用武汉市公里格网数据进行分析，格网数据为面域数据，除去部分水域外，每个格网代表一个1 km²的方格，每个格网有人口数据等若干属性。边界则使用武汉区一级行政区数据。实验使用格网的人口数作为阈值属性，城市用地面积作为差异属性，分区结果可用于对与土地利用相关的社会学应用提供初始的分区。在实际应用中，需要结合应用特点来选择合适的阈值属性和多个差异属性。本文设定人口阈值为15万，设定初始解生成数量为100，不确定度权值为10。启发式计算的划分结果如图 4所示。表 1列出了计算结果的目标函数、区域个数及各阶段消耗时间。

可以看出，模型根据边界限制在各区内进行划分，在人口密集的地方如洪山等主城区划分的区域则较多较密，而人口较少的地方则产生较大的区域，人口数达不到阈值的区(汉南区)则与邻近区发生了合并。最大p区域问题的初衷在于针对研究者无法确定分区个数时，提供一种通过设置阈值条件，根据区域差异值最小化目标得到分区结果的自动化方法，区域个数由最大化区域数和减少差异值目标下结合阈值的设置而产生，这虽不能满足所有研究情况，但提供了一种有效自动确定分区数目的方法，其产生的区域个数合理性需要结合应用领域进一步实证分析。本文扩展的边界约束最大p区域模型则进一步加入边界约束及融合空间单元与边界交叉产生的不确定度。实验结果表明本文算法能够适用于真实数据集，对于空间分析和实践业务所需要的研究或实验区域划分有实际意义。实验结果也表明，对于此类真实数据，启发式算法的执行特别是局部搜索阶段的效率仍有待提升，未来需要与其他局部搜索算法对比分析，并研究集群环境下的分布式并行算法，从而更可能地接近最优解^[20]。

4 结语

本文考虑了最大p区域问题在现实运用中的一个重要限制条件，即边界约束，并同时保证阈值约束和融入空间单元从属边界的不确定度的加权差异性最小化目标。在此基础上提出了边界约束最大p区域问题，并给出其模型形式化定义及其启发式算法。这为研究者和业务部门提供了相较最大p区域问题更为灵活的控制方式，在设计研究区域，如人口、犯罪、交通及各类研究区域等提供了一种定量工具。在考虑边界上，仍有许多可以扩展之处，如研究者有时需要考虑区域之间的物理隔段，包括河流、围墙、高速公路等，使得研究区域不能连接，因此在模型约束条件中需要考虑隔段周围如何分隔区域。