【GIS实践】基于改进区域生长的临近斑块分类方法

基于改进区域生长的临近斑块分类方法

---

一、背景描述

根据Tobler第一定理，距离相近的地物有着更加紧密的关系，因此，在对城市低效用地的重新评估过程中，往往需要考虑当前低效用地斑块与周边斑块的关系。低效用地斑块的整改也会对周边非低效用地产生影响，为了降低二次规划整治成本，提升区块景观美学价值，最大化经济社会效益，本文将考虑低效用地斑块与周边地物的关联程度，并将满足条件的非低效用地与低效用地聚合。

现有数据中，斑块用地是大小不一、形状各异的矢量面数据，传统的空间分类方式较难考虑到矢量面的大小、拓扑关系，容易将小低效用地与大面积非低效用地联合考虑，这样的结果并不符合最大化社会经济效益的目标。为了基于以低效斑块用地为核心的拓扑关系进行区域分类，本文拟采用区域生长法来实现相似斑块的合并。

---

二、算法理论

区域生长算法是图像处理领域中常见的图像分割方法，其基本思想是从图像的一个像素或邻域中出发，比较相邻部分的特征，按照相似原则不断合并，直至不能再合并的分割方式。区域生长算法按照邻域和相似性准则不同，可分为单一型、质心型和混合型三种区域增长方法，其核心任务在于：

• 确定区域的数目
• 选择有意义的特征
• 确定相似性准则

由于斑块之间是面积不一、没有规范空间分布的矢量面数据，在本文的工作中，将采用混合型增长的方式进行合并。

2.1 特征权重确定

为了确保选取的特征最具代表性，可以考虑基于信息增益的方式，利用熵权法(Entropy Weight Method)对现有特征进行权重分配。熵权法是客观赋权的一种方式，对应的主观赋权有专家打分法，熵权法只注重数据本身，不会受到人为因素的影响。熵权法是利用信息稳定程度而提出的方法，一般来说，某列属性越稳定，它的信息就越可信，那么在实际的权重也应当越高。

假设当前每个斑块具有$n$个特征，并且一共有$m$个斑块，其中低效用地斑块为$k$个，非低效用地为$m-k$个。

特征$i$的信息熵写作：

$$E_i=\sum-\log(P(X_{ij}))X_{ij}$$

其中，$X_{ij}$表示特征$i$的随机变量$j$。

为了去除各个特征的量纲影响，需要先对特征进行归一化。

$$Y_{ij}=\frac{X_{ij}-min(X_i)}{max(X_i)-min(X_i)}$$

其中，$Y_{ij}$表示特征$i$的随机变量$j$的归一化值。接着，基于频率替代概率，并计算归一化信息熵：

$$p_{ij}=\frac{Y_{ij}}{\sum Y_{ij}} \\ H_i=-\frac{1}{ln(m)}\sum p_{ij}ln(p_{ij})$$

其中，$H_i$表示特征$i$的归一化信息熵，最后再通过计算各个特征的归一化信息熵占比来确定特征的权重：

$$W_i=\frac{1-H_i}{k-\sum H_j}$$

经过选取后的特征权重为$W_i, i\in [1,n]$。

---

2.2 相似性准则确定

2.2.1 缓冲区距离

此时，我们需要考虑如何设置相似性准则，并且实现合并。

区块的合并需要考虑空间可达性，一般认为，设施==半径D[查文献]==(步行距离1200/800/400)内的区块具有良好的可达性。

因此，我们将对低效用地构建半径为$D$的缓冲区，并计算该区域内的相邻地块特征。

改进： 是否可以采用自适应缓冲区？

2.2.2 相似度度量

分类结果的具体质量应该取决于相似度度量算法，以及它的实现。如果一个分类还能挖掘一些或者全部的隐藏模式，我们也可以认为它是一个高质量的分类方法。通常，数据科学家们使用相异度/相似度矩阵 （Dissimilarity/Similarity metric）来对不同数据特征之间的相似度进行度量，一般是通过距离公式进行。

