1. 项目概述当TSP遇上中国城市地理数据最近在做一个关于路径优化的项目核心是经典的旅行商问题。为了测试算法的实际效果我需要一个贴近真实场景的数据集。直接使用网上现成的TSPLIB标准库比如berlin52、eil51总觉得差点意思那些大多是欧洲或美国的城市坐标缺乏一点“本土感”。于是我决定自己动手构建一个“中国34个城市”的TSP数据集。这个数据集不仅要包含城市的经纬度WGS84坐标系还要将其转换为平面坐标如高斯-克吕格投影以便进行不同维度和算法下的对比实验。这听起来像是个简单的数据收集工作但实际操作中从坐标获取、格式转换到数据清洗每一步都有不少门道和容易踩的坑。今天就把这个完整的数据集构建过程、核心原理以及我趟过的那些“雷”分享出来希望能给同样需要处理地理空间数据做算法验证的朋友们一些参考。这个数据集特别适合以下几类朋友一是正在学习或研究组合优化、运筹学、智能算法的同学需要一个有现实地理意义的TSP算例二是从事物流路径规划、车辆调度等领域的工程师想用真实中国城市数据测试模型三是任何对地理信息系统数据处理感兴趣想了解经纬度与平面坐标转换原理的开发者。通过这个项目你不仅能获得一个即拿即用的数据集更能掌握一套从原始地理信息到规整算法输入数据的标准处理流程。2. 数据核心城市选取与原始经纬度获取构建数据集的第一步是确定城市名单并获取其准确的经纬度坐标。我选择了中国34个省级行政区的省会/首府城市及部分重要城市这既保证了地理分布的广泛性覆盖全国又使问题规模34个点适中既不会因点数太少而失去优化意义也不会因点数太多而让一些基础算法难以演示。2.1 城市名单确定与数据源选择我最终确定的34个城市名单如下北京、上海、广州、深圳、成都、杭州、武汉、重庆、南京、西安、长沙、沈阳、青岛、郑州、济南、哈尔滨、长春、石家庄、天津、太原、合肥、南昌、福州、昆明、南宁、海口、贵阳、兰州、西宁、银川、乌鲁木齐、呼和浩特、拉萨、香港、澳门。选择这些城市主要基于其行政地位、经济重要性和地理代表性。获取经纬度坐标最直接的方法是使用公开的地理编码API如高德地图或百度地图的开放平台。这些API通常提供免费的额度足以满足个人项目需求。以高德地图地理编码API为例你需要注册开发者账号获取一个Key然后通过构造HTTP请求将城市名作为地址参数发送API会返回包含经纬度的JSON数据。这里有一个关键点务必确认API返回的坐标系。高德地图使用的是GCJ-02坐标系俗称“火星坐标”而百度地图使用的是BD-09坐标系。这两种都是在WGS-84GPS使用的全球标准坐标系基础上进行了加密偏移的。对于学术研究或算法测试为了数据的通用性和可对比性我强烈建议使用WGS-84坐标。注意直接使用高德/百度API返回的坐标进行国际论文中的算法对比可能会因为坐标系不同而导致距离计算出现偏差通常几百米到几公里。如果条件允许应寻找直接提供WGS-84坐标的数据源或使用可靠的坐标转换库进行逆转换。2.2 实操批量获取与数据清洗手动一个个查34个城市的坐标显然不现实。我写了一个简单的Python脚本利用requests库批量调用高德地理编码API。脚本的核心循环是遍历城市列表为每个城市名拼接请求URL。这里有个细节为了增加成功率可以在城市名后加上“市”字如“北京市”、“上海市”。import requests import pandas as pd import time # 你的高德API Key api_key your_amap_api_key_here base_url https://restapi.amap.com/v3/geocode/geo cities [北京, 上海, 广州, ...] # 34个城市列表 data [] for city in cities: params { key: api_key, address: f{city}市 } resp requests.get(base_url, paramsparams) result resp.json() if result[status] 1 and len(result[geocodes]) 0: location result[geocodes][0][location] # 格式经度,纬度 lon, lat location.split(,) data.append({city: city, lon: float(lon), lat: float(lat)}) else: print(fFailed to geocode {city}: {result.get(info, Unknown error)}) time.sleep(0.2) # 礼貌性延时避免请求过快 df pd.DataFrame(data) df.to_csv(china_cities_gcj02.csv, indexFalse)运行脚本后你会得到一个包含城市名、经度、纬度的CSV文件。但如前所述这是GCJ-02坐标。我个人的做法是使用一个经过广泛验证的Python库如coord-convert或gcoord将这些坐标批量转换为WGS-84坐标。转换存在一定精度损失但对于城市级尺度的TSP问题这个误差在可接受范围内。最终我得到了一个名为china_cities_wgs84.csv的纯净WGS-84坐标数据集。3. 坐标转换从球面经纬度到平面投影TSP问题的经典定义是在欧几里得平面上计算点与点之间的直线距离。然而地球是一个球体更准确地说是个椭球体其表面的经纬度是球面坐标。直接使用经纬度计算“直线距离”得到的是球面大圆距离如Haversine公式计算的距离这与平面上的欧氏距离有本质不同。为了在真正的欧氏平面上进行TSP求解和可视化我们需要将经纬度投影到二维平面上。3.1 为什么必须进行投影转换假设我们有两个城市A和B。用它们的经纬度(lon1, lat1)和(lon2, lat2)通过Haversine公式可以计算出它们在地球表面的实际最短路径长度大圆距离。但如果我们错误地将(lon, lat)直接当作平面直角坐标系的(x, y)来计算欧氏距离sqrt((lon2-lon1)^2 (lat2-lat1)^2)这个结果在数值和量纲上都是毫无意义的因为一度经度在不同纬度上的实际长度是不同的。因此要进行基于欧氏距离的TSP求解投影转换是必须的。投影会将地球表面的一块区域比如整个中国映射到一个平面上并尽可能保持距离、面积或形状中的某些特性。对于路径优化问题我们更关心距离的相对准确性。3.2 投影选择高斯-克吕格投影与Web墨卡托对于中国范围内的城市高斯-克吕格投影是一个行业标准选择。它是一种横轴墨卡托投影具有等角特性在小范围内能很好地保持形状和距离比例。中国常用的国家大地坐标系如CGCS2000就基于高斯-克吕格投影并分带如3度带、6度带来减小变形。将WGS-84经纬度转换到CGCS2000 3度带平面坐标能获得很高的精度。然而高斯-克吕格投影的计算相对复杂需要知道带号、中央子午线等参数。对于快速原型验证或希望与在线地图如Leaflet、OpenStreetMap对齐的可视化Web墨卡托投影是一个更方便的选择。它是谷歌地图、必应地图等网络地图服务使用的标准投影其EPSG代码为3857。它将地球模拟为一个球体并进行投影计算简单但在高纬度地区面积变形较大。不过对于主要位于中纬度的中国城市群其距离相对关系仍然具有较好的参考价值。在这个项目中我决定同时提供两种坐标WGS-84经纬度原始球面坐标用于计算真实地理距离。Web墨卡托平面坐标方便进行基于欧氏距离的算法测试和Web可视化。3.3 实操使用PyProj进行批量投影转换Python的pyproj库是处理地理坐标转换的利器。以下是如何将WGS-84经纬度批量转换为Web墨卡托坐标的示例代码import pandas as pd from pyproj import Transformer # 读取之前保存的WGS-84数据 df pd.read_csv(china_cities_wgs84.csv) # 定义坐标转换器从WGS84 (EPSG:4326) 到 Web墨卡托 (EPSG:3857) transformer Transformer.from_crs(EPSG:4326, EPSG:3857, always_xyTrue) # 应用转换 # 注意transformer.transform 期望顺序是 (经度, 纬度) df[x_web], df[y_web] zip(*df.apply(lambda row: transformer.transform(row[lon], row[lat]), axis1)) # 保存结果 df.to_csv(china_cities_coordinates_full.csv, indexFalse) print(df.head())转换后x_web和y_web就是以米为单位的Web墨卡托平面坐标。它们的数值通常会很大例如北京的坐标大约是(1.3e7, 4.8e6)。在用于TSP计算时为了数值稳定和可视化方便我们通常会对这些坐标进行归一化或缩放平移使其分布在一个合理的范围内比如0到1000之间但这不会改变点与点之间的相对距离关系因此不影响优化结果。实操心得pyproj库的版本需要注意。在旧版本中你可能需要使用Proj类而新版本推荐使用Transformer类。always_xyTrue参数非常重要它明确了输入输出顺序是(经度, 纬度)/(x, y)避免了常见的经纬度顺序混淆问题。4. 数据集的构建与格式化有了转换好的坐标下一步就是将它们格式化成TSP标准格式和算法友好的格式方便直接使用。4.1 TSPLIB格式输出TSPLIB是TSP问题的一个通用数据格式标准。一个典型的TSPLIB文件.tsp包含问题名称、类型、维度、节点坐标等信息。虽然我们的坐标是投影后的平面坐标但TSPLIB格式本身支持EUC_2D二维欧氏距离类型。我们需要将归一化后的平面坐标写入文件。首先对Web墨卡托坐标进行归一化处理# 继续之前的df # 归一化到 [0, 1000] 区间保持纵横比 x_min, x_max df[x_web].min(), df[x_web].max() y_min, y_max df[y_web].min(), df[y_web].max() range_x x_max - x_min range_y y_max - y_min scale 1000 / max(range_x, range_y) # 以最大范围缩放到1000 df[x_norm] (df[x_web] - x_min) * scale df[y_norm] (df[y_web] - y_min) * scale然后生成TSPLIB文件tsp_content fNAME: China_34_Cities TYPE: TSP COMMENT: 34 major cities in China (Web Mercator Projection, Normalized) DIMENSION: {len(df)} EDGE_WEIGHT_TYPE: EUC_2D NODE_COORD_SECTION for idx, row in df.iterrows(): tsp_content f{idx1} {row[x_norm]:.2f} {row[y_norm]:.2f}\n tsp_content EOF with open(china34.tsp, w) as f: f.write(tsp_content)这样生成的china34.tsp文件可以被大多数TSP求解器如Concorde, LKH或研究代码直接读取。4.2 算法友好格式JSON与Python字典对于在Python环境中快速进行算法实验JSON或Python原生数据结构更方便。我通常会生成一个包含所有信息的JSON文件import json data_dict { cities: [], coordinates: { wgs84: [], # [[lon, lat], ...] web_mercator: [], # [[x, y], ...] normalized: [] # [[x_norm, y_norm], ...] } } for _, row in df.iterrows(): data_dict[cities].append(row[city]) data_dict[coordinates][wgs84].append([row[lon], row[lat]]) data_dict[coordinates][web_mercator].append([row[x_web], row[y_web]]) data_dict[coordinates][normalized].append([row[x_norm], row[y_norm]]) with open(china_34_cities_data.json, w, encodingutf-8) as f: json.dump(data_dict, f, ensure_asciiFalse, indent2)这个JSON文件结构清晰包含了城市名、三种坐标非常适合在Jupyter Notebook或脚本中加载和使用。5. 距离矩阵计算球面距离 vs 平面距离TSP问题的核心是距离矩阵。基于我们已有的数据可以计算两种距离矩阵用于不同目的。5.1 真实地理距离矩阵球面使用WGS-84经纬度通过Haversine公式计算两两城市之间的球面大圆距离单位公里。这个距离矩阵反映了城市间真实的旅行距离假设地球是完美球体。import numpy as np from math import radians, sin, cos, sqrt, atan2 def haversine_distance(lon1, lat1, lon2, lat2): 计算WGS84坐标下两点间的大圆距离公里 R 6371.0 # 地球平均半径公里 lon1, lat1, lon2, lat2 map(radians, [lon1, lat1, lon2, lat2]) dlon lon2 - lon1 dlat lat2 - lat1 a sin(dlat/2)**2 cos(lat1) * cos(lat2) * sin(dlon/2)**2 c 2 * atan2(sqrt(a), sqrt(1-a)) return R * c # 假设coords_wgs84是N*2的数组每一行是[经度, 纬度] coords_wgs84 df[[lon, lat]].values n len(coords_wgs84) dist_geo np.zeros((n, n)) for i in range(n): for j in range(i1, n): d haversine_distance(coords_wgs84[i,0], coords_wgs84[i,1], coords_wgs84[j,0], coords_wgs84[j,1]) dist_geo[i, j] d dist_geo[j, i] d np.savetxt(distance_matrix_geo_km.csv, dist_geo, delimiter,, fmt%.3f)5.2 欧氏距离矩阵平面使用归一化后的平面坐标(x_norm, y_norm)直接计算欧氏距离。这个距离矩阵用于大多数基于欧氏距离假设的TSP算法如经典的2-opt、Lin-Kernighan启发式算法。coords_norm df[[x_norm, y_norm]].values # 利用NumPy广播高效计算欧氏距离矩阵 diff coords_norm[:, np.newaxis, :] - coords_norm[np.newaxis, :, :] dist_euclidean np.sqrt(np.sum(diff**2, axis-1)) np.savetxt(distance_matrix_euclidean_norm.csv, dist_euclidean, delimiter,, fmt.3f)关键对比比较dist_geo和dist_euclidean矩阵你会发现它们的数值比例和城市间相对远近关系大体一致但并非线性对应。例如乌鲁木齐到哈尔滨的真实地理距离非常远在平面投影上也被拉开而长三角城市群内部的距离两种矩阵都显示得很近。这验证了投影转换在保持局部相对距离关系上的有效性。6. 算法测试以自组织映射为例有了标准格式的数据集和距离矩阵就可以愉快地进行各种TSP算法测试了。这里以“自组织映射”这个热搜词相关的算法——神经气体网络的一个变种常用于路径规划——为例演示如何使用我们的数据集。自组织映射用于TSP的基本思想是在城市分布的平面上初始化一个环状的神经元链通过竞争学习让这个链逐渐“吸附”到各个城市点上最终神经元的连接顺序就构成了一个哈密顿回路TSP解。我们使用归一化的平面坐标进行演示因为SOM通常工作在欧氏空间。这里使用minisom库的一个简单实现import numpy as np from minisom import MiniSom import matplotlib.pyplot as plt # 准备数据归一化后的城市坐标 cities df[city].tolist() X df[[x_norm, y_norm]].values # 初始化SOM1xN的拓扑N等于城市数形成一个环 som_shape (1, len(X)) som MiniSom(som_shape[0], som_shape[1], input_len2, sigma1.0, learning_rate0.5, neighborhood_functionbubble, random_seed10) som.random_weights_init(X) # 训练 num_iter 10000 som.train_random(X, num_iter, verboseFalse) # 获取获胜神经元每个城市对应的最近神经元索引 winners np.array([som.winner(x) for x in X]).squeeze() # 形状 (34, 2) # 由于我们的SOM是1维环我们只关心第二维的索引在环上的位置 winning_indices winners[:, 1] # 根据神经元索引排序得到城市访问顺序 tour_order np.argsort(winning_indices) tour [cities[i] for i in tour_order] print(SOM初步求解的路径顺序起点和终点未连接) print( - .join(tour)) # 可视化 plt.figure(figsize(12, 8)) plt.scatter(X[:, 0], X[:, 1], cred, s100, labelCities, zorder5) for i, city in enumerate(cities): plt.annotate(city, (X[i, 0], X[i, 1]), fontsize9, haright) # 按顺序连接城市 tour_coords X[tour_order] plt.plot(tour_coords[:, 0], tour_coords[:, 1], b-, linewidth1, alpha0.6, labelSOM Path) plt.plot([tour_coords[-1, 0], tour_coords[0, 0]], [tour_coords[-1, 1], tour_coords[0, 1]], b-, linewidth1, alpha0.6) # 闭合回路 plt.legend() plt.title(TSP Solution using Self-Organizing Map (34 China Cities)) plt.xlabel(Normalized X (Web Mercator)) plt.ylabel(Normalized Y (Web Mercator)) plt.grid(True, alpha0.3) plt.show()这段代码会输出一个初步的路径顺序并可视化。需要指出的是基础的SOM用于TSP效果可能不如专门的启发式算法如LKH且可能产生交叉路径。但它提供了一个非常直观的“神经网络学习路径”的演示。得到的路径顺序可以作为更高级优化算法如2-opt局部搜索的初始解从而快速得到一个质量不错的解。7. 常见问题与实战避坑指南在整个数据集构建和算法测试过程中我遇到了不少典型问题。这里总结一下希望能帮你省去不少折腾的时间。7.1 坐标获取与转换中的坑问题1API返回的坐标是“火星坐标”怎么办正如前面提到的国内主流地图API返回的是GCJ-02坐标。如果你需要严格的WGS-84坐标有以下几个选择使用专业转换库如Python的coord-convert库它提供了gcj02_to_wgs84函数。务必使用star数多、维护活跃的库。寻找替代数据源一些开源地理数据集如Natural Earth或学术机构发布的数据可能直接提供WGS-84坐标。OSMOpenStreetMap的Nominatim API也可以查询但其坐标精度和覆盖率需要验证。明确标注如果只是用于算法对比测试且所有对比算法都使用同一套GCJ-02坐标那么相对距离关系仍然是成立的可以在论文中明确说明坐标系即可。问题2投影转换后坐标值巨大导致算法数值不稳定这是使用Web墨卡托坐标的典型现象。务必进行归一化或标准化。除了之前提到的线性缩放归一化到[0, 1000]也可以使用sklearn.preprocessing.StandardScaler进行标准化均值为0方差为1这能改善一些基于梯度或距离的算法的收敛性。问题3转换后的平面坐标东西方向的距离看起来被拉长了这是Web墨卡托投影在高纬度地区的特性。它保持了形状等角但牺牲了面积和距离的均匀性。越往北东西向拉伸越明显。如果你追求更高精度的平面距离应该使用高斯-克吕格投影分带。可以使用pyproj进行更精确的转换例如转换到CGCS2000 / 3-degree Gauss-Kruger zone (如EPSG:4547对应Zone 39)。但这需要你知道每个城市所在的投影带处理起来更复杂。7.2 数据格式与算法适配问题问题4生成的TSPLIB文件无法被求解器识别检查以下几点文件后缀确保是.tsp。格式严格性TSPLIB对空格、换行、章节名称非常敏感。确保NODE_COORD_SECTION之后每行是“序号 X Y”序号从1开始。最好用经典算例如eil51.tsp的格式对照。EDGE_WEIGHT_TYPE如果你提供的是归一化后的平面坐标使用EUC_2D。如果你提供的是经纬度并希望求解器计算球面距离应使用GEO或ATT特殊的伪欧氏距离但这并非所有求解器都支持。最稳妥的方式还是提供平面坐标用EUC_2D。问题5自组织映射SOM或其他算法得到的路径有交叉或明显不合理这很正常SOM本身不保证生成无交叉路径。务必后处理2-opt局部搜索这是解决TSP最经典、最有效的局部优化方法之一。它可以快速消除路径交叉大幅缩短总距离。几乎在任何启发式算法得到初始解后都应该跑一遍2-opt。多次随机初始化SOM对初始权重敏感。可以多次运行取最好的一次结果作为初始解再进行2-opt优化。调整SOM参数sigma邻域半径和learning_rate需要调整。通常训练初期用较大的sigma和learning rate后期逐渐衰减模拟退火思想。7.3 可视化与结果分析问题6如何在地图上可视化TSP路径如果你有城市的经纬度强烈推荐使用folium库进行交互式地图可视化。import folium # 以第一个城市为地图中心 center_lat df.iloc[0][lat] center_lon df.iloc[0][lon] m folium.Map(location[center_lat, center_lon], zoom_start4) # 添加城市标记 for _, row in df.iterrows(): folium.CircleMarker( location[row[lat], row[lon]], radius5, popuprow[city], colorred, fillTrue ).add_to(m) # 假设 best_tour_indices 是最优路径的城市索引列表 best_tour_coords df.iloc[best_tour_indices][[lat, lon]].values # 将路径连接起来注意首尾相连形成回路 folium.PolyLine(locationsbest_tour_coords.tolist() [best_tour_coords[0].tolist()], colorblue, weight2.5, opacity0.8).add_to(m) m.save(tsp_solution_map.html)生成的HTML文件可以在浏览器中打开你会看到一个可缩放、可拖动的中国地图上面用蓝线清晰地标出了TSP的优化路径效果非常直观。问题7如何评估解的质量对于34个城市我们可以用著名的Concorde精确求解器如果可用求出最优解作为基准。或者使用Christofides算法等近似算法得到一个理论上的上界。更实际的做法是用多个启发式算法如LKH, 蚁群算法遗传算法在同一数据集上运行多次比较它们的平均解和最佳解。记录总距离根据你使用的距离矩阵和计算时间。对于我们自己构建的数据集由于没有已知的最优解不同算法之间的相对比较和路径的视觉合理性就成为重要的评估依据。构建这个数据集的过程让我对地理空间数据处理的细节有了更深的认识。从坐标系的纷繁复杂到投影变换的数学原理再到最终为算法提供干净、可靠的输入每一步都需要仔细考量。最终生成的数据集和全套代码我已经整理好放在了GitHub上。当你需要测试一个路径规划算法又希望数据有点“中国味”的时候希望这个“中国34城”数据集能帮上忙。数据处理中最大的体会就是明确需求统一标准。想清楚你的算法到底需要哪种距离真实地理距离还是平面欧氏距离然后从数据源头到最终格式始终保持坐标系和计算标准的一致这样才能得到可靠、可复现的结果。