手机信令数据分析（人口分析与城市发展规划）

---

研究区域说明

本次研究区域为北京市，研究范围为东经115°25′至117°30′，北纬39°26′至41°03′之间，包括东城区、西城区、朝阳区、丰台区、石景山区、海淀区等16个区。
本研究将通过对信令数据的深入挖掘和分析，揭示北京市的人口分布、流动规律以及空间特征，基于分析结果，结合城市发展规划的实际需求，对人口分布优化、交通设施布局、公共服务设施配置等提出针对性的规划建议。

---

信令数据格式

信令数据是由移动通信网络记录的用户通信活动相关信息的集合，记录了用户与通信基站之间的交互信息。
为了方便后续的算法处理，将得到的所有信令数据存入数据库，这里选用MySQL数据库。

列名	数据类型	注释
uid	int	用户编号
time	datetime	信令时间
longitude	decimal(10,6)	经度
latitude	decimal(10,6)	纬度

这里只展示一部分信令数据：

用户	时间	经度	纬度
2	20240927 23:47:46	116.254753	39.887227
134	20240927 23:47:47	116.29467	39.85779
1716	20240927 23:47:51	116.281807	39.948059
4873	20240927 23:47:54	116.185365	39.921688
1572	20240927 23:47:55	116.336388	39.790555
767	20240927 23:47:55	116.514214	39.743594
5638	20240927 23:47:56	116.254349	39.901931
489	20240927 23:47:56	116.255012	40.141082
3172	20240927 23:47:56	116.31892	40.00735
3176	20240927 23:47:57	116.278114	40.169339

一共有1w用户，id为从0~9999，信令时间为使用手机进行了一次通信的发生时间，经纬度坐标为基站的坐标位置。

通过对上述这些信令数据的挖掘和分析，我们能够深入了解城市内部的人口分布情况，揭示出不同区域的人口密度和活动热点，为城市资源配置提供科学依据。其次，基于信令数据的研究可以帮助我们更好地理解城市人口的流动模式，分析人口迁移和通勤等行为规律，为交通规划和基础设施建设提供指导。最重要的是，通过深度分析信令数据，我们能够把握城市发展的趋势，为制定科学合理的城市规划和发展策略提供有力支持。

---

数据预处理

数据预处理的过程如下：
（1）缺失数据
缺失数据是指由于通信系统记录有误，导致用户的记录存在某一个或多个关键字段的内容缺失。缺失的内容会影响对用户有效信息的提取，因此需要剔除。
（2）重复数据
重复冗余数据是由于用户在同一基站范围内活动，导致同一用户在同一时间连接同一个基站。这样的数据各个字段的值完全相同，只需要保留一条即可。
（3）乒乓切换
在多个邻近基站交叉覆盖区域，两者信号强度相似，手机会在两者之间频繁进行切换。本文对信令记录进行统计，以基站对的形式考察各对之间的切换数量，短时间内切换数量过多则认为其存在乒乓效应。以乒乓效应发生期间连接次数最多的基站为主基站，剔除期间主基站外的记录。
（4）漂移数据
通信用户突然从邻近基站切换至远处基站一段时间后又切回邻近基站的情况，这种数据可以通过计算用户的移动速度来识别，剔除移动速度超过阈值的信令数据。

伪代码如下：

1.遇到缺失数据，直接剔除，不进行下一步处理
for `row` in `csv_reader`:
	if [uid] and [time] and [longitude] and [latitude]:
     	TODO NEXT

2.遇到完全重复数据，直接跳过
for `row` in `csv_reader`:
	if {uid, time, longitude, latitude} in all_data :
     	continue

3.删除乒乓切换信令数据
# 根据id/time排序 → 统计切换次数 → 阈值识别 → 确定主基站
all_data.sort_by(['uid', 'time'])
for_each(all_data): 
	swap({curr_longitude, curr_latitude}, {next_longitude, next_latitude})
	switch_count(count+1 if swap in all_switch_data)
if switch_count > threshold:
	ping_pong_periods()
