一种城市居住区划配套设施完备情况的评估模型
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(\theta) = \frac{1-\cos\theta}{2} = sin^2(\frac{\theta}{2})\) 。

​ 半正矢公式指出：考虑将经纬度坐标分别为 \((\varphi_1,\theta_1), (\varphi_2,\theta_2)\) 的两点，有

\[hav(\frac{d}{R}) = hav(\varphi_2-\varphi_1)+\cos\varphi_1 \cos\varphi_2\cdot hav(\theta_2-\theta_1) \]

​ 其中， \(d\) 即所求弧长； \(R\) 是球面高度，这里我们取海平面地球半径。

​ 为了便于计算，我们可以将此式化为：

\[d = 2R\arcsin(\sqrt{\sin^2{(\frac{\varphi_2-\varphi_1}{2})}} + \cos\varphi_1 \cos\varphi_2\cdot \sin^2(\frac{\theta_2-\theta_1}{2})) \]

​ 完成了这项工作，我们执行如下的 K-means 完成聚类分析：

​ 首先在线性空间内随机地取一系列点 \(\mu_1, \dots, \mu_k\) 作为数据点聚落的中心点，接下来：

1. 对于每一个 \(i\)，计算 \(k = \arg \min_j ||P_i - \mu_j||\)，并拟将点 \(P_i\) 归于 \(\mu_k\) 所代表的聚落中去。

2. 接着考虑更新每个聚落的中心点，具体地，对于每个聚类 \(j\) 中的点 \(P_{v_1}, \dots, P_{v_{L_j}}\)，其中 \(L_j\) 是该聚类中数据点的个数。我们取新的中心点为：

   \[\mu_j = \left( \sum_{k=1}^{L_j} P_{v_k} \right) \bigg/ L_j \]

   在笛卡尔坐标系中，这一过程当然是求取各个点的算数平均的过程，也就是找到一个到所有点的距离均等的聚落“重心”。那么在球坐标系中，我们需要如何定义 \(\sum\) 的符号含义，才能达到相同的效果呢？容易证明，我们只需简单的取 \((\varphi_1 + \varphi_2, \theta_1 + \theta_2) = (\varphi_1, \theta_1) + (\varphi_2, \theta_2)\) 即可。

​ 如果新的中心点在计算意义上与旧中心点重合，则算法结束，否则迭代执行 1 过程。具体地，我们规定距离小于 \(10^{-12} \rm{km}\) 时的两个新旧中心点视为重合点。

​ 实现该算法时，向量模长的计算需要多次求三角函数值，开销较大。因此即使处理的数据点数量较小，也需要消耗大量的时间。因此我在具体实现时使用了并发计算方法。K-means 的执行过程中，每个数据点 \(P_i\) 的具体归类过程不会相关影响，因此可以考虑在几个进程中分别处理一些数据点。实现上创建四个进程 \(\rm{thd}_{1, 2, 3, 4}\) 并用 \(\rm{thd}_i\) 处理连续的点 \(P_{\frac{n}{4} \times (i - 1) + 1}, \dots, P_{\min\{\frac{n}{4} \times i, n\}}\) 进行运行加速。同时，我修改了数据写回的顺序，集中处理共享数据避免多次进程加锁的开销，又进一步提高了模型的效率。

​ 完成这里的基本模型构建后，我们在高德地图开发者平台上申请一个 Web 服务应用，获得对应的 key 以访问数据库。接着查阅数据库提供的应用程序接口（Application Programming Interface, API）文档，可以编写相应程序进行数据采集。本文只采集了 POI 的名称与经纬度坐标信息，筛除了一些不需要的信息。我们特别关注了三种类型的 POI：也即 120000 商务住宅，140000 科教文化场所。其中商务住宅的聚类情况直接反映了城市的实际居住区划。另两种 POI 的分布情况则反映了具体类型的配套设施分布情况。

​ 得到数据后，我们就可以进行模型运算，求取上文提到的两项指标。

​ 第一项指标考虑设施点聚落与居住区聚落的中心距离，距离越近，聚落分布的一致性就越高。这一指标设计上是以住宅区为基点出发，反映我们与所关注的 POI 之间的距离。在结论上，这一指标人本地反映了城市居民获取对应资源的便利程度。我们称这项指标为「疏远度」。

​ 第二个指标，是将设施点与居住区点数据糅合之后进行聚落分析，再考虑各个聚落的设施点占比。这一指标设计上是从配套设施出发，考虑其向周围辐射的能力。在结论上，这一指标物质地反映了城市资源设施向居住区覆盖的程度高低。我们称这项指标为「覆盖率」。

