把 5G 天线波束“画”到地图上：PyVista + Trimesh 波束赋形三维演示

做无线网络规划的同事经常要回答一类问题：这根天线电下倾再压 2°，覆盖会缩多少？把 8×8 阵列换成 4×16，波束会变成什么形状？加了幅度加权，旁瓣到底还剩多少？这些问题在仿真软件里有标准答案，但很少有人能“看见”它——波束是三维的，而我们平时看的覆盖图是二维的。

这次的实验就是把它们放到同一个场景里：铁塔、阵列天线、三维波束、地面信号强度和空间信号体素，全部叠在真实城市（贵阳市中心，OpenStreetMap 建筑）的地图上，调一个参数，后端立刻重新计算并重绘。

• 在线演示：https://study.geobridge.cn/web-mapping/antenna-beam/
• 计算与模型：PyVista + Trimesh（Python）
• Web 服务：FastAPI
• 前端：MapLibre GL + deck.gl

上图是默认场景：一座 30 m 挂高的三扇区站（方位角 0°/120°/240°，每扇区 8×8 阵列，3.5 GHz），半透明的花瓣是每个扇区的三维波束包络，橙色方块是强信号体素，地面是用户终端高度（1.5 m）的 RSRP 热力图，虚线是 −80/−90/−100/−110 dBm 等值线。注意地面那些沿街道放射出去的高亮“条带”——那是射线检测算出来的视距（LOS）街道峡谷，建筑背后的阴影区信号会掉 20～30 dB。

一、整体架构

浏览器（MapLibre 底图 + 建筑拉伸 / deck.gl 三维网格）
   │  POST /api/compute  {站点 + 3 个扇区参数 + 显示选项}
   ▼
FastAPI
   ├─ antenna.py   3GPP 38.901 阵元方向图 + 均匀矩形阵列（URA）波束赋形
   ├─ scene.py     Trimesh：OSM 建筑拉伸、程序化铁塔模型、embree 射线遮挡
   └─ engine.py    PyVista：波束包络曲面、地面等值线、三维体素阈值提取
   │  ← JSON（网格顶点/法线/索引/颜色用 base64 二进制）+ PNG 覆盖图 + GeoJSON 等值线

分工很简单：数学在 NumPy，几何在 Trimesh，网格/体数据处理在 PyVista，渲染在浏览器。后端不画图，只产出“可以直接上 GPU 的网格”，前端也不做任何电磁计算。

一次完整计算（三扇区波束网格 + 9 万点地面覆盖 + 9.7 万体素）在笔记本上大约 0.4～0.6 s；建筑遮挡的射线结果按天线位置缓存，所以拖动下倾、偏转这类参数时只需要重算增益，体感上就是“松手即出图”。

二、天线模型：阵元方向图 × 阵列因子

2.1 单个阵元

阵元采用 3GPP TR 38.901 表 7.3-1 的方向图（θ 是天线坐标系下的天顶角，φ 是方位角）：

• 垂直：A_V(θ) = −min{12·((θ−90°)/65°)², 30 dB}
• 水平：A_H(φ) = −min{12·(φ/65°)², 30 dB}
• 合成：A(θ,φ) = −min{−(A_V + A_H), 30 dB} + 8 dBi

这就是常说的 65° 半功率波宽、30 dB 前后比的面板阵元。

2.2 阵列与波束赋形

把 M 行 × N 列阵元按 d_v、d_h（以波长为单位）排成平面阵。要把波束指向电下倾 θ_s、水平偏转 φ_s，就给每个阵元乘一个相位补偿权值：

w(m,n) = a_m · b_n · exp(−j·2π·(n·d_h·sinθ_s·sinφ_s + m·d_v·cosθ_s))

a、b 是幅度加权（均匀 / Hamming / Taylor −30 dB）。因为权值在行、列两个方向可分离，阵列因子可以拆成两个一维线阵相乘，计算量从 P×M×N 降到 P×(M+N)：

def _axis_factor(psi, psi_s, amp):
    k = np.arange(amp.size)
    phase = np.exp(1j * np.outer(psi - psi_s, k))   # (P, K)
    return np.abs(phase @ amp) ** 2 / np.sum(amp ** 2)

array_gain = af_h * af_v                         # 两个方向的线阵相乘
gain_db = element_gain_db(theta, phi) + 10 * np.log10(array_gain)