if count < max(switch_count):
	delete

4.剔除漂移数据
calculate_distance({curr_longitude, curr_latitude}, {next_longitude, next_latitude})
calculate_time({curr_time}, {next_time})
# 漂移速度过大，或者短时间内两次信令跨越多个基站，删除
If (distance / time > threshold) or isInclude(lotPoints):
	delete

---

人口分布特征分析

这里主要是根据信令数据，划分出人口的工作区域和生活区域。
首先需要做的是，根据经纬度来得到用户所在的北京市区域情况，这里用python解析实现：

# 读取JSON文件，该文件为北京市区域地理位置文件
with open('./files/110000.json', 'r', encoding='utf-8') as file:
    area_local = json.load(file)

# 构建区域边界的shapely Polygon对象
polygons = []
for feature in area_local['features']:
    polygon = shape(feature['geometry'])
    polygons.append((feature['properties']['name'], polygon))

# 解析xinlin.csv信令文件数据, 数据格式[{}, {}...]
data_list = []
with open('./files/xinlin.csv', 'r') as file:
    csv_reader = csv.reader(file)
    for row in csv_reader:
        row_info = {}
        row_info['uid'] = int(row[0])
        row_info['time'] = datetime.strptime(row[1], '%Y-%m-%d %H:%M:%S').strftime('%Y-%m-%d %H:%M:%S')
        row_info['longitude'] = Decimal(row[2])
        row_info['latitude'] = Decimal(row[3])
        row_info['area'] = ''
        point = Point(row_info['longitude'], row_info['latitude'])
        for name, polygon in polygons:
            if polygon.contains(point):
                row_info['area'] = name
                break
        data_list.append(row_info)
    #     print(row_info)
    # print(data_list)

根据上述代码，就可以解析出用户所在北京市的哪个区。现在我们的表结构字段如下：

列名	数据类型	注释
uid	int	用户编号
time	datetime	信令时间
longitude	decimal(10,6)	经度
latitude	decimal(10,6)	纬度
area	varchar(10)	所属区域

然后开始分析用户生活和工作分别在哪些区域。

工作区域人口分布

我们将大多数人的日常工作时间段定为早上9点到下午6点。在这个时间段内，我们统计出每个用户通信次数最多的区域，合理地推断为该用户的工作区域，进而可以统计出每个区域的人口数量。
提取处理的核心算法如下所示：

SELECT `uid`, `area`, count(id) FROM case_all_data
WHERE DATE_FORMAT(time, '%H:%M:%S') BETWEEN {start} AND {end}
GROUP BY `uid`, `area`;

# 根据查询结果，分离出每个用户在9:00-18.00出现次数最多的区域
for info in id_area_frequency:
    if id not in id_area.keys():
        id_area[id] = {
            "area": area,
            "frequency": frequency
        }
    else:
        if frequency > id_area[id]["frequency"]:
            id_area[id] = {
                "area": area,
                "frequency": frequency
            }

通过可视化手段，我们绘制了人口工作区域分布图，如下图：

生活区域人口分布

我们将大多数人回到居住区域或者休息的时间段定为晚上8点到凌晨4点，将每个用户在此时段内最后一次出现的区域确定为该用户的居住区域，统计出每个区域的人口数量。提取算法如下所示：

WITH ranked_data AS (
    SELECT `uid`, `area`, time AS last_occurrence_time,
        ROW_NUMBER() OVER (PARTITION BY `uid` ORDER BY time DESC) AS rn
    FROM case_all_data
    WHERE (time >= '2024-04-08 20:00:00' OR time <= '2024-04-08 04:00:00')
)
SELECT `uid`, `area`, last_occurrence_time
FROM ranked_data
WHERE rn = 1;

通过可视化手段，我们绘制了人口居住区域分布图，如下图：

---

人口流动规律分析

通过提取用户的移动轨迹，记录OD数据，我们可以清晰地了解到城市内不同区域之间的人口流动情况，包括流动的数量、方向以及时间分布等。

