一种城市居住区划配套设施完备情况的评估模型
An Evaluation Model for How a City's Residential Area is Equiped with Some Non-trival Facilities
摘要： 本文提出了一种基于网络地图数据库，使用 K-means 算法与可视化技术的智慧城市评估模型。该模型可以用于评估城市居住区划的相关设施配套情况。我们具体地以南京市与苏州市为例，比较评价了两城市的居住区划的科教文化场所配备情况。
Abstrat: The essay introduced an evaluation model based on web map data base, K-means algorithm and visualization techniques, which can be quantitatively useful when inspecting how a city's residential area is equiped with some non-trival facilities. For example, we apply this method to compare Nanjing and Suzhou on how their residential area is equiped with science, education and culture places.
正文
本文的模型根据高德地图提供的兴趣点（Point of Interests, POI）数据库接口建立。我们首先采集城市的商务住宅位置信息，对这些信息进行聚类分析，就能大体上得到城市的实际居住区划模式。对原始数据进行聚类分析的意义在于，若城市某个地带有零星几处孤立的住宅公寓，其相关配套设施情况不会影响城市总体水平。而考虑这零星几处的住宅公寓再多一栋，再多两栋呢？这里就涉及到一个边界不清的量变与质变问题。因此我们采用了机器学习的手段进行聚类分析，就能有效地处理这一问题。
本文提出了一种基于网络地图数据库，使用 K-means 算法与可视化技术的智慧城市评估模型。该模型可以用于评估城市居住区划的相关设施配套情况。我们具体地以南京市与苏州市为例，比较评价了两城市的居住区划的科教文化场所配备情况。
接着，我们可以考虑科教文化场所等配套设施的情况。需要特别注意的是，这些配套设施未必只服务于某个特定的居住区，而恰恰可能有相当广的辐射范围。于是，本文主张用两个指标评估城市居住区总体的设施配套情况。第一个指标，是单独对设施点进行聚落分析，再考虑设施点聚落与居住区聚落的中心距离；第二个指标，是将设施点与居住区点数据糅合之后进行聚落分析，再考虑各个聚落的设施点占比。接下来本文将说明这两个指标的意义与合理性。
在实践上，我们首先通过 Python 语言解析数据库 JSON 文件，并将采集到的数据组织在一个向量结构中。又使用 C++ 语言实现了一个并发式计算的 K-means 算法（在效率上远优于其他语言提供的标准算法），用于城市功能区划的聚类分析与中心抽象。基于这个过程得到的结果，本文接着给出了直观的可视化表示方法，又进一步使用分析了上文提到的数值指标。
本文提供的模型以 K-means 算法为主体。抽象地讲，K-means 算法研究一个这样的问题：给定 (n) 个线性代数意义上的点 (P_1, dots, P_n)，按距离远近将这些点分类到 (k) 个不同的聚类中去。我们规定 (P_i, P_j) 两点间的距离为对应向量的模长 (||v||), (v = P_i - P_j)。
在本文的模型中，这里的点也就是高德地图提供的 POI，点的坐标就定义成数据库提供的经纬度信息。在这个意义下，一个首要的问题就是如何定义向量的模长。经纬度在本质上即地球所在的球坐标系中的前两坐标，因此我们简化地将向量模长定义为海平高度球面上的弧长。本文采用半正矢公式（Haversine Formula）进行计算：
首先引入半正矢函数 (hav( heta) = frac{1-cos heta}{2} = sin^2(frac{ heta}{2})) 。
半正矢公式指出：考虑将经纬度坐标分别为 ((varphi_1, heta_1), (varphi_2, heta_2)) 的两点，有
其中， (d) 即所求弧长； (R) 是球面高度，这里我们取海平面地球半径。
为了便于计算，我们可以将此式化为：
完成了这项工作，我们执行如下的 K-means 完成聚类分析：
首先在线性空间内随机地取一系列点 (mu_1, dots, mu_k) 作为数据点聚落的中心点，接下来：
-
对于每一个 (i)，计算 (k = arg min_j ||P_i - mu_j||)，并拟将点 (P_i) 归于 (mu_k) 所代表的聚落中去。
-
接着考虑更新每个聚落的中心点，具体地，对于每个聚类 (j) 中的点 (P_{v_1}, dots, P_{v_{L_j}})，其中 (L_j) 是该聚类中数据点的个数。我们取新的中心点为：
[mu_j = left( sum_{k=1}^{L_j} P_{v_k} ight) igg/ L_j ]在笛卡尔坐标系中，这一过程当然是求取各个点的算数平均的过程，也就是找到一个到所有点的距离均等的聚落“重心”。那么在球坐标系中，我们需要如何定义 (sum) 的符号含义，才能达到相同的效果呢？容易证明，我们只需简单的取 ((varphi_1 + varphi_2, heta_1 + heta_2) = (varphi_1, heta_1) + (varphi_2, heta_2)) 即可。
如果新的中心点在计算意义上与旧中心点重合，则算法结束，否则迭代执行 1 过程。具体地，我们规定距离小于 (10^{-12} m{km}) 时的两个新旧中心点视为重合点。
实现该算法时，向量模长的计算需要多次求三角函数值，开销较大。因此即使处理的数据点数量较小，也需要消耗大量的时间。因此我在具体实现时使用了并发计算方法。K-means 的执行过程中，每个数据点 (P_i) 的具体归类过程不会相关影响，因此可以考虑在几个进程中分别处理一些数据点。实现上创建四个进程 ( m{thd}_{1, 2, 3, 4}) 并用 ( m{thd}_i) 处理连续的点 (P_{frac{n}{4} imes (i - 1) + 1}, dots, P_{min{frac{n}{4} imes i, n}}) 进行运行加速。同时，我修改了数据写回的顺序，集中处理共享数据避免多次进程加锁的开销，又进一步提高了模型的效率。
完成这里的基本模型构建后，我们在高德地图开发者平台上申请一个 Web 服务应用，获得对应的 key 以访问数据库。接着查阅数据库提供的应用程序接口（Application Programming Interface, API）文档，可以编写相应程序进行数据采集。本文只采集了 POI 的名称与经纬度坐标信息，筛除了一些不需要的信息。我们特别关注了三种类型的 POI：也即 120000 商务住宅，140000 科教文化场所。其中商务住宅的聚类情况直接反映了城市的实际居住区划。另两种 POI 的分布情况则反映了具体类型的配套设施分布情况。
得到数据后，我们就可以进行模型运算，求取上文提到的两项指标。
第一项指标考虑设施点聚落与居住区聚落的中心距离，距离越近，聚落分布的一致性就越高。这一指标设计上是以住宅区为基点出发，反映我们与所关注的 POI 之间的距离。在结论上，这一指标人本地反映了城市居民获取对应资源的便利程度。我们称这项指标为「疏远度」。
第二个指标，是将设施点与居住区点数据糅合之后进行聚落分析，再考虑各个聚落的设施点占比。这一指标设计上是从配套设施出发，考虑其向周围辐射的能力。在结论上，这一指标物质地反映了城市资源设施向居住区覆盖的程度高低。我们称这项指标为「覆盖率」。
这里首先求疏远度，考虑设施点聚落与居住区聚落的中心距离。我们分别采集科教文场所与居住区聚落的信息，运行 K-means 算法后利用高德实验室提供的可视化接口，可以得到：
|
|
 |
 |