单元测试里核对了几个“教科书数字”：8×8 均匀阵列峰值增益 = 8 + 10·lg64 ≈ 26.1 dBi；16 列、半波长间距的水平半功率波宽约 6.3°；加 Hamming 窗后波束变宽——这些都和理论一致，方便后面放心地看图。

2.3 坐标系

世界坐标用以场景原点为中心的局部 ENU（东-北-天，单位米）。每个扇区面板有自己的坐标系：x′ 指向波束主轴，y′ 水平向左，z′ 沿面板向上。方位角（罗盘方向，顺时针）和机械下倾决定这三根轴在 ENU 中的方向，任何一个世界方向向量投影到这三根轴上，就得到天线坐标系下的 (θ, φ)。

上图是直接用 PyVista 离屏渲染的三种阵列方向图（半径 = 相对峰值的增益，色标为 dBi）：均匀加权有明显的一圈旁瓣；Taylor −30 dB 加权几乎把旁瓣抹平，代价是主瓣变胖、峰值增益降约 1.4 dB；4 行 × 16 列的阵列水平方向更窄、垂直方向更宽，偏转 30° 后主瓣明显变“扁”，并且出现了指向另一侧的栅瓣苗头。

三、PyVista：把方向图变成地图上的“波束”

方向图是“相对增益”，直接画在地图上没有尺度。我把它换成一个物理上有意义的量：等信号包络面——沿每个方向，算出 RSRP 恰好等于阈值（默认 −100 dBm）的距离：

RSRP_tx + G(θ,φ) − PL(d) = 阈值   ⇒   d(θ,φ) = 10^((RSRP_tx + G − 阈值 − 13.54 − 20·lg f) / 39.08)

PL 取 38.901 UMa 的城区 NLOS 斜率，这样主瓣伸到几百米、旁瓣短一截、零陷收缩到天线附近——形状和长度都是“真”的。

在 (θ, φ) 网格上（1.5° 步长，121×241 个采样）得到每个方向的半径后，交给 PyVista：

surf = pv.PolyData.from_regular_faces(world.reshape(-1, 3), tris)  # 每个格子两个三角形
surf.point_data["gain"] = g.ravel()
surf = surf.clean(tolerance=1e-6)            # 增益过低、半径为 0 的点焊接到天线中心
surf = surf.clip(normal=(0, 0, 1), origin=(0, 0, 0.3), invert=False)  # 切掉地面以下
surf = surf.triangulate()

两个细节：

1. clean 焊点：把峰值以下 35 dB 的方向半径设为 0，它们全部收缩到同一个点，clean 一次就能合并，网格从 2.9 万点降到 2 千点左右，传输和渲染都轻很多。
2. clip 切地面：下倾的主瓣会“扎”进地下，直接按 z=0.3 m 平面裁掉，得到波束打到地面的切口——这正是地面覆盖最强的那一圈。

一开始我直接用 StructuredGrid.extract_surface().triangulate()，三个扇区要 600 ms，triangulate 占了九成；改成自己用 NumPy 拼三角形索引后降到每扇区 25 ms。

近景里可以看到铁塔模型和三个扇区各自的“花瓣”：亮的是主瓣，暗的是旁瓣（顶点亮度按增益映射，前端再乘上扇区颜色），细线是每个扇区当前的波束指向。

四、Trimesh：建筑、铁塔和遮挡射线

4.1 建筑

263 栋建筑来自 OpenStreetMap（Overpass API 一次性导出）。OSM 里只有 23 栋带 height / building:levels，其余按占地面积和 OSM id 生成确定性的高度（小楼 9～24 m，大楼 24～100 m），保证每次运行场景一致。每栋楼用 trimesh.creation.extrude_polygon 拉伸，最后 concatenate 成一个 6 千面的网格。

4.2 铁塔

铁塔是用 Trimesh 程序化拼出来的：四根收分的主腿、每 6 m 一圈横杆和 X 形斜撑（全部是 6 棱柱 cylinder(segment=...)）、天线平台、避雷针，以及三块朝向各自方位角的天线面板。挂高一变，塔身节数和收分自动跟着变。

4.3 视距 / 非视距：射线检测

传播损耗是否走 NLOS，取决于天线到接收点的连线是否穿过建筑。这一步完全交给 Trimesh：