OD数据（Origin-Destination Data）是描述出行流量的一种形式，其核心作用在于记录不同起点（Origin）与终点（Destination）之间的出行关系

为揭示北京人口在不同时间段内的流动规律，我们首先对处理好的北京市信令数据按照时间顺序，将一天内的24小时划分为12个时间段，每个时间段为两个小时。然后，我们提取了每个时间段内每个用户的移动轨迹数据，并根据用户的起始点和目的地提取了OD数据，统计在不同时间段内每一对OD的移动路径和数量。
这里拿其中一个时间段（08:00–10:00）举例，代码如下：

sql_select = """
    SELECT `uid`, `area` FROM case_all_data WHERE `time` BETWEEN %s AND %s
"""
try:
    cursor.execute(sql_select, ("2024-09-27 08:00:00", "2024-09-27 10:00:00"))
except Exception as err:
    print(f"数据库查询失败: ", err)
id_area = cursor.fetchall()

# 1.计算每个用户在08:00--10:00这一时间段的区域移动轨迹
user_area = {}
for info in id_area:
    if info[0] not in user_area.keys():
        user_area[info[0]] = [info[1]]
    else:
        if user_area[info[0]][-1] != info[1]:
            user_area[info[0]].append(info[1])

# 2.计算人口跨区域流动次数            
beginEnd_count = {}
for uid, area_list in user_area.items():
    for i in range(len(area_list)-1):
        key = area_list[i] + '-' + area_list[i+1]
        key_reverse = area_list[i+1] + '-' + area_list[i]
        if key not in beginEnd_count.keys():
            if key_reverse not in beginEnd_count.keys():
                beginEnd_count[key] = 1
            else:
                beginEnd_count[key_reverse] += 1
        else:
            beginEnd_count[key] += 1

# isAll选择是否把没有流动的区域也加入置为0
if isAll:
    all_area = ["怀柔区", "海淀区", "朝阳区", "丰台区", "房山区", "大兴区", 
                "昌平区", "西城区", "东城区", "通州区", "石景山区","顺义区",
                "平谷区", "密云区", "门头沟区", "延庆区"]
    for begin in all_area:
        for end in all_area:
            if begin != end:
                begin_end = begin + '-' + end
                end_begin = end + '-' + begin
                if begin_end not in beginEnd_count.keys() and end_begin not in beginEnd_count.keys():
                    beginEnd_count[begin_end] = 0

# 写入 CSV 文件
file_start_time = datetime.strptime(start_time, '%Y-%m-%d %H:%M:%S').strftime('%H')
file_end_time = datetime.strptime(end_time, '%Y-%m-%d %H:%M:%S').strftime('%H')
fileName = f'./output/move_data_{file_start_time}_{file_end_time}.csv'
with open(fileName, mode='w', newline='', encoding='utf-8') as file:
    writer = csv.writer(file)
    writer.writerow(['origin', 'destination', 'count'])
    for key, value in beginEnd_count.items():
        begin = key.split('-')[0]
        end = key.split('-')[1]
        row_list = [begin, end, value]
        if begin and end:
            writer.writerow(row_list)

cursor.close()
db.close()

为了更直观地展示这些数据，我们绘制了全天12个时段的OD图。这些图显示了不同区域间的人口流动情况，连线的颜色和粗细代表此条OD的人数。因为我们绘制的OD图中是没有区分方向的，所以，我们为每个时间段绘制了一个并排柱状图，其中包括了各个区域对应时间段内的流出人口、流入人口和内部流动人口的数量。
如下是时间段（08:00–10:00）的人口流动图：


接下来我们从另一个角度来分析北京市居民在不同时间段内的流动规律。我们按照上述工作区/生活区的划分方法，将9:00至18:00时间段内用户出现次数最多的区域识别为工作区域，将20:00至4:00时段内最后一次出现的区域识别为居住区域。我们设定条件：早晨（7:00-10:00）的出行中，起点在生活区且终点在工作区的OD对被认为是通勤OD；傍晚（17:00-20:00）的出行中，起点在工作区且终点在生活区的OD对被认为是通勤OD。我们根据这些条件从出行数据中筛选出符合条件的OD对，并匹配它们的起点和终点，形成完整的通勤路径。核心算法如下所示：