​ 这里首先求疏远度，考虑设施点聚落与居住区聚落的中心距离。我们分别采集科教文场所与居住区聚落的信息，运行 K-means 算法后利用高德实验室提供的可视化接口，可以得到：

图1: 南京 POI 聚落分布 图2：苏州 POI 聚落分布

​ 两组图片中，左侧均为科教文场所聚类分布情况，右侧均为住宅区域分布情况。

​ 我们赋予不同的点互异的颜色，来代表它们从属的不同聚落。直观地看，我们可以发现两类 POI 的聚落分布是基本一致的。为了精准地进行量化评估，我们可以对数据做进一步的抽象。首先，考虑 K-means 算法得到的一系列聚落中心点 \(\mu_1, \dots, \mu_k\)，根据聚落内点的个数确定其相对大小。我们将这些数据送由 Python 可视化处理，就能得到这样的图像：

图3：城市 POI 聚落叠加分布图，左为南京，右为苏州

​ 点的意义如图例所示，其相对尺寸又反映了对应聚落的大小。我们记居住区聚落点为 \(\mu_1, \dots, \mu_k\)，我们研究的 POI 聚落点为 \(\varphi_1, \dots, \varphi_k\)，点 \(p\) 的大小记为 \(|p|\)。于是，我们则可以定义疏远度 \(\mathscr A\) 为：

\[\mathscr{A} = \sum_{i = 1}^k \min_{1 \leq j \leq k} ( ||\mu_i - \varphi_j|| \cdot \frac{|\mu_i|}{|\varphi_j|}) \]

​ 这里我们通过之前定义的向量模长表征两个聚落之间的距离，从而反映其关联程度，又以尺寸比值为权求距离的加权最小值，这样就能在相当程度上反映城市总体上居住区向特定 POI 的疏远度。编写程序进行运算，可以得到 \(\large\mathscr{A}_{\text{Suzhou}} = \small 16.0698,~~~\large\mathscr{A}_{\text{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

表1：南京总体 POI 聚落点构成表

	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

表2：苏州总体 POI 聚落点构成表

​ 在一些城市中的大学城、科创城等区域中，科教文场所与居住区域的比值会显著高于其他区域。这些区域聚类中的科教文场所辐射范围也有限，很难代表城市的总体水平，但也不能完全忽略。基于这样的考量，我们考虑将每个聚落抽象成点，将两类 POI 数量分别为横纵坐标，再通过线性回归拟合出一条过原点的直线。那么这条直线的斜率就反映了该城市 POI 覆盖的总体水平。我们不妨即定义覆盖率指标为该直线的斜率。

​ 于是，记各聚落对应的点分别为 \((x_1,y_1), \dots, (x_k, y_k)\)，通过最小二乘法计算，可以得到覆盖率满足公式：

\[\mathscr{C} = \left({\sum_{i=1}^k x_i \cdot y_i}\right) \bigg/ \left({\sum_{i=1}^k x_i^2}\right) \]

​ 编写程序，我们可以算出对应结果，并给出一个较好的可视化表示：

图4：城市 POI 覆盖率回归图，左为苏州，右为南京

​ 可视化表示之外，我们也可以在数值上发现 \(\large\mathscr{C}_{\text{Suzhou}} = \small 0.840029,~~~\large\mathscr{C}_{\text{Nanjing}} = \small 0.560124\)。这一数值对比提示，物质地看，苏州市的科教文化场所覆盖居住区划程度更高。

结论

​ 本文基于高德地图的数据库，提供了一系列数学与计算机方法构建了一个评估模型。该模型可以用于评价城市居住区划的设施配套情况。模型基于机器学习的聚类分析方法，有效避免对城市边缘地区或特殊地段对结果的影响，也能较客观地反映城市的实际居住区划分野。

​ 这一模型关注两个重要指标，疏远度 \(\mathscr{A}\) 与 覆盖率 \(\mathscr{C}\)。其中，疏远度是一个更人本的指标，关注城市居民由居住区划出发获取对应设施资源的便利度；覆盖率是一个更物质的指标，关注城市特定设施覆盖居住区划的有效程度。两项指标都是我们评价城市居住区划配套设施完备情况的重要参考，不可偏废。

​ 具体地，我们考虑城市居住区划的科教文设施配套情况，并分别在苏州市与南京市运行了模型，得到了一些初步的结果。我们发现疏远度计算结果为 \(\large\mathscr{A}_{\text{Suzhou}} = \small 16.0698,~~~\large\mathscr{A}_{\text{Nanjing}} = \small 43.1623\)，而覆盖率计算结果为\(\large\mathscr{C}_{\text{Suzhou}} = \small 0.840029,~~~\large\mathscr{C}_{\text{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.