locs, ray_idx, _ = mesh.ray.intersects_location(origins, dirs, multiple_hits=False)
hit = np.linalg.norm(locs - origin, axis=1)
blocked[ray_idx] = hit < dist[ray_idx] - 0.5

安装了 embreex 时 Trimesh 会自动用 Intel Embree 做加速，9 万条射线对 6 千个三角形也就一百多毫秒。射线结果只和天线位置、网格有关，和下倾、偏转、功率都无关，所以用 lru_cache 按天线坐标缓存——调参时完全不需要再打射线。接收点落在建筑轮廓内时（用 Shapely 的 STRtree 判断）按室内处理，额外加 20 dB 穿透损耗。

4.4 楼顶站

“放置铁塔”时如果点在建筑上，铁塔会自动立在楼顶（基座高度 = 屋顶高度），挂高变成“楼顶以上的高度”。下图把站放到一栋 84 m 的高楼楼顶，再架 12 m 支架——视距街道明显变多，覆盖面积（≥ −100 dBm）从 1.5 km² 涨到 2.3 km²，LOS 比例从 29% 涨到 47%：

五、传播与覆盖：RSRP、等值线

• 路径损耗：3GPP TR 38.901 UMa，LOS：PL = 28 + 22·lg d + 20·lg f；NLOS：取 max(PL_LOS, 13.54 + 39.08·lg d + 20·lg f)。
• RSRP：RSRP 是“每个资源粒子”的接收功率。100 MHz、30 kHz 子载波的 NR 载波有 273×12 = 3276 个子载波，所以 46 dBm（40 W）的扇区总功率折算到每个 RE 是 46 − 10·lg3276 ≈ 10.8 dBm。不做这一步的话，算出来的“信号强度”会比现网读数高 35 dB。
• 最佳服务扇区：每个点取三个扇区中 RSRP 最大的一个。

地面网格默认 8 m 分辨率、±1.2 km，共 90 601 个点。RSRP 用 turbo 色带编码成 PNG（低于 −115 dBm 透明），前端作为 MapLibre 的 image source 贴到地面。等值线用 PyVista 的 marching squares：

img = pv.ImageData(dimensions=(nx, ny, 1), spacing=(res, res, 1), origin=(x0, y0, 0))
img.point_data["rsrp"] = best.ravel()
iso = img.contour([-90], scalars="rsrp")           # 折线段
lines = iso.lines.reshape(-1, 3)[:, 1:]             # 每段两个点 → GeoJSON MultiLineString

六、信号体素：三维的“覆盖图”

地面覆盖只回答“1.5 m 高度有没有信号”。无人机、高层住宅、楼宇室分更关心空间里的信号分布，所以第三层是体素：

grid = pv.ImageData(dimensions=(nx+1, nx+1, nz+1), spacing=(20, 20, 10),
                    origin=(x0, y0, 0))                  # 一个 cell = 一个体素
centers = grid.cell_centers().points                    # 97 200 个体素中心
grid.cell_data["rsrp"] = best_rsrp(centers)             # 同一套天线 + 传播 + 射线模型
kept = grid.threshold(value=-78, scalars="rsrp")        # 只留强信号体素
kept = kept.threshold(0.5, scalars="inside", invert=True)  # 去掉建筑内部

PyVista 的 ImageData 天然就是规则体素网格，threshold 直接返回满足条件的 cell，取 cell_centers 就是要画的体素位置。返回给前端的只有“中心坐标 + RSRP + 主服务扇区”三个数组，前端用 deck.gl 的 SimpleMeshLayer 实例化同一个立方体网格去画，几万个体素依然流畅。

体素有两种看法：全部体素（半透明，看整体空间形状），或者某个高度的水平切片（不透明，看具体数值）。着色也可以从 RSRP 切换为“主服务扇区”，三个扇区的分界一目了然：

七、前端：MapLibre + deck.gl 交错渲染