1.统计人口工作区-生活区(生活区-工作区)的跨区域流动情况
sql_select = """
    SELECT 
        live.area, work.area 
    FROM 
        case_live_area as live, case_work_area as work 
    where 
        live.uid = work.uid
"""
cursor.execute(sql_select)
live_work = cursor.fetchall()

2.计算人员工作区-生活区(生活区-工作区)的跨区域流动数量
for info in live_work:
    if info[0] != info[1]:
        key = info[0] + '-' + info[1]
        if key not in beginEnd_count.keys():
            beginEnd_count[key] = 1
        else:
            beginEnd_count[key] += 1

我们对所有的通勤OD对进行统计，然后使用地理信息系统软件将通勤OD的路径可视化展示在地图上，通过有向箭头连接两个时间段通勤的起点和终点，箭头的颜色和粗细表示通勤人数的多少，如下图所示：

---

人口空间特征分析

通过信令数据提取人口的空间特征，如聚居区、商圈、交通枢纽等。分析这些特征对城市规划的影响，并提出相应的优化建议。以下主要是针对根据信令数据分析提取出交通枢纽。
对于交通枢纽，本文设计算法对交通枢纽进行识别。
步骤1：计算每个用户下一个位置的基站经纬度和时间，即根据用户id和时间进行分组和排序，得到每个用户按时间顺序的轨迹。
步骤2：根据相差经纬度计算地理距离，相差时间算出时间间隔，求出每段的相对用户速度，根据阈值去掉漂移数据。这里选取100m/s作为速度阈值，因为高铁的最高时速为100m/s，所以认为用户最大的移动速度必须小于阈值。
步骤3：对所有用户的起点经纬度、终点经纬度分组，求每个起点-终点组合的数量，只保留高频移动路径。
步骤4：统计所有高频移动路径的起点次数、终点次数，同时统计每个点总的经过次数。设定高于50次为阈值，大于阈值的点则被判定为交通枢纽。

又新增了字段，现在我们的数据格式如下：

列名	数据类型	注释
uid	int	用户编号
time	datetime	信令时间
longitude	decimal(10,6)	经度
latitude	decimal(10,6)	纬度
next_longtitude	decimal(10,6)	下一时刻经度
next_latitude	decimal(10,6)	下一时刻纬度
next_time	datetime	下一信令时刻
distance	decimal(10,4)	距离
speed	decimal(10,2)	速度
area	varchar(10)	所属区域

本文所设计的核心算法如下所示：

1.计算用户下一个时刻的位置和时间
INSERT INTO temp_case_unique_data 
    (uid, time, next_longitude, next_latitude, next_time)
SELECT
    uid,
    time,
    LEAD(longitude) OVER (PARTITION BY uid ORDER BY time) AS next_longitude,
    LEAD(latitude) OVER (PARTITION BY uid ORDER BY time) AS next_latitude,
    LEAD(time) OVER (PARTITION BY uid ORDER BY time) AS next_time
FROM case_unique_data;
UPDATE case_unique_data cud1
    JOIN temp_case_unique_data cud2
    ON cud1.uid = cud2.uid AND cud1.time = cud2.time
SET cud1.next_longitude = cud2.next_longitude,
    cud1.next_latitude = cud2.next_latitude,
    cud1.next_time = cud2.next_time;