在我们的工作中，将采用余弦相似度对数据进行度量：

$$sim(x,y)=\frac{x\cdot y}{||x||\times||y||}$$

其中，$x,y$为斑块$X$与斑块$Y$的特征向量，这些特征向量的构成方式为：

$$x=[W_1\times X_1,...,W_i\times X_i]^T$$

本文将采用非负实数$\epsilon$作为相似度阈值，当$sim(x,y)\ge \epsilon$时，视为二者具有较高的相似度。本文将分别设置$\epsilon为0.6,0.7,0.8$，并探究其对最终结果的影响。、

同时，为了考虑空间衰减，可能需要根据缓冲区的类型来进行空间衰减，距离越接近的区块，$\epsilon$的取值越低。

2.2.3 面积关系

实际上，不同的区块有着一个至关重要的属性：面积。即使两个斑块的特征十分接近，但若低效用地面积远小于非低效用地，那么对于他们之间的合并是有待商榷的。为了尽可能保证经济效益，应当避免将面积过大的非低效用地考虑进面积较小的低效用地的联合整改队列中。

考虑到低效用地$i$与非低效用地$j$的面积$S_i$与$S_j$，现有如下分段措施：

$$\begin{cases} \epsilon=1,\quad S_j\ge2S_i \\\epsilon=\frac{1}{1+e^{-\alpha\frac{S_j}{S_i}}},\quad S_j \in (0.5S_i,2S_i) \\\epsilon=0.5,\quad S_j\le0.5S_i \end{cases}$$

其中，$\alpha$为一个非负参数，用于控制$\epsilon$的变化速率。本实验将$\alpha$的值设置为$3$。

2.2.4 空间重叠

实际上，由于缓冲区本身的特性，有些非低效用地可能会被多个低效用地斑块所考虑，此时，我们将多个低效用地视为一个面积更大的区域，该区域的面积为：

$$S=\sum S_i$$

其中，$S$表示大区域的面积，$S_i$表示可达该非低效用地的第$i$个低效斑块。

实际上 ，当发生空间叠置时，即使是相同类型的低效斑块，其属性也有可能有较大不同。

因此，在这种情况下，需要先评估重叠的几个低效斑块的相似度，再决定能否合成一个大的区块。

若低效斑块相似度较高，则将该非低效用地与几个低效斑块用地进行相似度比较，选择最高值来判断是否高于相似度阈值。

若低效相似度较低，则需要分别考虑每个低效斑块与非低效斑块的相似度，而不是将其视为一个更大的源区块。

---

三、技术流程图

graph TB
	id1[所有斑块属性导出]-->id2(基于熵权法确定特征权重)
	id2-->id3(构建特征向量)
	
	id4[低效用地]-->id5(构建缓冲区)
	id5-->id6(确定缓冲区内的非低效用地)
	
	id6-->id7(相似度计算)
	id3-->id7
	
	id7-->id8(判断空间关系)
	id8-->id10(确定低效斑块集群面积)
	id10-->id11(更新阈值参数)
	id11-->id12{相似度是否大于阈值}
	id12--Y -->id13[将非低效用地视为联合整改部分]
	id12--N-->id14[不做修改]

---

四、实现代码

4.1 数据预处理

import pandas as pd
import numpy as np
import math
import os
import copy
import geopandas as gpd
import random
import collections

path=r"G:\bankuai"
data=gpd.read_file(path+r"\zonghe2.shp")

# visualizing the distribution of blocks in ChangSha 
data.plot()

# exploratory data analysis
D=copy.deepcopy(data)
data.shape
# (3700, 14)
data.head(5)

	OBJECTID	lb	Shape_Leng	Shape_Area	idd	OBJECTID_1	idd_1	COUNT	AREA	MEAN	MEAN_1	MEAN_12	MEAN_12_13	geometry