• 底图与建筑：OpenFreeMap 暗色矢量底图（可切换 Esri 卫星影像）；建筑用 MapLibre 原生 fill-extrusion，和后端射线检测用的是同一份 GeoJSON。
• 三维网格：deck.gl MapboxOverlay({ interleaved: true })，deck 图层和 MapLibre 图层共享同一个 WebGL 上下文和深度缓冲，所以波束会被建筑正确遮挡。
• 坐标：所有网格都用 COORDINATE_SYSTEM.METER_OFFSETS，以场景原点为 coordinateOrigin，后端算出来的 ENU 米坐标原样上屏，不需要在 Python 里做经纬度转换。
• 二进制传输：网格不用 JSON 数组，而是定点数的 base64：顶点坐标 int16（单位分米，±3.2 km 足够）、法线 int8、索引 uint16、每顶点 1 字节亮度；体素中心和 RSRP 也用 int16（0.1 m / 0.1 dB）。前端解码后直接交给 luma.gl 的 Geometry。等值线则先用 vtkStripper（strip(join=True)）把 marching squares 的零碎线段连成折线，再用 Shapely 按 0.4 个格子做 Douglas-Peucker 简化。一次完整结果从最初的 2 MB 降到 460 KB，gzip 后 230 KB。
• 拖动时只要“轻量结果”：滑块拖动过程中的请求不带体素和等值线（只要波束 + 覆盖图），松手后再请求一次完整结果，屏幕上暂时保留上一份体素和等值线。
• 动画：波束网格以天线为原点下发，所以只要每帧改 getScale（0→1）和透明度，就能做出“从天线往外发射”的脉冲；站点底部有一圈扩散的信标；“波束扫描”按钮会让当前扇区的水平偏转在 ±45° 之间往返，每一步都是一次真实的后端计算（扫描时关闭体素以保证流畅）。
• 请求调度：60 ms 防抖，同一时刻只保留一个在途请求，新的参数排队、旧的结果直接覆盖——拖得再快也不会堆积请求。

八、拿它做几个实验

8.1 幅度加权：旁瓣去哪了

只开 α 扇区、侧视。左边均匀加权，主瓣周围有一圈旁瓣（第一旁瓣约 −13 dB）；右边 Taylor −30 dB 加权，旁瓣几乎消失，峰值增益从 26.0 降到 24.6 dBi，水平波宽从 12.6° 变成 15.8°。在城区，旁瓣意味着对邻区的干扰，这就是赋形里常见的“增益换干净”。

均匀加权	Taylor −30 dB

8.2 电下倾：覆盖往回收

同一扇区、机械下倾 0°，电下倾 2° 和 14° 的俯视对比。2° 时主瓣落地点在几百米外，近处靠旁瓣覆盖，等值线呈长条；14° 时能量集中在铁塔附近，≥ −90 dBm 的面积从 0.55 km² 缩到 0.34 km²，换来的是近处更强、更“收得住”的信号——这正是“下倾控制覆盖、减少越区”的直观效果。

电下倾 2°	电下倾 14°

8.3 水平波束偏转：同一块面板，指向不同街道

12 列阵列，水平偏转 −40° 与 +40°。机械方位角没动，波束靠相位权值转到了另一条街上。偏得越多，阵元方向图在该方向的增益越低（峰值从 27.8 降到 23.2 dBi），主瓣也会变宽——这是相控阵扫描角不宜过大的原因。

偏转 −40°	偏转 +40°

九、导出 GLB：同一份网格，换个地方看

面板上的“导出 GLB”会把建筑、铁塔和波束包络打包成一个 glTF 二进制文件（trimesh.Scene.export(file_type="glb")，顺带把 Z 轴朝上转成 glTF 的 Y 轴朝上），可以直接拖进 Blender、three.js editor 或任何 glTF 查看器。下面这张图则是在 Python 里用 PyVista 离屏渲染的同一个场景（建筑 + Trimesh 铁塔 + 波束包络 + −82 dBm 以上的体素），用来核对前端的结果：

十、局限与下一步

这是一个“直觉演示”，不是规划工具，几处简化需要说明：

• 传播模型是 38.901 UMa 的确定性均值，没有阴影衰落、绕射和反射，NLOS 只靠“连线是否穿楼”判断；
• 建筑高度大多是估计值，地形按平地处理；
• 波束包络用的是 NLOS 斜率下的等信号距离，便于在地图上看形状，和真实空间覆盖并不等价（真实的空间覆盖请看体素层）；
• 没有考虑干扰（SINR）与多用户 MU-MIMO 的动态赋形。

下一步比较想做的是：把 SINR 加进去（三扇区互为干扰），以及用真实的 DEM 地形替换平地，让山城贵阳的高差也参与遮挡计算。

---

建筑数据 © OpenStreetMap contributors；底图 © OpenFreeMap / OpenMapTiles，影像 © Esri。