2.计算用户移动的相对地理距离、时间差以及速度
DELIMITER $$
CREATE FUNCTION calculate_distance(lat1 DECIMAL(10, 6), lon1 DECIMAL(10, 6), lat2 DECIMAL(10, 6), lon2 DECIMAL(10, 6))
RETURNS DECIMAL(10, 4)
DETERMINISTIC
BEGIN
    DECLARE R INT DEFAULT 6371000;
    DECLARE dLat DECIMAL(20, 10);
    DECLARE dLon DECIMAL(20, 10);
    DECLARE a DECIMAL(20, 10);
    DECLARE c DECIMAL(20, 10);
    SET dLat = RADIANS(lat2 - lat1);
    SET dLon = RADIANS(lon2 - lon1);
    SET a = SIN(dLat / 2) * SIN(dLat / 2) + COS(RADIANS(lat1)) * COS(RADIANS(lat2)) * SIN(dLon / 2) * SIN(dLon / 2);
    SET c = 2 * ATAN2(SQRT(a), SQRT(1 - a));
    RETURN R * c;
END $$
DELIMITER ;

3.统计高频移动路径的起点/终点次数
CREATE TABLE case_start_points/case_end_points AS
SELECT longitude, latitude, COUNT(*) AS start_count
FROM case_path_counts
GROUP BY longitude/next_longitude, latitude/next_latitude;

4.合并起点和终点统计数据
CREATE TABLE case_merged_points AS
SELECT csp.longitude, csp.latitude, csp.start_count, IFNULL(cep.end_count, 0),
(csp.start_count + IFNULL(cep.end_count, 0)) AS total_count
FROM case_start_points csp
LEFT JOIN case_end_points cep
ON csp.longitude = cep.longitude AND csp.latitude = cep.latitude
UNION
SELECT cep.longitude, cep.latitude, IFNULL(csp.start_count, 0), cep.end_count,
(IFNULL(csp.start_count, 0) + cep.end_count) AS total_count
FROM case_start_points csp
RIGHT JOIN case_end_points cep
ON csp.longitude = cep.longitude AND csp.latitude = cep.latitude
WHERE csp.longitude IS NULL;

根据以上步骤，绘制出了北京市交通枢纽图，如下图所示：


如上图所示，图中圈代表着识别出来的交通枢纽点，红色圆圈的大小与人口流动的数量成正比。圆圈越大，人口流动数量越多；圆圈越小，人口流动数量越少。

为了验证本文所提取交通枢纽的有效性，我们使用以下代码在高德地图开放平台爬取相应坐标点半径范围内的交通枢纽点。由于数据集中所展示为基站经纬度，与用户实际经纬度存在偏差，所以选择指定半径范围内的交通枢纽进行匹配。通过这种方式，可以对比分析提取出的交通枢纽点与实际交通枢纽的匹配度，以确保提取结果的准确性和可靠性。代码如下：

url = f"https://restapi.amap.com/v3/place/around?key={key}&location={longitude},{latitude}&radius={radius}&keywords={keywords}&output=json"
response = requests.get(url)
return data = response.json()

通过验证发现，当选择半径为500米时，所有提取出的点在该半径范围内都存在实际的交通枢纽，准确率达到100%。为了进一步提高准确度，我们将半径缩小至200米，此时91%的提取点在该半径范围内存在实际交通枢纽，这表明我们的交通枢纽识别算法非常有效。下表展示了所识别的部分对应点的经纬度，以及200米范围内实际交通枢纽的分布情况。由于数据量较大，仅展示其中10条记录。

经度	纬度	实际交通枢纽
116.303258	39.986733	海淀桥北(公交站)，万泉河桥(地铁站)
116.29826	39.98953	海淀公园(公交站), 万泉河桥(地铁站)
116.192888	39.912758	古城公园(公交站), 古城东街(公交站)
116.862688	40.372084	保利站(公交站)
116.472627	39.908643	万达广场(公交站), 大望路(地铁站), 核心区(地铁站)
116.65564	40.13787	双兴小区(公交站), 胜利小区(公交站)
116.333921	39.731071	大兴区医院(公交站), 大兴区医院(黄村西大街)(公交站), 黄村西大街(地铁站)
116.686567	40.106039	港馨小学(公交站)
116.638632	39.905164	通州北苑路口东(公交站), 北苑路口站(公交站), 通州北苑(地铁站)
116.265611	40.15105	北京师范大学(公交站), 北沙河西三路南站(公交站), 满井路东口(公交站)