两组图片中，左侧均为科教文场所聚类分布情况，右侧均为住宅区域分布情况。
我们赋予不同的点互异的颜色，来代表它们从属的不同聚落。直观地看，我们可以发现两类 POI 的聚落分布是基本一致的。为了精准地进行量化评估，我们可以对数据做进一步的抽象。首先，考虑 K-means 算法得到的一系列聚落中心点 (mu_1, dots, mu_k)，根据聚落内点的个数确定其相对大小。我们将这些数据送由 Python 可视化处理，就能得到这样的图像：
点的意义如图例所示，其相对尺寸又反映了对应聚落的大小。我们记居住区聚落点为 (mu_1, dots, mu_k)，我们研究的 POI 聚落点为 (varphi_1, dots, varphi_k)，点 (p) 的大小记为 (|p|)。于是，我们则可以定义疏远度 (mathscr A) 为：
这里我们通过之前定义的向量模长表征两个聚落之间的距离，从而反映其关联程度，又以尺寸比值为权求距离的加权最小值，这样就能在相当程度上反映城市总体上居住区向特定 POI 的疏远度。编写程序进行运算，可以得到 (largemathscr{A}_{ ext{Suzhou}} = small 16.0698,~~~largemathscr{A}_{ ext{Nanjing}} = small 43.1623)。
这一数值对比提示：人本地看，苏州城市居民由城市居住区划出发，获取科教文资源更加方便。
接着我们计算覆盖率指标。首先将已经采集到的各类 POI 数据糅杂，重新输入 K-means 模型进行聚类分析。可以得到下面的表格：
| A | B | C | D | E | F | G | H | I | |
|---|---|---|---|---|---|---|---|---|---|
| Science, Education and Culture Places | 188 | 64 | 23 | 92 | 35 | 76 | 9 | 30 | 34 |
| Residential Area | 343 | 94 | 26 | 190 | 45 | 108 | 16 | 31 | 47 |
| Proportion | 0.54 | 0.68 | 0.88 | 0.48 | 0.78 | 0.70 | 0.56 | 0.97 | 0.72 |
| A | B | C | D | E | F | G | H | I | |
|---|---|---|---|---|---|---|---|---|---|
| Science, Education and Culture Places | 214 | 73 | 55 | 155 | 98 | 32 | 49 | 102 | 97 |
| Residential Area | 274 | 6 | 22 | 210 | 53 | 19 | 36 | 140 | 115 |
| Proportion | 0.78 | 12.16 | 2.50 | 0.73 | 1.84 | 1.68 | 1.36 | 0.72 | 0.84 |
在一些城市中的大学城、科创城等区域中，科教文场所与居住区域的比值会显著高于其他区域。这些区域聚类中的科教文场所辐射范围也有限，很难代表城市的总体水平，但也不能完全忽略。基于这样的考量，我们考虑将每个聚落抽象成点，将两类 POI 数量分别为横纵坐标，再通过线性回归拟合出一条过原点的直线。那么这条直线的斜率就反映了该城市 POI 覆盖的总体水平。我们不妨即定义覆盖率指标为该直线的斜率。
于是，记各聚落对应的点分别为 ((x_1,y_1), dots, (x_k, y_k))，通过最小二乘法计算，可以得到覆盖率满足公式：
编写程序，我们可以算出对应结果，并给出一个较好的可视化表示：
可视化表示之外，我们也可以在数值上发现 (largemathscr{C}_{ ext{Suzhou}} = small 0.840029,~~~largemathscr{C}_{ ext{Nanjing}} = small 0.560124)。这一数值对比提示，物质地看，苏州市的科教文化场所覆盖居住区划程度更高。
结论
本文基于高德地图的数据库，提供了一系列数学与计算机方法构建了一个评估模型。该模型可以用于评价城市居住区划的设施配套情况。模型基于机器学习的聚类分析方法，有效避免对城市边缘地区或特殊地段对结果的影响，也能较客观地反映城市的实际居住区划分野。
这一模型关注两个重要指标，疏远度 (mathscr{A}) 与 覆盖率 (mathscr{C})。其中，疏远度是一个更人本的指标，关注城市居民由居住区划出发获取对应设施资源的便利度；覆盖率是一个更物质的指标，关注城市特定设施覆盖居住区划的有效程度。两项指标都是我们评价城市居住区划配套设施完备情况的重要参考，不可偏废。
具体地，我们考虑城市居住区划的科教文设施配套情况，并分别在苏州市与南京市运行了模型，得到了一些初步的结果。我们发现疏远度计算结果为 (largemathscr{A}_{ ext{Suzhou}} = small 16.0698,~~~largemathscr{A}_{ ext{Nanjing}} = small 43.1623)，而覆盖率计算结果为(largemathscr{C}_{ ext{Suzhou}} = small 0.840029,~~~largemathscr{C}_{ ext{Nanjing}} = small 0.560124)。两项指标均显示苏州市的居住区划配备科教文设施情况更优。
总而言之，本文提供的模型与方法有较好的应用意义与可推广性：使用者只需修改一个或两个参数，就可以转而研究国内的任一城市的相关情况。本文通过一系列可视化方法，使得结果具有很好的直观性；又通过公式推导与程序计算，使得结果有较好的精确性。可以说，本文提出的模型与方法对于智慧城市规划，是很有一定作用和意义的。
附录：关键模型代码节选
# Python Code Blocks 1#: Cluster Visualization
plt.axis('equal')
ax = plt.subplot()
func = lambda x: min(x * 3, 500)
ax.scatter(eduX, eduY, s=list(map(func, eduSiz)), c = 'red',
alpha = 0.6, label = "Science, Education and Culture Places")
ax.scatter(resX, resY, s=list(map(func, resSiz)), c = 'green',
alpha = 0.6, label = "Residential Area")
plt.show()
// C++ Code Blocks 2#: Multithreading K-means Algorithm
double CalcuDis(const Point& a, const Point& b) {
static const double p = M_PI / 180;
return 12742 * asin(sqrt(0.5 - cos((b.lat - a.lat) * p) / 2
+ cos(a.lat * p) * cos(b.lat * p) * (1 - cos((b.lon - a.lon) * p)) / 2));
}
std::mutex mtx;
void runThrough(std::vector<int>& clusterSize,
std::vector<int>& assignment,
const std::vector<Point>& m_points,
const std::vector<Point>& m_centers,
int lef, int rig, int m_numCenters) {
std::vector<int> minVec;
for (int i = lef; i != rig; ++i) {
int minInd = 0;
double minDis = 0x7f7f7f7f, curDis;
for (short j = 0; j != m_numCenters; ++j) {
curDis = CalcuDis(m_points[i], m_centers[j]);
if (curDis < minDis) {
minInd = j;
minDis = curDis;
}
}
minVec.push_back(minInd);
}
std::lock_guard<std::mutex> lock(mtx);
for (int i = lef; i != rig; ++i) {
--clusterSize[assignment[i]];
assignment[i] = minVec[i - lef];
++clusterSize[minVec[i - lef]];
}
}
std::vector<index_t> Kmeans::Run(int maxIterations) {
std::vector<index_t> assignment(m_numPoints, 0);
int currIteration = 0;
const double eps = 1e-12;
double newX[m_numCenters], newY[m_numCenters];
int interval = m_numPoints / 4;
std::vector<int> clusterSize(m_numCenters, 0);
clusterSize[0] = m_numPoints;
while (currIteration < maxIterations) {
std::vector<std::thread> thdVec;
for (int ind = 0; ind < m_numPoints; ind += interval) {
thdVec.push_back(std::thread(runThrough,
std::ref(clusterSize), std::ref(assignment),
std::ref(m_points), std::ref(m_centers),
ind, std::min(m_numPoints, ind + interval), m_numCenters));
}
for (auto& thd: thdVec) thd.join();
for (short i = 0; i != m_numCenters; ++i) {
newX[i] = newY[i] = 0;
}
for (int i = 0; i != m_numPoints; ++i) {
newX[assignment[i]] += m_points[i].lat;
newY[assignment[i]] += m_points[i].lon;
}
bool flag = false;
for (short i = 0; i != m_numCenters; ++i) {
newX[i] /= clusterSize[i];
newY[i] /= clusterSize[i];
Point newCentroid(newX[i], newY[i], 0);
if (CalcuDis(newCentroid, m_centers[i]) > eps) {
m_centers[i] = newCentroid;
flag = true;
}
}
if (!flag) break;
++currIteration;
}
return assignment;
}
参考文献
[1] Robusto C C . Classroom Notes: The Cosine-Haversine Formula.[J]. Amer.math.monthly, 1957(1):38-40.
[2] Wong J A H A . Algorithm AS 136: A K-Means Clustering Algorithm[J]. Journal of the Royal Statistical Society, 1979, 28(1):100-108.
[3] Geladi P , Kowalski B R . Partial Least-Squares Regression: A Tutorial[J]. Analytica Chimica Acta, 1986, 185(1):1-17.