0	1	0	462.698517	7842.037150	1	1	1	634	10138.566790	4.885384	419.828282	583860.130028	1.097882e+06	POLYGON Z ((694423.228 3118300.460 0.000, 6943...
1	2	0	1693.463425	115652.876853	2	2	2	9333	149248.018695	2.210137	49.167335	177414.284438	7.246673e+05	POLYGON Z ((698632.417 3106524.626 0.000, 6987...
2	3	11	1100.566517	76822.340277	3	3	3	6191	99002.944792	2.194382	55.816331	83970.746248	6.581549e+05	POLYGON Z ((698213.537 3106892.204 0.000, 6981...
3	4	0	813.027834	33316.074792	4	4	4	2683	42905.007410	3.220493	60.284221	138651.257203	7.675609e+05	POLYGON Z ((700789.640 3106976.923 0.000, 7006...
4	5	0	931.012175	54895.664372	5	5	5	4426	70778.070368	1.702563	34.498813	145172.013136	7.156062e+05	POLYGON Z ((699126.327 3107255.440 0.000, 6989...

这些Mean由城市开敞度、绿度、人口密度等属性构成。

data['class']=[0 if data.iloc[i]['lb']!=0  else 1 for i in range(data.shape[0])] # assign the type according to the lb attribute. 
# lb!=0 means inefficient land-use block

# create buffer
data["geometry"]=data.buffer(1200)
data['geometry'].plot()

# construct the adjacency list between two types.
Df=collections.defaultdict(list)
A=data[data['class']==0] # inefficient
B=copy.deepcopy(data) 
isConsidered=collections.defaultdict(list)

for i in range(B.shape[0]):
    # isConsidered：[Is it used?the closest object]
    isConsidered[B.iloc[i].OBJECTID]=[False,B.iloc[i].OBJECTID]

    
# Determining adjacency
import alive_progress
Df=collections.defaultdict(list)
InvertDf=collections.defaultdict(list)
with alive_progress.alive_bar(A.shape[0],force_tty=True) as bar:
    for i in range(A.shape[0]):
        for j in range(B.shape[0]):
            if  (c:=A.iloc[i])["geometry"].intersects((v:=B.iloc[j]).geometry):
                Df[c.OBJECTID].append(v.OBJECTID)
                InvertDf[v.OBJECTID].append(c.OBJECTID)
        bar()

4.2 特征向量构建

# Using Entropy Weight Method
# By the way, perhaps Coefficient of Variation Method is better

def EWM(data):
    t=(data-data.min(axis=0))/(data.max(axis=0)-data.min(axis=0))
    rt=t
    t=t/t.sum(axis=0)
    t[t<0.0001]=0.0001
    entropy=-1/np.log(t.shape[0])*np.sum(t*np.log(t))
    if type(entropy)!=list:
        entropy=[entropy]
    return [(1-i)/(len(entropy)-sum(entropy)) for i in entropy],rt

EWMD=data[['MEAN',"MEAN_1","MEAN_12","MEAN_12_13"]]

features=EWM(EWMD)

# constructing feature vectors
s=s=features[1]*features[0][0]
nf=pd.DataFrame(s)
nf['oid']=data['idd']
nf

4.3 相似度计算

# Using cosine similarity
def sim(x,y):
    return np.dot(x,y)/(np.linalg.norm(x)*np.linalg.norm(y)+1e-6)

如果不考虑空间异质性和时空同现什么的，也许并查集会是一个更高效的数据结构。

class us(object):
    def __init__(self,n):
        self.n=n
        self.init()
    def find(self,x):
        if self.parents[x]!=x:self.parents[x]=self.find(self.parents[x])
        return self.parents[x]
    def union(self,a,b):
        if (ra:=self.find(a))!=(rb:=self.find(b)):
            self.parents[ra]=rb
            self.areas[rb]+=self.areas[ra]
    def init(self):
        self.parents=dict(zip((v:=[self.n.iloc[i].OBJECTID for i in range(self.n.shape[0])]),v))
        self.areas=dict(zip(v,[self.n.iloc[i].Shape_Area for i in range(self.n.shape[0])]))
    def getArea(self,x):
        return self.areas[x]

def InitUS(epsilon=0.9999):
    u=us(A)
    for i in range(A.shape[0]):
        idxI=A.iloc[i].OBJECTID
        for j in range(i+1,A.shape[0]):
            idxJ=A.iloc[j].OBJECTID
            xi,xj=nf[nf['oid']==idxI].values.reshape(-1)[:-1],nf[nf['oid']==idxJ].values.reshape(-1)[:-1]
            if sim(xi,xj)>=epsilon:
                u.union(idxI,idxJ)
    return u
    
u=InitUS() 

很可惜，如果想要考虑空间衰减，并查集需要做一些改进。这里处于时间关系，就直接使用邻接表进行了。

# Iterate over InvertDf to determine the similarity of the source points
def SourceSim(T:list,epsilon=0.9):
    s=A[A['OBJECTID']==T[0]]['Shape_Area'].values[0]
    R=[T[0]]
    if (v:=len(T))==1:
        return s,R

    for n in range(v):
      
        ts=A[A['OBJECTID']==T[n]]['Shape_Area'].values[0]
        record=[T[n]] 
        for m in range(n+1,v):
            xn,xm=nf[nf['oid']==T[n]].values.reshape(-1)[:-1],nf[nf['oid']==T[m]].values.reshape(-1)[:-1]
            if sim(xn,xm)>=epsilon:
                ts+=A[A['OBJECTID']==T[m]]['Shape_Area'].values[0]
            record.append(T[m])
        if s<ts:
            s=ts
            R=record
    return s,R

def bone_area_withInvert(alpha,isConsidered):
    ic=copy.deepcopy(isConsidered)
    newDf=copy.deepcopy(Df)
    with alive_progress.alive_bar(len(newDf),force_tty=True) as bar:
        for i,j in newDf.items():
            x=nf[nf['oid']==i].values.reshape(-1)[:-1]
            for _ in j:
                y=nf[nf['oid']==_].values.reshape(-1)[:-1]
                if ic[_][0]:continue
                S_y=B[B['OBJECTID']==_]['Shape_Area'].values[0]
                S_x,records=SourceSim(InvertDf[_],0.7)
                epsilon=0.5 if S_y<=(0.5*S_x) else 1/(1+np.exp(-alpha*(S_y/S_x)))
                s=sim(x,y)
                for k in records:
                    if (S:=sim(nf[nf['oid']==k].values.reshape(-1)[:-1],y))>=epsilon and S>=s:

                        ic[_][0]=True
                        ic[_][1]=k
            bar()
                        
    return ic


# Assign attributes
data['IsUpdate']=[False if ((v:=data.iloc[i])['class']==0 or IC[v.OBJECTID][0]==False) else True  for i in range(data.shape[0])]
data['Source']=[data.iloc[i].OBJECTID if (v:=data.iloc[i])["class"]==0 else IC[v.OBJECTID][1] for i in range(data.shape[0])]
data['Final']=[0 if (data.iloc[i]['class']==0 or data.iloc[i]['IsUpdate']) else 1 for i in range(data.shape[0])]

# undo buffer
data.geometry=D.geometry

# save
data.to_file(path+r"\\Recluster_with4features_.shp",encoding='utf-8')

五、结果展示

分类结果

部分分类结果

---

六、进一步工作

接下来将抽取10%的重分类斑块由专家进行评估，获得是否需要重建、部分整改、全部整改的标签(A,B,C)

并基于机器学习方法，对剩下的90%的斑块进行计算

基于结构分析模型对其中的隐变量关系进行探究