第03章 - 空地一体时空融合

无人机在地图上划一条航线只要三秒钟，摄像头拍一段视频也只要按一下录制，可要把这两样东西真对上，你会发现它们压根不在同一个世界里：视频是像素，航线是经纬度，飞控说的是上电后的毫秒计数，摄像机说的是 PTS。本章就解决这件事——把空中平台的数据与地面 GIS 底座放进同一套坐标、同一条时间轴、同一个精度口径，然后才谈得上配准、融合、联动。

读完这一章，你能自己跑通一条最小的空地一体链路：把像素坐标换算到经纬度、把 POS 与视频帧对齐到 20 ms 以内、用 2.5D 代价地图规划一条带安全裕度的巡检航线、用平面标记引导 3 m 内自主降落、把实时视频按真视角贴到三维地球上、最后用 ffmpeg + COLMAP（或 VGGT 前馈）+ OpenDroneMap 出正射、用同一套位姿训 3DGS。三维底座按 3D Tiles 1.1 与带 LOD 的 splat 分层，精度按底图来源分级承诺。所有环节都给出可复制的命令与验收数字，而不是概念图。

学习目标：掌握七个坐标系与一条转换链，能用 pyproj + numpy 独立实现并做往返校验；能算出 GSD、重叠率、航高与航时；能写出 solvePnP 位姿解算与视觉伺服状态机；能在「地面四边形贴图」和「按 DEM 逐像素投影」之间做取舍、能按 3D Tiles 1.1 与 splat LOD 分层搭三维底座；能用 DeepSeek Harness（下称 DSH）批量造数、生成回归测试、定位日志异常时段。 前置知识：Python 3.10+、numpy、OpenCV 基础；知道旋转矩阵与欧拉角的关系；听过 WGS84、ECEF、RTK、POS 这几个词；会写最简单的 bash。

3.1 空地一体：概念与坐标体系

3.1.1 一句话定义：同基准才谈得上融合

空地一体不是「把无人机视频接到大屏上」。它的定义是：把空中平台（无人机、自动机场、吊舱）与地面设施（固定相机、RTK 基站、GIS 底座、指挥席）的数据，放进统一的时空框架里做配准、融合与联动。拆开只有三件事：

1. 统一坐标：像素、机体系、局部 ENU、ECEF、经纬度之间可逆互转，误差可量化。
2. 统一时间：UTC、GPS 周秒、PTS、设备上电计数能互相对齐，对齐残差有统计。
3. 统一度量：GSD、定位精度、时间误差统一到「米」和「秒」，而不是「差不多清楚」「基本同步」。

缺任何一条，后端的融合都会变成玄学：大屏上视频飘、航线偏、降落点抖，而你无法判断是标定错了、时间错了还是底图基准错了。

                 ┌─────────────────────────────┐
                 │  空基：无人机 / 机场 / 吊舱   │
                 │  POS · IMU · 视频 · 载荷数据  │
                 └──────────────┬──────────────┘
                                │ 坐标系：ECEF / ENU / 机体系
┌───────────────────────────────▼───────────────────────────────┐
│ 时空基准层（本章重点，必须先做对）                              │
│  · 坐标：WGS84 ↔ CGCS2000 ↔ 高斯投影 ↔ 地方独立系 ↔ ENU        │
│  · 时间：UTC ↔ GPS 周秒 ↔ 视频 PTS ↔ 飞控上电计数               │
│  · 度量：GSD / 定位精度 / 时间误差 → 统一折算成米               │
└───────────────────────────────┬───────────────────────────────┘
                                │ 统一后的 (t, lon, lat, h, rpy)
┌───────────────────────────────▼───────────────────────────────┐
│ 地基：GIS 底座 · 三维场景 · 成果库 · 指挥席                     │
└───────────────────────────────────────────────────────────────┘

把融合按「在流水线哪一段汇合」分成四层，各层的验收口径完全不同：

大多数团队一上来就做「决策层大屏」，结果每一层都在给基准层还债。正确的顺序是反的：先把基准层做对，再做数据层，特征层和决策层才有意义。

3.1.2 七个坐标系与一条转换链

同一个空间点，在一条数据链路上会以七种身份出现。记住这张表，比记住任何公式书都管用：

转换链必须一条走通，中间任何一跳搞错都会在最后表现为「偏移几十米」：

像素 (u,v)
  │ ① 去畸变 + 内参逆：x_cam = K^-1 · [u, v, 1]^T   (归一化像平面)
  ▼
相机坐标系 (x 右, y 下, z 前)
  │ ② 外参：T_body_cam（手眼标定得到，含 R 与 t）
  ▼
机体坐标系 (FRD: 前右 下 / FLU: 前左 上，注意符号)
  │ ③ 姿态：R_enu_body = Rz(yaw)·Ry(pitch)·Rx(roll)
  ▼
ENU 局部切平面 (原点为 home 点)
  │ ④ 平移：加原点 ECEF，再用 ENU 旋转矩阵转 ECEF
  ▼
ECEF (X, Y, Z)
  │ ⑤ 椭球反算（Bowring 或迭代）
  ▼
经纬高 (lon, lat, h)

用矩阵写出来只有一句：T_enu_pixel = T_enu_body · T_body_cam · T_cam_pixel。难的不是乘法，而是每个矩阵的方向约定和时间戳归属——这个矩阵描述的是哪一时刻的姿态，比它算得准不准更致命。

坑 1：把「像素到经纬度」当成一个二维单应。只要相机有俯仰、地面有起伏（>2 m），单应就会在画面边缘产生数米误差。斜视场景必须走三维射线与地面求交，不能用单应近似。

坑 2：机体坐标系方向搞反。PX4 与 ArduPilot 的姿态输出是 FRD（z 向下），ROS/OpenCV 习惯 FLU 与「z 向前」。差一个符号就是「飞机往左倾时画面往上飘」。

坑 3：把 EPSG:4326 的轴序当 (x, y)。OGC 规范里 4326 是 lat, lon，而大多数代码库默认 lon, lat。pyproj 里务必显式写 always_xy=True，否则整条链路静默转置。

3.1.3 时间基准：0.5 s 的错位就是 5 m 的位移

坐标对齐之后，第二个必须钉死的是时间。先把常见时间源列清楚：

数量级换算一下：多旋翼巡检速度 10 m/s，0.5 s 的时间错位就是 5 m 的位置错位；测绘航线 20 m/s 时，100 ms 错位等于 2 m。所以时间对齐的目标不是「看起来同步」，而是「残差有数、漂移可测」。

用最小二乘同时估出「偏移」和「漂移率」，比手工改偏移靠谱得多：

# time_align.py —— 用一次飞行的视频 PTS 与飞控 UTC 时间戳做线性拟合
import numpy as np

# 真实场景：t_video 来自 ffprobe 逐帧 pts_time，t_utc 来自飞控日志的 GPS 时间
# 这里用合成数据复现「25 fps + 30 ppm 时钟漂移 + 0.8 s 固定偏移 + 8 ms 抖动」
t_video = np.arange(0.0, 300.0, 0.04)
rng = np.random.default_rng(0)
t_utc = 1735689600.0 + t_video * (1 + 30e-6) - 0.8 + rng.normal(0.0, 0.008, t_video.size)

k, b = np.polyfit(t_video, t_utc, 1)
resid = t_utc - (k * t_video + b)
print(f"scale   = {k:.9f}")
print(f"drift   = {(k - 1) * 1e6:+.1f} ppm")
print(f"offset  = {b - 1735689600.0:+.3f} s")
print(f"rms     = {resid.std() * 1000:.1f} ms")
print(f"max     = {np.abs(resid).max() * 1000:.1f} ms")
print(f"30ppm 漂移在 300 s 内累计 {30e-6 * 300 * 1000:.0f} ms，约 {30e-6 * 300 * 10:.2f} m（10 m/s）")

输出应接近 drift = +30 ppm、offset = -0.800 s、rms ≈ 8 ms。如果残差里出现「阶梯」，通常是丢帧或时间基换算错了（比如 mp4 的 timescale 是 12800 而不是 90000）；如果残差是单调增大的曲线，说明两边时钟有温漂，需要用分段线性或二次拟合。

验收口径建议直接写进 CI：

3.1.4 用 pyproj 与 numpy 打通转换链

先解决投影。国内项目最常踩的坑是「WGS84 直接当 CGCS2000 用」和「地方独立坐标系没有参数只有四参数文本」。前者在 cm 级任务（视觉降落、形变监测）不能忽略，后者必须把参数变成 PROJ 字符串或写进 proj.db。

# crs_convert.py —— WGS84 ↔ CGCS2000 高斯投影（3 度带，中央子午线 114E）
from pyproj import Transformer, CRS

WGS84 = CRS.from_epsg(4326)          # 经纬度，度
CGCS2000 = CRS.from_epsg(4490)       # 经纬度，度
GK_CM114 = CRS.from_epsg(4547)       # CGCS2000 / 3-degree Gauss-Kruger CM 114E

to_gk = Transformer.from_crs(WGS84, GK_CM114, always_xy=True)
to_wgs = Transformer.from_crs(GK_CM114, WGS84, always_xy=True)

lon, lat = 113.900000, 22.500000
x, y = to_gk.transform(lon, lat)     # x = 东坐标，y = 北坐标，单位 m
lon2, lat2 = to_wgs.transform(x, y)
print(f"GK      : E={x:.3f} N={y:.3f}")
print(f"round   : {lon2:.9f} {lat2:.9f}")
print(f"d_lon   : {abs(lon2 - lon) * 1e9:.3f} 纳度")

# WGS84 与 CGCS2000 的差异：两系都对齐到 ITRF 系列，国内范围差异在 cm 级
t_2k = Transformer.from_crs(WGS84, CGCS2000, always_xy=True)
print(f"WGS84==CGCS2000 差异（同一点）: {t_2k.transform(lon, lat)}")

注意 EPSG:4534~4554 是「CM」系列，东坐标不含带号前缀（500 km 假东之后的 6 位数）；而 4513~4533 是「zone」系列，东坐标带两位带号。两者看似只差一个数字，混用会让整景数据平移 30 km 以上。带号与 CM 的对应关系可以用一条命令验证：

# 查看某个 EPSG 的 PROJ 字符串与适用范围
projinfo EPSG:4547 -o PROJ
projinfo EPSG:4547 -o WKT2:2019 | head -20
# 打印该投影带内一点的实际坐标
echo "113.9 22.5" | cs2cs EPSG:4326 EPSG:4547 -f "%.3f"

再解决局部切平面。ENU 是所有避障、编队、降落控制的公共语言，必须自己实现一遍并做往返校验，不要依赖「某个库应该有」：

# enu.py —— 大地坐标 ↔ ECEF ↔ ENU，全部用 numpy 手写，便于单元测试
import math
import numpy as np

A = 6378137.0                 # WGS84 长半轴
F = 1 / 298.257223563         # 扁率
E2 = F * (2 - F)              # 第一偏心率平方

def geodetic_to_ecef(lon, lat, h):
    lam, phi = math.radians(lon), math.radians(lat)
    sin_phi = math.sin(phi)
    n = A / math.sqrt(1 - E2 * sin_phi * sin_phi)
    x = (n + h) * math.cos(phi) * math.cos(lam)
    y = (n + h) * math.cos(phi) * math.sin(lam)
    z = (n * (1 - E2) + h) * sin_phi
    return np.array([x, y, z])

def ecef_to_geodetic(xyz, tol=1e-12, max_iter=15):
    x, y, z = xyz
    lam = math.atan2(y, x)
    p = math.hypot(x, y)
    phi = math.atan2(z, p * (1 - E2))
    h = 0.0
    for _ in range(max_iter):
        n = A / math.sqrt(1 - E2 * math.sin(phi) ** 2)
        h = p / math.cos(phi) - n
        phi_new = math.atan2(z, p * (1 - E2 * n / (n + h)))
        if abs(phi_new - phi) < tol:
            phi = phi_new
            break
        phi = phi_new
    n = A / math.sqrt(1 - E2 * math.sin(phi) ** 2)
    h = p / math.cos(phi) - n
    return math.degrees(lam), math.degrees(phi), h

def enu_rotation(lon0, lat0):
    lam, phi = math.radians(lon0), math.radians(lat0)
    sl, cl, sp, cp = math.sin(lam), math.cos(lam), math.sin(phi), math.cos(phi)
    return np.array([
        [-sl, -sp * cl, cp * cl],
        [cl, -sp * sl, cp * sl],
        [0.0, cp, sp],
    ])   # 每行对应 E / N / U 方向在 ECEF 中的分量

def geodetic_to_enu(lon, lat, h, lon0, lat0, h0):
    d_ecef = geodetic_to_ecef(lon, lat, h) - geodetic_to_ecef(lon0, lat0, h0)
    return enu_rotation(lon0, lat0) @ d_ecef

def enu_to_geodetic(e, n, u, lon0, lat0, h0):
    d_ecef = enu_rotation(lon0, lat0).T @ np.array([e, n, u])
    return ecef_to_geodetic(geodetic_to_ecef(lon0, lat0, h0) + d_ecef)

if __name__ == "__main__":
    lon0, lat0, h0 = 113.9, 22.5, 20.0
    lon, lat, h = 113.901, 22.5008, 120.0
    e, n, u = geodetic_to_enu(lon, lat, h, lon0, lat0, h0)
    print(f"ENU     : E={e:.3f} N={n:.3f} U={u:.3f}")
    lon2, lat2, h2 = enu_to_geodetic(e, n, u, lon0, lat0, h0)
    print(f"back    : {lon2:.9f} {lat2:.9f} {h2:.4f}")
    print(f"残差    : {abs(lon2 - lon) * 1e9:.3f} 纳度, {abs(h2 - h):.6f} m")

往返残差应该在 1e-8 度量级（即毫米以下），如果跑出来是「米」，先看 enu_rotation 的转置用错了没有：从 ECEF 到 ENU 是 R·d，从 ENU 到 ECEF 是 R^T·d。

还剩一个高频事故点：高程基准。上面的 h 是椭球高，而国内地图、DEM、机场标高用的是 1985 国家高程基准的正常高，两者差一个高程异常（大地水准面差距），国内量级在 -10 m 到 -40 m 之间，必须用 EGM96/EGM2008 或省级似大地水准面精化模型换算。换算错了不会报错，只会让飞机「按地图高度在山顶下方 20 m 飞」。

# 用 PROJ 做一次大地水准面换算（示例：WGS84 椭球高 → EGM96 正常高）
# 需要 proj-data 提供 geoid 网格；国内高精度请换用省级精化模型
echo "113.9 22.5 120.0" | cs2cs EPSG:4979 EPSG:4326+5773 -f "%.3f"

坑 4：always_xy 漏写。pyproj 的 Transformer 默认按 CRS 定义轴序（4326 是 lat, lon），漏写等于把经纬度转置，投影结果还能算出来，只是全错。

坑 5：ENU 原点随意选。原点必须选在任务区中心（或 home 点）并在整条链路上固定，中途换原点等于把所有历史轨迹平移。

坑 6：把 ENU 的 U 当成海拔。U 是相对原点的椭球高差；要报给业务方的「海拔」得先加原点高程，再扣高程异常。

3.1.5 实景三维中国与地方基准（2026）

前面四小节讲的是「怎么转」。这一节讲「转过去之后信谁」：项目里真正决定精度的不是你的转换代码，而是底图本身。2026 年的现实口径是这样——实景三维中国已经从大规模建设转入「十五五」提质增效阶段，国、省、市、县四级协同建设，全国 DEM 在东部与西部重点区域优于 2 m 格网、其余优于 5 m 格网，重要地理实体按年度增量更新。这三句话决定了三件事：公共底座是分级精度的、更新是按年计的、成果形式是分层的。进展口径以自然资源部公开信息为准，写方案时注明引用年份。

把「底图从哪来」列成一张能直接拿去做技术方案的对照表：

这张表最容易误读的是第一行：「2 m 格网」说的是采样间距，不是高程中误差。2 m 格网的 DEM 在丘陵与城区，高程中误差做到 0.3~1 m 是常态，拿它做贴地飞行规划，安全裕度必须按最坏值加，而不是按「2 m 精度」加。

承诺精度要按业务场景写，不要按「我们用了实景三维」写：

地方基准是第二个容易翻车的地方。每个项目开工前必须把五样东西拿到手并写进项目的基准说明：地方坐标系统名与中央子午线、四参数或七参数及其符号约定、用的是 1985 国家高程基准还是地方高程系、似大地水准面精化模型的版本、CORS/网络 RTK 的源节点与覆盖范围。少一样，后面就会有人用「差不多」的参数凑出一套偏移几十米的成果。

底图到手先做一次「户口检查」，用 GDAL 把 CRS、格网、范围与元数据打印出来，不要看文件名猜：

# 底图三件套：CRS、格网、元数据（-json 输出便于进 CI 做版本比对）
gdalsrsinfo -o proj4 base/dem_2m.tif
gdalinfo -stats -json base/dem_2m.tif | head -40

# 地方四参数（Δx, Δy, 旋转 θ, 尺度 k）写成 PROJ 串再校验，别只存 Excel
# +theta 单位与旋转方向（position_vector / coordinate_frame）由测绘主管单位给定
projinfo "+proj=helmert +x=-1234.567 +y=2345.678 +s=1.2 +theta=0.000573 +convention=position_vector"

再给底图算一次「已知点都差多少」，把绝对精度测出来，而不是靠宣称：

# base_check.py —— 用已知点检查底图/成果的平面偏移，输出残差表
# known_points.csv 列：lon,lat,e_local,n_local（已知点来自 CORS 静态观测或实测控制点）
import csv

import numpy as np
from pyproj import Transformer

T = Transformer.from_crs("EPSG:4326", "EPSG:4547", always_xy=True)  # 换成项目的地方 EPSG

rows = []
with open("known_points.csv", newline="", encoding="utf-8") as f:
    for r in csv.DictReader(f):
        e, n = T.transform(float(r["lon"]), float(r["lat"]))
        rows.append((float(r["e_local"]) - e, float(r["n_local"]) - n))

d = np.array(rows)
print(f"点数 = {len(d)}")
print("平面残差 dE/dN (m)：")
print(np.array2string(d, precision=3, suppress_small=True))
print(f"RMS 平面偏差 = {np.sqrt((d ** 2).sum(axis=1).mean()):.3f} m（RTK 条件下验收 < 0.05 m）")
print(f"最大偏差     = {np.abs(d).max():.3f} m")

已知点残差是唯一能写进验收报告的数字：它同时暴露底图误差、地方参数错误与转换符号错误。残差在米级时先查参数，降到厘米级才轮到讨论底图本身的精度。

注意：① 把「2 m 格网」当「2 m 精度」。格网间距是采样间隔，高程中误差要看成果报告；两者混用会直接把安全裕度算小。② 拿国家成果名义承诺 1:500。国家层面的实景三维成果不是按 1:500 图根精度验收的，1:500 只能自采 + 控制点 + 独立检查点。③ 地方参数只存文本、不带符号约定。四参数的旋转方向与单位（弧度/秒/ppm）错一位，成果整体平移或旋转，而且照样能算出数、不报错。

地方基准怎么入库、怎么在服务层统一、怎么与 3.1.2 的七坐标系对齐，见 第08章 空间数据底座与时空基准。

3.2 占用网格与航线规划

3.2.1 地图表示：先算内存，再谈算法

规划的第一件事不是选算法，是选地图表示。同样的 1 km² 任务区，换一种表示，内存能差两个数量级：

内存不用真跑一遍才知道：按「格点数 × 每格字节数」手算即可。1 km² 在 0.1 m 分辨率下就是 1 亿格，单字节也是 100 MB 量级，所以选表示之前先算这一笔。

结论很直接：0.1 m 的二维栅格只适合几百米见方的局部避障，一个 km² 的巡检任务应该用 0.5~1 m 的 2.5D 代价栅格做全局规划，局部再用 0.1~0.2 m 的局部栅格做重规划。三维占用网格只在「必须穿过建筑/树林间隙」时才值得上，否则纯属给飞控电脑加负载。

3.2.2 代价地图：把「不能飞」和「不想飞」分开

代价地图（costmap）的关键是把两类约束表达清楚：

• 硬约束（blocked）：禁飞区、限高面以下、障碍物物理安全半径内——这些格点直接不可通行。
• 软代价（cost）：离障碍太近、坡度过大、逆光拍摄、能耗高——这些格点可通行但代价高。

最常用的合成公式：

cost(n) = w_obs   · C_obs(n)      障碍斥力（距障碍越近越大，0~1）
        + w_nfz   · C_nfz(n)      禁飞区/限高（软约束部分，硬约束另用 blocked 标）
        + w_slope · C_slope(n)    坡度能耗（坡度/30° 截断到 1）
        + w_energy· C_energy(n)   电池余量折算的返航风险
        + w_quality·C_quality(n)  视觉质量（太阳高度角、GSD 偏差、拼接重叠）
path_cost   = Σ cost(n) × 步长(n)

不同任务给不同的权重，别一套参数打天下：

生成 2.5D 代价栅格的最小实现（依赖 numpy + scipy）：

# costmap.py —— 障碍膨胀 + 禁飞区 + 坡度，输出 cost 与 blocked 两张同尺寸栅格
import numpy as np
from scipy.ndimage import distance_transform_edt, sobel

def build_costmap(dem, occ_mask, nfz_mask, res_m, cfg):
    """
    dem      : 地面高程栅格 (m)，用于算坡度
    occ_mask : 障碍占位布尔栅格（建筑、树冠、塔吊）
    nfz_mask : 禁飞区布尔栅格
    res_m    : 栅格分辨率 (m)
    cfg      : 权重字典，如 {"safety_m":1.5,"w_obs":1.0,"w_slope":0.3}
    """
    safety = cfg["safety_m"]
    dist_m = distance_transform_edt(~occ_mask) * res_m     # 自由格到最近障碍的距离
    c_obs = np.clip(1.0 - dist_m / safety, 0.0, 1.0)        # 越近越大，0 表示足够远
    blocked = (dist_m < safety * 0.3) | nfz_mask            # 物理安全半径 + 禁飞区

    # Sobel 归一化系数 8 为 3x3 核常数，再除以分辨率得到无量纲坡度
    gx = sobel(dem.astype(np.float32), axis=1) / (8.0 * res_m)
    gy = sobel(dem.astype(np.float32), axis=0) / (8.0 * res_m)
    slope = np.hypot(gx, gy)
    c_slope = np.clip(slope / 0.5, 0.0, 1.0)                # 坡度 0.5（约 26.6 度）即满代价

    cost = cfg["w_obs"] * c_obs + cfg["w_slope"] * c_slope
    cost[blocked] = np.inf
    return cost, blocked

# 自测：一块 200x200 栅格，中心立一栋 20x20 的建筑
if __name__ == "__main__":
    dem = np.zeros((200, 200), dtype=np.float32)
    occ = np.zeros((200, 200), dtype=bool)
    occ[80:120, 90:110] = True
    nfz = np.zeros((200, 200), dtype=bool)
    nfz[:, 180:] = True
    cost, blocked = build_costmap(dem, occ, nfz, 0.5, {"safety_m": 1.5, "w_obs": 1.0, "w_slope": 0.3})
    print("blocked 格点数:", int(blocked.sum()), "/", blocked.size)
    print("cost 范围:", float(np.nanmin(cost)), float(np.nanmax(cost)))

注意 cost[blocked] = np.inf 这个写法：把硬约束塞进同一张图，A* 只需要判断 np.isfinite(cost[n])，不用再维护第二张掩膜。代价是不能拿它做插值或梯度——ESDF 要单独保留距离场。

坑 7：膨胀半径按「旋翼半径」算。实际安全裕度应该是旋翼半径 + RTK 水平误差（固定解 2~3 cm，浮动解 0.3~0.5 m）+ 阵风位移（5 级风约 0.3 m）+ 地图配准误差（0.2~0.5 m）。1.5 m 是巡检的常用值，不是随便拍的。

坑 8：DEM 与障碍物不是同一时刻。用两年前的 DOM/DEM 规划，新加的光伏板、通信塔、临时脚手架一个都不在图上。工程上要么起飞前做一次正射巡视，要么把「未知区域」默认加代价。

坑 9：坡度用高程差直接相减。相邻格点高差除以格距才是坡度，且必须先确认 DEM 的像元单位是度还是米（有些数据的 CRS 是经纬度，像元大小是 1e-5 度）。

3.2.3 规划算法：A* 打底，D* Lite 重规划，覆盖路径单独一套

A* 的可运行最小实现（八邻域 + 八方向启发式 + 地形代价）：

# astar.py —— 八邻域 A*，步长代价 = 几何距离 x (1 + 地形代价)
import heapq

import numpy as np

def astar(cost, start, goal):
    h, w = cost.shape
    usable = np.isfinite(cost)

    def heuristic(a, b):
        dx, dy = abs(a[0] - b[0]), abs(a[1] - b[1])
        return (dx + dy) + (np.sqrt(2.0) - 2.0) * min(dx, dy)   # 八方向精确启发

    if not (usable[start] and usable[goal]):
        return None

    open_heap = [(heuristic(start, goal), 0.0, start)]
    best = {start: 0.0}
    parent = {}
    closed = set()

    while open_heap:
        _, g_cur, cur = heapq.heappop(open_heap)
        if cur in closed:
            continue
        closed.add(cur)
        if cur == goal:
            path = [cur]
            while cur in parent:
                cur = parent[cur]
                path.append(cur)
            return path[::-1]
        for dx in (-1, 0, 1):
            for dy in (-1, 0, 1):
                if dx == 0 and dy == 0:
                    continue
                nb = (cur[0] + dx, cur[1] + dy)
                if not (0 <= nb[0] < h and 0 <= nb[1] < w) or not usable[nb]:
                    continue
                step = np.hypot(dx, dy) * (1.0 + float(cost[nb]))
                ng = g_cur + step
                if ng < best.get(nb, float("inf")):
                    best[nb] = ng
                    parent[nb] = cur
                    heapq.heappush(open_heap, (ng + heuristic(nb, goal), ng, nb))
    return None

if __name__ == "__main__":
    rng = np.random.default_rng(7)
    grid = rng.random((120, 160)).astype(np.float32) * 0.6
    grid[10:15, 20:140] = np.inf          # 一条横向禁飞带，留出缺口
    grid[10:15, 70:75] = 0.0
    path = astar(grid, (5, 5), (115, 155))
    print("路径点数:", None if path is None else len(path))
    print("前 5 个点:", None if path is None else path[:5])

覆盖率更高的做法是分层规划：全局用 1 m 栅格 A（或 D Lite）出航段，局部用 0.2 m 栅格 + 实时障碍做重规划，中间用「航点 + 速度剖面」衔接。把 A* 直接跑在 0.1 m 全图上，规划一次十几秒，实飞时根本来不及。

D* Lite 的价值在于「地图边走边更新」：它从终点反向搜索，障碍变化时只重算受影响的一段，实测比从头跑 A* 快一到两个数量级。工程上常用的简化是「A* + 局部重规划窗口」——每次只重算前方 50~100 m 的走廊，效果够用，实现量小很多。

覆盖路径的生成只需要三个量：单张旁向覆盖宽度 W = GSD × 短边像素数、航带间距 S = W × (1 − 旁向重叠率)、航点数 = 航带数 × 2。矩形区域按牛耕式逐带扫描，奇数带正向、偶数带反向即可。

三个常见变体：分区覆盖（大区域切成若干凹/凸多边形分别跑，边界处再加一条过渡航带）、螺旋覆盖（适合圆形区域与机场净空区）、沿边覆盖（适合河湖岸线、公路、输电线路）。选哪种取决于区域形状：矩形用牛耕，圆形用螺旋，带状用沿边加往返。

3.2.4 GSD、重叠率、航高：三个数互相锁死

地面分辨率（GSD，Ground Sample Distance）的定义是：影像上一个像元在地面上对应多远。用相似三角形直接推：

相机针孔模型：物方一点 P 与像方成像 p 满足相似三角形
  像元尺寸 p (m)         地面上对应的 GSD (m)
  ────────────────   =   ─────────────────────
  焦距 f (m)              航高 H (m)

=>  GSD = p × H / f             （最常用的一条公式，记牢）
=>  H   = GSD × f / p           （反算：想要某 GSD 该飞多高）
=>  单张覆盖宽 W = GSD × 像元数  （像元数取对应方向的像素数）

算例（1 英寸 20 MP 相机：像元 2.4 μm，焦距 8.8 mm，5472×3648）：

同一块地，GSD 从 2.7 cm 放到 10 cm，飞行时间从 42 min 降到 11 min——这就是为什么「先问成图比例尺，再定 GSD，再定航高」必须按顺序来。成图比例尺与 GSD 的经验换算：

重叠率直接决定空三能不能成功：

航时估算按「有效扫描宽度」算，而不是按航带条数拍脑袋：有效扫描宽度 S = 覆盖宽 × (1 − 旁向重叠率)，纯飞行时间 t ≈ 面积 / (速度 × S)，再按每条航带 8 s 左右加转弯损耗。

坑 10：GSD 用长边像素算。覆盖宽度必须用「垂直航向的那个边」的像素数，用错边会导致航带间距系统性偏大或偏小，成图出现条带状空洞。

坑 11：把相对航高当绝对航高。H 是「相机到地面」的距离，山区要用地形跟随，否则山脊处 GSD 变小、山谷处 GSD 变大，同一景里 GSD 能差 2 倍。

坑 12：忽略相机倾斜。五镜头倾斜摄影里，倾斜相机的 GSD 与正射相机不同，航线设计要按主相机算，倾斜相机只影响建模纹理。

3.2.5 安全裕度与返航策略

航线再漂亮，安全策略没写清楚就是隐患。给一套能直接落地成状态机的规则：

把返航当成一次普通的 A* 规划，只是终点换成 home 点、权重里加大能耗项；电量判断用「返航所需电量 × 1.5 ≤ 剩余电量」这条硬线，其中返航所需电量 = (巡航时间 + 爬升时间) × 单位时间耗电量。

坑 13：限高用海拔判断。120 m 是「真高」概念，山区按海拔判断会直接撞山，必须用地形跟随得到的 AGL。

坑 14：返航高度低于障碍物。返航直线穿越塔吊/高压线是事故高发点，返航路径也要过代价地图。

坑 15：电量模型用标称值。低温、大载重、逆风都会让能耗上升 20%~40%，返航阈值必须留够冗余，宁可提前收工。

3.3 视觉引导自主降落

3.3.1 四种方案对比：精度都是标称的，条件才是真的

降落是整个链路里对精度最敏感的环节：水平偏 1 m 就可能砸到护栏。四种主流方案的实际差异不在「标称精度」，而在「什么条件下还能保持精度」：

工程上的正确做法不是四选一，而是分层降级：远距离靠 RTK 直降到 8 m 以下，5 m 内切视觉标记精调，1 m 内用测距仪与视觉联合判定触地。纯视觉从 50 m 直接降的，十次里总有一两次因为标记太小检不出来。

3.3.2 标定：内参决定下限，手眼决定上限

相机内参标定是整条链路的地基。流程与验收标准要写死：

# 依赖
python -m pip install "opencv-contrib-python>=4.8" numpy

# 采集要求（直接影响能不能达标）
#  · 棋盘格 9x6 内角点，方格边长 25 mm（用卡尺实测，不要信打印标称）
#  · 20~30 张，覆盖画面四角与中心，倾角 15~45 度
#  · 光照均匀，避免反光与运动模糊（相机固定、板动）

# calib_intrinsic.py —— 棋盘格内参标定，输出 RMS 与内参
import glob

import cv2
import numpy as np

PATTERN = (9, 6)              # 内角点数（列, 行），不是方格数
SQUARE_M = 0.025              # 方格实际边长（米）

objp = np.zeros((PATTERN[0] * PATTERN[1], 3), np.float32)
objp[:, :2] = np.mgrid[0:PATTERN[0], 0:PATTERN[1]].T.reshape(-1, 2) * SQUARE_M

obj_points, img_points = [], []
for path in sorted(glob.glob("calib/*.jpg")):
    img = cv2.imread(path)
    gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
    ok, corners = cv2.findChessboardCorners(gray, PATTERN, None)
    if not ok:
        print("跳过（未检出角点）:", path)
        continue
    corners = cv2.cornerSubPix(
        gray, corners, (11, 11), (-1, -1),
        (cv2.TERM_CRITERIA_EPS + cv2.TERM_CRITERIA_MAX_ITER, 30, 0.001))
    obj_points.append(objp)
    img_points.append(corners)

rms, K, dist, rvecs, tvecs = cv2.calibrateCamera(
    obj_points, img_points, gray.shape[::-1], None, None)

print(f"重投影 RMS = {rms:.4f} px   （验收 < 0.5 px，优秀 < 0.3 px）")
print("K =")
print(np.array2string(K, precision=3, suppress_small=True))
print("畸变系数 =", np.array2string(dist.ravel(), precision=4))
print(f"每张图的平均重投影误差: {[f'{np.mean(np.abs(cv2.projectPoints(objp, rvecs[i], tvecs[i], K, dist)[0] - img_points[i])):.3f}' for i in range(min(3, len(rvecs)))]}")

手眼解算直接用 cv2.calibrateHandEye（eye-in-hand 用 TSAI 或 PARK 均可），输入是 N 组「机体在基站系下的位姿」与「标定板到相机的外参」。验收两条：代回全部样本后重投影误差 < 0.5 px、平移残差 < 5 mm；不达标先加姿态样本，别急着换算法。

坑 16：棋盘格打印后直接量纸。打印机会缩放 1%~3%，必须用卡尺量「黑白格交界到交界」的实际边长，否则内参系统性偏差，降落到最后 1 m 时误差被放大。

坑 17：手眼标定只采一组姿态。运动学退化（所有姿态近似共面）时解不唯一，必须做绕轴旋转激励，采集时把无人机架在三轴云台上转一圈最省事。

坑 18：标定与使用不是同一分辨率。标定用 1920×1080，推理时缩放成 640×360，内参必须按比例缩放（c_x、c_y、f_x、f_y 同时乘缩放比），否则 PnP 直接算歪。

3.3.3 solvePnP：四个点为什么够，平面为什么会退化

透视 n 点定位（PnP）解决的问题是：已知若干点的三维坐标与对应像素坐标，求相机相对这些点的位姿。理论最少点数：

平面标记的核心风险是平面位姿二义性：一个平面目标在图像上的投影，可能对应两组位姿（镜面对称解），在正对且距离较远时尤其明显。IPPE 会一次返回两组解，必须用「历史位姿连续性 + 重投影误差」选出真解：

# pnp_landing.py —— 平面标记位姿解算，含二义性消解
import cv2
import numpy as np

MARKER_M = 0.40          # 标记边长（米），实测值

def marker_object_points(size=MARKER_M):
    s = size / 2.0
    return np.array([[-s, -s, 0.0], [s, -s, 0.0], [s, s, 0.0], [-s, s, 0.0]],
                    dtype=np.float64)

def solve_pose(image_points, K, dist, prev_tvec=None):
    objp = marker_object_points()
    ok, rvecs, tvecs, errs = cv2.solvePnPGeneric(
        objp, image_points, K, dist, flags=cv2.SOLVEPNP_IPPE)
    if not ok:
        return None
    cands = []
    for rvec, tvec, err in zip(rvecs, tvecs, errs):
        rvec = rvec.reshape(3, 1)
        tvec = tvec.reshape(3, 1)
        proj, _ = cv2.projectPoints(objp, rvec, tvec, K, dist)
        reproj = float(np.mean(np.linalg.norm(
            proj.reshape(-1, 2) - image_points.reshape(-1, 2), axis=1)))
        # 距离先验：标记必须落在相机前方 0.3~10 m
        cands.append({"rvec": rvec, "tvec": tvec, "reproj": reproj,
                      "z": float(tvec[2])})
    cands = [c for c in cands if 0.2 < c["z"] < 12.0]
    if not cands:
        return None
    if prev_tvec is not None and len(cands) > 1:
        # 用上一帧位置选最近的解，抑制正对时的翻转
        for c in cands:
            c["jump"] = float(np.linalg.norm(c["tvec"].ravel() - np.ravel(prev_tvec)))
        cands.sort(key=lambda c: (c["jump"] > 0.5, c["reproj"]))
    else:
        cands.sort(key=lambda c: c["reproj"])
    return cands[0]

精度怎么估？在深度 Z、焦距 f 下，1 像素的角度误差约等于 arctan(1/f)，横向位置误差近似：

横向误差 ≈ Z / f  （单位：米每像素）
例：Z = 3 m, f = 900 px  =>  1 px 约 3.3 mm
    若重投影误差 1.5 px，则横向不确定度约 5 mm/轴
纵向（深度）误差 ≈ Z² / (f × 标记边长像素数)，随距离平方增长
=> 结论：1 m 内用 1 像素的误差也只有毫米级，10 m 外误差到分米级

所以降落策略必然是「越近越准」：8 m 以上别指望视觉给厘米级，老老实实用 RTK；4 m 以内视觉才真正跑赢 RTK。

坑 19：标记边长填标称值。打印 40 cm 的板子实际可能是 39.6 cm，导致深度整体偏 1%，5 m 处就是 5 cm 的纵向误差，最后表现为「降落时总是偏一点点」。

坑 20：只用 4 个角点。标记四个角点受限于角点检测精度（亚像素约 0.1~0.3 px），更好的做法是用标记内部的多个编码点，或对每个角点做多帧滤波。

坑 21：忽略镜头畸变。广角镜头边缘畸变 5%~10%，PnP 会把畸变当成位姿误差，表现为「标记在画面中心时准，偏到边缘就歪」。

3.3.4 视觉伺服状态机：五段式，每段都要有超时

降落不是一个 PID 控制器能搞定的事，它是一条带安全约束的状态机。五个阶段的进入/退出条件、超时与保护动作如下：

# landing_sm.py —— 五段式视觉伺服状态机骨架（可直接接 ROS2 / MAVSDK）
import time

class LandingSM:
    def __init__(self, kp_xy=2.5, kp_z=1.2, max_vxy=0.8, max_vz=0.6):
        self.state = "SEARCH"
        self.t_enter = time.time()
        self.kp_xy, self.kp_z = kp_xy, kp_z
        self.max_vxy, self.max_vz = max_vxy, max_vz
        self.lost = 0

    def _timeout(self):
        return {"SEARCH": 20.0, "ALIGN": 15.0, "DESCEND": 30.0,
                "FINAL": 20.0, "TOUCHDOWN": 5.0}[self.state]

    def step(self, det, alt_agl, t_now):
        """det: None 或 {"x":..,"y":..,"z":..}（机体系下标记位置，单位 m）"""
        dt = t_now - self.t_enter
        if dt > self._timeout():
            return self._on_timeout()

        if det is None:
            self.lost += 1
            if self.lost > 15 and self.state in ("ALIGN", "DESCEND", "FINAL"):
                self.state = "SEARCH"
                self.t_enter = t_now
            return {"state": self.state, "vz": 0.0 if alt_agl < 6 else -0.2,
                    "vxy": (0.0, 0.0), "note": "标记丢失"}
        self.lost = 0

        err = (det["x"], det["y"])
        vx = max(-self.max_vxy, min(self.max_vxy, self.kp_xy * err[0]))
        vy = max(-self.max_vxy, min(self.max_vxy, self.kp_xy * err[1]))

        if self.state == "SEARCH" and det is not None:
            self._goto("ALIGN", t_now)
        elif self.state == "ALIGN" and abs(err[0]) < 0.15 and abs(err[1]) < 0.15:
            self._goto("DESCEND", t_now)
        elif self.state == "DESCEND" and alt_agl < 1.2:
            self._goto("FINAL", t_now)
        elif self.state == "FINAL" and alt_agl < 0.15:
            self._goto("TOUCHDOWN", t_now)

        # 下降速度随高度收敛：近地更慢，避免触地冲击
        vz = -min(self.max_vz, max(0.15, self.kp_z * alt_agl * 0.5))
        if self.state in ("SEARCH", "ALIGN"):
            vz = -0.2 if self.state == "ALIGN" else 0.0
        return {"state": self.state, "vxy": (vx, vy), "vz": vz}

    def _goto(self, new_state, t_now):
        self.state = new_state
        self.t_enter = t_now

    def _on_timeout(self):
        if self.state == "SEARCH":
            return {"state": "ABORT", "vz": 0.0, "vxy": (0.0, 0.0), "note": "搜索超时返航"}
        return {"state": "HOVER", "vz": 0.0, "vxy": (0.0, 0.0), "note": "超时悬停重检"}

PID 整定的经验顺序：先位置环 P → 再加 D 压振荡 → 最后少量 I 消稳态误差，不要一上来三个参数一起调。

控制输出必须限幅并做一阶低通：视觉检测在 15~30 Hz，飞控在 50~400 Hz，中间不平滑会出现「打摆子」。一阶低通 v_k = a·v_meas + (1-a)·v_{k-1}，a 取 0.2~0.35 通常够用；对抖动更敏感的场景用 One Euro 滤波，能在静止时压噪、运动时保相位。

3.3.5 鲁棒性：失锁、抖动、风扰、光照

真实场景里，视觉降落的失败原因八成不是算错，而是「看不见了」：

视觉降落故障树
├─ 检测不到标记
│  ├─ 光照：逆光 / 强反光 / 夜间无补光 → 加偏振片、遮光罩、主动补光、曝光锁定
│  ├─ 尺度：标记 < 15 px → 放大标记 or 降低起始高度 or 换长焦
│  ├─ 角度：俯仰/滚转 > 30 度 → 标记透视压缩，改云台增稳
│  └─ 遮挡：树叶/电线/吊臂 → 多标记冗余 + RTK 兜底
├─ 检到了但位姿跳变
│  ├─ 平面二义性翻转 → 用上一帧做连续性约束
│  ├─ 角点抖动 → 亚像素 + 多帧卡尔曼平滑
│  └─ 时间错位 → 位姿按视频 PTS 插值，不用「最新姿态」
├─ 控制发散
│  ├─ 增益过高 → 降 kp、加 D、加低通
│  ├─ 风扰 → 前馈 + 积分限幅
│  └─ 传感器延迟 → 缩短链路、预测补偿
└─ 触地误判
   ├─ 测距仪噪声 → 中值滤波 + 连续 N 帧确认
   ├─ 地面杂物 → 融合 IMU 加速度突变
   └─ 地效区气压跳变 → 近地 1 m 内不采信气压高度

失锁处理要分级，不能一丢就爬升：

风扰的处理不是「加大增益」，而是加前馈：把风速估计（或 IMU 的比力变化）作为前馈项加到速度指令上，反馈只负责修正残差。积分项一定要限幅（如 ±0.3 m/s），否则风停之后会「窜一下」。

坑 22：用「最新姿态」而不是「视频帧采集时刻的姿态」。相机曝光到图像送入算法有 50~150 ms 延迟，这个延迟在 0.5 m/s 下降时就是 7.5 cm 的纵向位置差，表现为「越降越偏」。

坑 23：触地判定只看高度。气压计在地效区会跳变 0.5~1 m，正确做法是「测距仪连续 N 帧 < 阈值 + IMU 垂直加速度突变 + 视觉 Z 突然不变」三选二。

坑 24：不设视觉降落的「起降窗口」。夜间、雨天、6 级以上风必须直接禁用视觉，靠 RTK 或人工接管，这些条件要写进状态机的入口判断。

3.4 无人机视频 AR 投影

3.4.1 数学：把像素反投影成地面点

AR 投影要回答的问题是：画面上的第 (u, v) 个像素，对应地面上哪个点？针孔模型下就三步：

1) 像素 -> 相机射线（归一化像平面）
   d_cam = K^-1 · [u, v, 1]^T           K = [[fx,0,cx],[0,fy,cy],[0,0,1]]
2) 相机 -> ENU 方向
   d_enu = R_enu_body · R_body_cam · d_cam
3) 射线与地面求交（地面高程 h_g 已知，或与椭球/瓦片求交）
   t = (h_g - h_cam) / d_enu.z
   P_enu = C_enu + t · d_enu        （C_enu 为相机在 ENU 下的位置）
   若 d_enu.z >= 0 （射线朝上或水平）则无交点，该像素不可投影

这就是所有 AR 投影的共同内核。剩下的是工程取舍：算哪些像素、用哪张地面、怎么把结果画出来。

3.4.2 两条实现路线：贴四边形，还是逐像素投影

路线 A 的关键是记得它假设地面是一个平面。俯仰角越大、地形起伏越大，误差越离谱——因为透视四边形在平面上才能用仿射纹理映射，而真实地面是曲面，四边形贴上去必然拉伸。

Cesium 里路线 A 的最小写法（视频元素 + 视频材质 + 四个角点）：

// ar_overlay.js —— 把 <video> 按四角点贴到地球上
const videoEl = document.getElementById("uav-video");

// 1) 反投影四个角：注意 Cesium 的像素坐标原点是画布左上角
function cornerOnGround(viewer, x, y) {
  const ray = viewer.camera.getPickRay(new Cesium.Cartesian2(x, y));
  if (!ray) return null;
  return viewer.scene.globe.pick(ray, viewer.scene) || null;
}

const corners = [
  cornerOnGround(viewer, 0, 0),
  cornerOnGround(viewer, w, 0),
  cornerOnGround(viewer, w, h),
  cornerOnGround(viewer, 0, h),
];
if (corners.some((c) => !c)) {
  console.warn("有角点射线未命中地球，检查相机姿态或降低俯仰角");
}

// 2) 用 VideoMaterialProperty 贴到矩形；矩形只适合「近似正射」的场景
viewer.entities.add({
  rectangle: {
    coordinates: Cesium.Rectangle.fromDegrees(lonMin, latMin, lonMax, latMax),
    material: new Cesium.VideoMaterialProperty(videoEl, {
      loop: true,
      autoPlay: true,
      color: Cesium.Color.WHITE.withAlpha(0.85),
    }),
    classificationType: Cesium.ClassificationType.TERRAIN,
  },
});

注意上面的矩形写法本身仍是「轴对齐矩形」，只适合小范围且相机近似垂直向下。真正的透视四边形必须自己建网格（路线 B 的简化版）：把四边形切成 N×N 个小三角形，每个顶点独立求交，UV 用双线性插值，这样透视关系才正确。纹理矩阵（TextureMatrix）的本质就是把这个从 UV 到纹理坐标的映射作为 uniform 传给着色器；具体 API 名称随 Cesium 版本变化，写代码前用当前版本文档核对，不要照抄旧博客。

本章的 Cesium 口径统一按 CesiumJS 1.147（2026-11 发布） 核对。这个版本除了常规修复，还补了依赖与内存问题——对「连续渲染几天不出错」的指挥大屏来说，锁到一个已验证的小版本比追最新更划算。VideoMaterialProperty、TextureMatrix、Cesium3DTileset 的构造参数与属性名在 1.x 里改过多次，升级前先把用到的 API 对一遍 1.147 文档，再动代码。

3.4.3 坐标纠偏：GCJ-02 只在展示层出现

国内底图（高德、腾讯、部分天地图切片）用的是 GCJ-02，而无人机 RTK 输出是 WGS84，两者在同一位置能差几百米。这不是「精度问题」，是坐标系不同，必须显式转换：

# gcj02.py —— WGS84 -> GCJ-02（公开的近似算法，误差约 1~2 m）
import math

A = 6378245.0
EE = 0.00669342162296594323

def _out_of_china(lng, lat):
    return not (72.004 <= lng <= 137.8347 and 0.8293 <= lat <= 55.8271)

def _transform_lat(x, y):
    ret = -100.0 + 2.0 * x + 3.0 * y + 0.2 * y * y + 0.1 * x * y + 0.2 * math.sqrt(abs(x))
    ret += (20.0 * math.sin(6.0 * x * math.pi) + 20.0 * math.sin(2.0 * x * math.pi)) * 2.0 / 3.0
    ret += (20.0 * math.sin(y * math.pi) + 40.0 * math.sin(y / 3.0 * math.pi)) * 2.0 / 3.0
    ret += (160.0 * math.sin(y / 12.0 * math.pi) + 320.0 * math.sin(y * math.pi / 30.0)) * 2.0 / 3.0
    return ret

def _transform_lng(x, y):
    ret = 300.0 + x + 2.0 * y + 0.1 * x * x + 0.1 * x * y + 0.1 * math.sqrt(abs(x))
    ret += (20.0 * math.sin(6.0 * x * math.pi) + 20.0 * math.sin(2.0 * x * math.pi)) * 2.0 / 3.0
    ret += (20.0 * math.sin(x * math.pi) + 40.0 * math.sin(x / 3.0 * math.pi)) * 2.0 / 3.0
    ret += (150.0 * math.sin(x / 12.0 * math.pi) + 300.0 * math.sin(x / 30.0 * math.pi)) * 2.0 / 3.0
    return ret

def wgs84_to_gcj02(lng, lat):
    if _out_of_china(lng, lat):
        return lng, lat
    dlat = _transform_lat(lng - 105.0, lat - 35.0)
    dlng = _transform_lng(lng - 105.0, lat - 35.0)
    radlat = lat / 180.0 * math.pi
    magic = 1 - EE * math.sin(radlat) ** 2
    sqrtmagic = math.sqrt(magic)
    dlat = (dlat * 180.0) / ((A * (1 - EE)) / (magic * sqrtmagic) * math.pi)
    dlng = (dlng * 180.0) / (A / sqrtmagic * math.cos(radlat) * math.pi)
    return lng + dlng, lat + dlat

if __name__ == "__main__":
    print(wgs84_to_gcj02(116.404, 39.915))   # 北京，偏移量约 0.005 度量级

反向（GCJ-02 → WGS84）没有解析解，用迭代逼近，一般 3~5 次就收敛到 1e-9 度。工程原则只有一条：

坑 25：把 GCJ-02 转过的坐标写回数据库。二次转换会累积 1~2 m 误差，且失去可逆性；正确做法是只在前端渲染前转。

坑 26：以为「偏移是固定的」。GCJ-02 偏移量随位置非线性变化，国内不同城市方向与大小都不同，不能用一个常数平移。

坑 27：AR 投影用底图坐标算相机位置。相机位姿来自 WGS84，底图是 GCJ-02，两者混算的结果就是视频整体平移几百米还带旋转。

3.4.4 遮挡、拉伸与边缘：三个必须限制的参数

限制俯仰角是性价比最高的一条：绝大多数「投影难看」的问题，都是相机俯视角度太小导致的数值发散，把显示条件卡死在 70 度以内，问题就少一大半。

3.4.5 3D Tiles 1.1 与 splat LOD：三维底座的两种格式

巡检项目的三维地球从来不是「一团模型」。它要同时装两类东西：一类是测绘成果——倾斜模型、点云、DEM、BIM/CAD，要求可量测、带属性、有层级；另一类是近景真实感——塔基、站房、设备、走廊，要求「看着像现场」。2026 年的工程答案是把它们拆成两种格式、两层底座，而不是硬塞进同一种数据里。

第一种：3D Tiles 1.1。 相对 1.0，工程上真正用得着的是两个变化：

• metadata（要素级元数据）：构件的编号、材质、电压等级、上次巡检结论可以直接写在瓦片里，前端点选即可读出，不用再维护一套「属性服务 + ID 关联」。
• implicit tiling（隐式瓦片）：用模板 URI 加 subtree 文件描述四叉/八叉层级，tileset.json 从几十 MB 降到几百 KB，加载器按需推导子瓦片。千万级瓦片才谈得上运维。
• 另外 1.1 允许 glTF 直接作为瓦片内容、一个瓦片挂多份内容；旧的 b3dm/i3dm/pnts 仍能读，但新数据建议直接按 glTF + 结构化元数据组织。

第二种：带层次 LOD 的 3D Gaussian Splats。 Cesium 在 2026-04 推出了「3D Tiles + 层次 LOD」的 splat 方案，CesiumJS 1.147 可以加载带 LOD 的 splat tileset。要点是：splat 不能当一个大 .ply 整体扔进场景，必须按空间切块、每块做多级细节（不同点数与分辨率），运行时按视距切。原因很物理——加载器每一帧都要对可见高斯排序，点数上去了 GPU 排序开销是硬成本，LOD 不是可选优化。

两种格式的分工：

加载骨架如下（属性名以 CesiumJS 1.147 文档为准；splat 的构造入口与发布格式强相关，换格式先查文档）：

// base_dual.js —— 双底座：测绘成果走 3D Tiles 1.1，近景走带 LOD 的 splat
const viewer = new Cesium.Viewer("cesiumContainer", {
  terrain: Cesium.Terrain.fromWorldTerrain(),
  animation: false,        // 指挥大屏常年挂机，关掉用不到的动画/时间轴
  timeline: false,
});

async function loadDualBase(v) {
  // 1) 测绘成果：SSE 决定瓦片调度密度，显存给足但别指望无限加
  const survey = await Cesium.Cesium3DTileset.fromUrl(
    "https://example.com/tiles/survey/tileset.json",
    {
      maximumScreenSpaceError: 16,   // 默认 16；判读时可降 8，代价是带宽与显存
      maximumMemoryUsage: 1024,      // MB，默认 512
      dynamicScreenSpaceError: true, // 高空俯瞰时放宽，减少无效瓦片
    });
  v.scene.primitives.add(survey);

  // 2) 近景真实感：splat tileset 同样按瓦片调度，不要整包加载
  const splat = await Cesium.Cesium3DTileset.fromUrl(
    "https://example.com/tiles/splat/tileset.json",
    { maximumScreenSpaceError: 8, maximumMemoryUsage: 768 });
  v.scene.primitives.add(splat);
  return { survey, splat };
}

loadDualBase(viewer).then(({ survey, splat }) => {
  // 3) 两块底座必须落在同一套 CRS/ECEF 上：先对齐再叠加
  console.log(survey.boundingSphere, splat.boundingSphere);
});

性能问题基本都集中在这几处：

注意：① metadata 不是「随便塞属性」。属性名与类型要按结构化元数据 schema 定义并随 tileset 发布，前端才能按 key 查询；把自定义 JSON 直接塞进瓦片，加载器读不出来。② splat 的几何是高斯拟合，不能当测绘成果用——要量测、要角点、要净空分析，回 3D Tiles 的显式几何，splat 只负责「看起来对」。③ 位姿没做米制恢复就叠 splat，近景会整体偏大或偏小，贴到塔基上出现缝隙或穿插，第一反应常误判成「模型精度不够」。

3.4.6 视频贴地 + 模型属性联动（CesiumJS 1.147）

3.4.2 的路线 A/B 解决「视频贴在哪里」，这一节解决「贴上去之后能点开什么」。巡检指挥席上最常见的交互是：视频按真视角铺在塔基附近的地面上，用户点一下画面里的塔，右侧弹出这座塔的属性——电压等级、所属线路、上次巡检结论，再顺手调出对应机位的实时视频。这要求三维场景里的模型不只是「几何」，还要带「要素」与「属性」。

CesiumJS 1.147 补上的一块是 Cesium ion Element Service 支持：BIM/CAD 模型的属性查询与计算、BIM 风格的构件树视图。把它和 3.4.5 的 3D Tiles 1.1 metadata 放在一起看，边界很清楚——metadata 解决「属性存在瓦片里、离线可查」，Element Service 解决「属性怎么组织成可浏览的构件树、怎么做计算」。两条路怎么选：

联动代码的骨架只有三步：拾取要素、读属性、按业务 ID 查数据。

// link_video_model.js —— 点选模型构件 -> 读属性 -> 联动视频与工单面板
const handler = new Cesium.ScreenSpaceEventHandler(viewer.scene.canvas);

handler.setInputAction((click) => {
  const picked = viewer.scene.pick(click.position);
  if (!Cesium.defined(picked)) return;

  if (picked instanceof Cesium.Cesium3DTileFeature) {
    // 1) 3D Tiles 1.1：要素属性直接从瓦片 metadata 读
    const props = {};
    picked.getPropertyIds().forEach((id) => { props[id] = picked.getProperty(id); });

    // 2) 业务联动：用构件/设备编号去查视频通道、工单、台账
    const deviceId = props.device_id || props.element_id;
    console.log("构件属性:", props, "业务 ID:", deviceId);
    // fetch(`/api/devices/${encodeURIComponent(deviceId)}/streams`)
  }
}, Cesium.ScreenSpaceEventType.LEFT_CLICK);

三个实现约束：

1. 拾取层级要定义清楚。视频贴地一般用分类（ClassificationType.TERRAIN）挂在地形上，3D Tiles 模型在它上面；同一屏幕坐标可能同时命中视频 entity、地形和模型。业务上要先判模型、再判视频，并把「点视频点不到模型」当成可复现的 bug 修，而不是让用户换个角度点。
2. 属性查询不要逐帧做。一次点击查一次是毫秒级；但「视频里框出 50 个目标同时查属性」要改用空间索引或 metadata 批量查询，否则每帧几十次拾取，帧率直接掉一半。
3. 关联字段要版本化。device_id、element_id 这类关联键写进数据契约（见 3.6.1 的字段规范），模型侧与业务侧字段名对不上时，联调时间比开发还长。

注意：① Element Service 是 ion 侧能力，内网项目不能默认它可用；私有化场景把属性写进 3D Tiles 1.1 metadata，前端自己建树。② 视频贴地 entity 与模型的深度关系受 classificationType、depthTest 影响，贴地不等于「盖在模型上」，要叠加先确认渲染顺序，别用透明度硬凑。③ Cesium3DTileFeature 的方法名与 metadata 读取方式在 1.x 里调整过，升级按 1.147 文档核对，别只改版本号。

3.5 视频轨迹与实时拼接

3.5.1 轨迹与拼接是两件事

新手最常混的一对概念。它们输入相似、输出完全不同：

把这两件事塞进一个服务里，最后一定会因为「轨迹要秒级刷新」和「拼接要滑窗优化」抢资源而崩掉，拆成两个服务更省事。

3.5.2 轨迹：先对齐时间，再抽稀，最后才渲染

轨迹链路三步：POS 与视频对齐（3.1.3 已解决）→ 抽稀 → 三维渲染。抽稀必须在局部平面坐标里做，直接对经纬度做 Douglas-Peucker 会因为经线收敛而失真：

# track.py —— 投影到 ENU 后做 Douglas-Peucker 抽稀
import numpy as np

def perpendicular_distance(p, a, b):
    if np.allclose(a, b):
        return float(np.linalg.norm(p - a))
    return float(np.abs(np.cross(b - a, p - a)) / np.linalg.norm(b - a))

def douglas_peucker(points, tol):
    """points: (N, k) numpy 数组，tol 为容差（米）"""
    if len(points) < 3:
        return points
    keep = np.zeros(len(points), dtype=bool)
    keep[0] = keep[-1] = True
    stack = [(0, len(points) - 1)]
    while stack:
        i, j = stack.pop()
        if j - i < 2:
            continue
        seg = points[i + 1:j]
        d = np.array([perpendicular_distance(p, points[i], points[j]) for p in seg])
        k = int(np.argmax(d))
        if d[k] > tol:
            idx = i + 1 + k
            keep[idx] = True
            stack += [(i, idx), (idx, j)]
    return points[keep]

if __name__ == "__main__":
    # 合成一段 10 Hz、带 RTK 抖动的圆形航迹
    t = np.arange(0, 120, 0.1)
    r = 80.0
    xy = np.c_[r * np.cos(t / 20.0), r * np.sin(t / 20.0)]
    xy += np.random.default_rng(1).normal(0, 0.03, xy.shape)   # 3 cm 噪声
    for tol in (0.1, 0.5, 1.0, 2.0):
        keep = douglas_peucker(xy, tol)
        print(f"tol={tol:4.2f} m 保留 {len(keep):4d} / {len(xy)} 点 ({len(keep) / len(xy):.1%})")

经验参数：航迹回放取 tol = 0.5 m（人眼在 1:1000 比例下看不出差别），飞行品质分析取 tol = 0.1 m（要保留抖动特征），只做统计报表取 2~5 m。抽稀前先做一次 5~9 点滑动平均，否则噪声点会被 DP 当成「特征点」全留下来。

渲染侧用 Cesium 的 polyline 即可，注意两点：一是 clampToGround: false（要贴地就配 heightReference），二是点数超过 1 万时改用 PolylineCollection 或按时间分段加载，否则首屏卡顿。

3.5.3 实时拼接：特征 → 单应 → 位姿图

实时拼接的骨架是四步，每一帧都在跑：

帧 k-1, k
  │ ① ORB 特征提取（每帧 1500~2500 个）
  ▼
  │ ② 特征匹配 + ratio test（0.75）+ RANSAC 求单应 H
  ▼
  │ ③ 内点数 < 阈值 或 位移 > 阈值 → 标记为关键帧
  ▼
  │ ④ 关键帧之间做滑窗位姿图优化（g2o / Ceres），非关键帧直接沿用累积 H
  ▼
  │ ⑤ 用 POS 初值或地面控制点做地理配准，把像素坐标拉到经纬度
  ▼
  拼接图（带地理参考）

工程实现上，关键帧判定是三个条件的「或」：内点数 < 30、帧间位移 > 0.3 倍视野、旋转 > 5 度。只有关键帧才更新累积单应与位姿图，非关键帧直接复用上一帧结果，这是把 CPU 占用压下来的关键。

「实时」与「精度」的取舍体现在三个旋钮上：

实测参考：1080p 输入、关键帧间隔 0.3 FOV、滑窗 10 帧，在 8 核 CPU 上能跑 8~12 fps 的拼接刷新，端到端延迟 1.5~2.5 s；这个延迟对「实时看」够用，对「实时控制」不够。

3.5.4 离线高精度：ffmpeg 抽帧 → COLMAP → ODM → GDAL

要出可交付的正射影像，别在实时链路上加码，直接走离线 SfM 流水线：

# 1) 抽帧：视频 25 fps，抽成 1 fps（重叠率仍够），保留 EXIF/时间戳
mkdir -p frames
ffmpeg -i flight.mp4 -vf "fps=1" -q:v 2 frames/%05d.jpg

# 2) 若视频没有地理标签，用 POS 生成 COLMAP 的 priors（简化做法）
python tools/pos_to_colmap.py --pos pos.csv --frames frames --out colmap_priors.txt

# 3) 空三
colmap feature_extractor --database_path colmap.db --image_path frames \
    --ImageReader.single_camera 1 --ImageReader.camera_model OPENCV
colmap exhaustive_matcher --database_path colmap.db
mkdir -p sparse
colmap mapper --database_path colmap.db --image_path frames --output_path sparse
colmap model_analyzer --path sparse/0

# 4) 稠密重建 + 正射（OpenDroneMap 一条命令；也可 OpenMVS 手动跑）
docker run -ti --rm -v "$PWD/dataset:/datasets/project" opendronemap/odm \
    --project-path /datasets project --dsm --dtm --orthophoto-resolution 2 --pc-quality high

# 5) 转 COG 并核验
gdal_translate -of COG -co COMPRESS=DEFLATE -co BLOCKSIZE=512 \
    odm_orthophoto/odm_orthophoto.tif ortho_cog.tif
gdalinfo -json ortho_cog.tif | head -40

抽帧率的经验公式：帧间位移应小于单张覆盖宽度的 30%，即 fps >= v / (0.3 × W)。10 m/s、覆盖 150 m 时约 0.22 fps，取 1 fps 已经很冗余；反过来，视频 25 fps 全抽会得到 25 倍冗余影像，空三时间爆炸而精度几乎不涨。

3.5.5 质量评估：接缝、鬼影、几何精度

检查点布置比数量更重要：均匀覆盖测区四角与中心、避开建筑阴影与水面、与控制点至少相距 50 m。全部堆在测区一角，RMSE 数字再漂亮也说明不了问题。

坑 28：用控制点算 RMSE。控制点参与了空三解算，误差被「吸收」了，算出来必然偏小；必须留出独立的检查点。

坑 29：正射分辨率拍脑袋定。--orthophoto-resolution 2 表示 2 cm/px，比 GSD 还小并不会变清楚，只会让文件变大、处理变慢。

坑 30：忽略匀色。多架次拼接时不同时间的光照差异会形成明显色块，必须先做色差校正（ODM 的 --texturing-skip-global-seam-leveling 反过来说就是别关它）。

3.5.6 2026 路线：从「帧间单应拼接」升级到「位姿 + 双产出」

3.5.3 那套「ORB 特征 → 单应 → 滑窗位姿图」是 2020 年前后的主流做法：延迟低、实现短，但误差沿链路累积，只适合做「看起来连续」的预览。它的问题不在算法，而在目标定错了——把预览产物当交付产物用，接缝与漂移永远解释不清。2026 年的工程路线把它拆成两条通道：

视频 + POS
   │
   ├── 实时通道（预览）：ORB + 单应/稀疏位姿 + POS 约束
   │      → 拼接预览，端到端 1.5~2.5 s，米级，只看覆盖到哪、有无漏拍
   │
   └── 离线通道（交付）：抽帧 → 位姿（COLMAP SfM 或 VGGT 前馈）
          ├── 3DGS 训练（gsplat / brush）→ 带 LOD 的 splat tileset
          └── 稠密重建 + 正射（OpenDroneMap）→ DOM / DSM → COG

两条通道的公共中间层是「位姿 + 帧索引」，交接物用 COLMAP 文本格式最省事：实时侧不依赖它，离线侧所有工具都吃它，第06章已经把格式转换与踩坑走过一遍。把这一层定死之后，实时与离线不再互相绑架——预览卡了不影响交付，交付慢也不影响现场判断。

位姿来源有三条路，选型按这张表：

VGGT 的 2026 口径要写清楚：官方在 2026-05-15 修复了显存问题，同等显存可跑约 2~3 倍帧数；2026-05-18 发布了后继版本 VGGT-Omega。它对「先出个结果看看」特别有用，但三条限制没变：输出是相对尺度、长序列会漂移、没有全局 BA。商用项目还要注意授权——2025-07-29 起官方提供需申请的商用许可 checkpoint（VGGT-1B-Commercial），学术权重不能直接用于商用交付，条款以官方仓库为准。

离线双产出的流水线可以一条命令串起来（抽帧率、参数都以各工具当前 --help 为准）：

# 一次抽帧，两条产物：正射/DSM 走 ODM，splat 走 3DGS
mkdir -p frames work dataset/images
ffmpeg -i flight.mp4 -vf "fps=2" -q:v 2 frames/%05d.jpg

# ① 位姿：优先 COLMAP；注册率 < 80% 或时间不够时切 VGGT 前馈
colmap feature_extractor --database_path work/colmap.db --image_path frames \
    --ImageReader.single_camera 1 --ImageReader.camera_model OPENCV
colmap exhaustive_matcher --database_path work/colmap.db
colmap mapper --database_path work/colmap.db --image_path frames --output_path work/sparse

# ② 交付产物：正射 + DSM（ODM 要求 project/images 下有影像）
cp frames/*.jpg dataset/images/
docker run -ti --rm -v "$PWD/dataset:/datasets/project" opendronemap/odm \
    --project-path /datasets project --dsm --dtm --orthophoto-resolution 2 --pc-quality high

# ③ 展示产物：同一套 COLMAP 位姿喂给 3DGS 训练器（brush 示例）
brush work/dataset --total-train-iters 30000 --export-path work/splat/

两条通道的验收口径必须分开写，否则一定有人拿预览的延迟去比交付的精度：

注意：① 拿单应拼接的预览图当交付正射。接缝处的漂移在预览里看不出来，测图时才发现整景偏移，返工成本远高于直接走离线通道。② VGGT 输出直接当米制用。它的尺度来自训练分布，不做恢复就叠到 3D Tiles 或正射上，几何整体缩放，检查点一测就露馅。③ 双产出各训一遍位姿。正射与 splat 必须共用同一套位姿与同一份米制恢复参数，否则两个产品在同一场景里对不上，前端叠加时看起来像「模型漂移」。

3.6 系统集成：一条最小闭环

3.6.1 组件划分与接口契约

┌──────────────┐   RTMP/WHIP   ┌──────────────┐   gRPC     ┌──────────────┐
│ 设备接入服务  │ ────────────▶ │ 时空基准服务  │ ─────────▶ │ 视觉服务      │
│ 无人机/机场   │   MQTT POS    │ 坐标+时间转换 │  位姿+帧   │ 检测/位姿解算 │
└──────┬───────┘               └──────┬───────┘            └──────┬───────┘
       │ MQTT telemetry               │ 统一时空对象                │ 结构化结果
       ▼                              ▼                            ▼
┌──────────────────────────────────────────────────────────────────────┐
│ 消息总线（MQTT 测控 / Kafka 事件）→ 三维可视化（Cesium）→ 成果库      │
└──────────────────────────────────────────────────────────────────────┘

接口设计的一条硬规矩：每个跨服务消息都带 ts_utc_ms 与 crs 两个字段。没有时间戳的位姿、没有 CRS 的坐标，一律当脏数据处理。

3.6.2 时间同步：NTP 够用，PTP 才够稳

偏差检测不要只看「时钟同步了没有」，要看「同一条数据在两个通道里的时间差」的统计量：offset 取中位数抗异常值，MAD 判断抖动是否稳定。两者都稳定才说明固定偏移可补偿；MAD 随时间漂移说明时钟在跑偏，先解决同步再谈融合。

3.6.3 消息总线：MQTT 管测控，Kafka 管事件

数据契约用 JSON Schema 固定下来，前后端都按它校验：

{
  "schema_version": "1.2",
  "type": "object",
  "required": ["device_id", "ts_utc_ms", "crs", "seq", "pose"],
  "properties": {
    "device_id": {"type": "string"},
    "ts_utc_ms": {"type": "integer", "minimum": 1600000000000},
    "crs": {"type": "string", "enum": ["WGS84", "CGCS2000"]},
    "seq": {"type": "integer", "minimum": 0},
    "pose": {
      "type": "object",
      "required": ["lon", "lat", "alt_ellipsoid"],
      "properties": {
        "lon": {"type": "number", "minimum": -180, "maximum": 180},
        "lat": {"type": "number", "minimum": -90, "maximum": 90},
        "alt_ellipsoid": {"type": "number"},
        "alt_agl": {"type": "number"},
        "yaw_deg": {"type": "number"}
      }
    }
  }
}

3.6.4 联调顺序：先回放，再半实物，最后实飞

坑 31：跳过回放直接实飞。90% 的坐标和时间错误在离线回放阶段就能暴露，且不花钱。

坑 32：总线里塞视频原始帧。视频走流媒体通道（RTMP/WebRTC），总线只传元数据，否则单条消息几百 KB，MQTT broker 直接被打挂。

3.7 用 DeepSeek Harness 做联调与仿真

3.7.1 造数据：让 DSH 生成回放素材

真实的 POS、视频、DEM 不好凑，先用 DSH 造一套可复现的。把任务书写成「可验收」的形式：

# 在项目根目录开一个 headless 会话，任务书要带验收标准
dsh --profile headless "在 tools/ 下生成两个脚本：
1) gen_fake_pos.py：输出 10 Hz 的 POS CSV，字段 time_utc,lon,lat,alt_ellipsoid,roll,pitch,yaw；
   生成 120 s 的圆形航迹（半径 80 m），叠加 3 cm 高斯噪声与 0.8 s 时间偏移；
2) gen_fake_dem.py：用 numpy 生成 512x512 的 fBm 地形，高程 0~120 m，
   用 rasterio 写成 EPSG:4547 的 GeoTIFF，像元 1 m。
验收：两个脚本 python -m py_compile 通过；生成文件能被 pandas/rasterio 读回；
      航迹闭合误差 < 0.5 m；DEM 的 min/max 与打印值一致。"

自己写也就二十行：采样率、航迹半径、噪声水平是仅有的三个参数。把它固化成脚本而不是一次性代码，后面每次改算法都能重跑同一套数据。

3.7.2 写测试：坐标转换要用「已知点对」锁死

坐标转换的测试不能只测「跑通」，要测「数值对不对」。四类用例覆盖足够：

让 DSH 跑回归并把失败用例整理成最小复现：

dsh --profile headless "运行 pytest tests/test_crs.py -q 和 tests/test_time_align.py；
如果有失败，只输出：失败用例名、实际值 vs 期望值、最小复现代码片段（不超过 20 行）；
不要改任何源码，只给诊断。" --json > test_report.json

3.7.3 日志分析：先机器筛窗口，再让人/模型判根因

飞控日志动辄几十万行，直接丢给模型是浪费。先用阈值把「可疑时段」筛成几十行的窗口（视觉置信度、位姿跳变、IMU 方差），再把窗口交给 DSH 判根因，命中率高得多。

dsh --profile headless "vision_log 里 confidence 低于 0.5 且连续超过 1 s 的窗口有 3 段（附 CSV）。
同期飞控日志的 IMU 加速度方差上升 4 倍。判断最可能原因，给出：
1) 排除顺序；2) 现场复现步骤；3) 代码里应加的保护逻辑。不要改文件。"

3.7.4 沉淀成技能：一次调通，次次复用

把「坐标链路验收」这类流程写进 $DSH_HOME/skills/<name>/SKILL.md，下次直接加载，不用重新描述：

---
name: crs-check
description: 校验一条空地数据链的坐标与时间基准，输出偏差表
---

1. 读取 pos.csv 与 video_meta.json，检查字段与 CRS 声明
2. 跑 tests/test_crs.py、tests/test_time_align.py
3. 对 3 个实测控制点做往返比对，输出平面/高程偏差
4. 时间残差 > 20 ms 或坐标偏差 > 2 cm 时，给出排查顺序（不改代码）

放进 CI 就是一条命令的事：

# 只读模式跑审查任务，产出机器可读事件流
dsh --profile headless "按 skills/crs-check 的步骤跑一遍，输出偏差表 markdown" \
    --permission-mode read-only --json > crs_report.json

坑 33：让 DSH 直接改生产代码。先让它出诊断和补丁建议，人审之后再落地；权限预设用 read-only 或 workspace-write 分级。

坑 34：造数据不带随机种子。每次生成的 POS 不一致，测试就不可复现，回归等于没做。

3.8 本章小结

1. 空地一体的本质是统一基准：坐标、时间、度量三件事没做对，后面所有融合都是在还债。
2. 判断一条链路是否可信，先看四个数：坐标往返残差、时间残差 RMS、漂移率、帧—位姿配对率。
3. 高程基准是最高频的隐形事故：椭球高与 1985 国家高程基准相差 10~40 m，必须用格网模型换算。
4. 地图表示先算内存：1 km² 用 0.5~1 m 的 2.5D 栅格做全局，0.1~0.2 m 局部栅格做重规划，3D 占用网格只在必须穿障时上。
5. 硬约束与软代价必须分开：禁飞区、限高、膨胀半径进 blocked，坡度与视觉质量进 cost，A* 只判断有限性。
6. 航高、GSD、重叠率三者互锁：先定成图比例尺，再定 GSD、航高、航带间距与航时。
7. 视觉降落必须分层降级：8 m 以上靠 RTK，5 m 内切平面标记，1 m 内融合测距与触地判定。
8. 平面标记的位姿存在二义性，必须用上一帧连续性约束选解；边长实测误差会直接放大成深度误差。
9. AR 投影的精度上限由地面模型决定：平地可用透视四边形贴图，起伏地形必须按 DEM 网格求交。
10. 轨迹与拼接拆成两个服务：轨迹秒级刷新，拼接滑窗优化，混在一起必然抢资源。
11. 验收要有独立检查点：平面 RMSE ≤ 1 × GSD、高程 RMSE ≤ 2 × GSD，检查点与控制点不能共用。
12. 联调顺序是离线回放 → 半实物 → 系留 → 实飞，九成的坐标与时间错误在回放阶段就能暴露。
13. 底图精度按来源分级承诺：公共实景三维成果给的是 2~5 m 格网与年度增量，1:500 只能自采 + 控制点 + 独立检查点，方案里写检查点 RMSE 而不是标称精度。
14. 三维底座分两层：3D Tiles 1.1 承载可量测的测绘成果与 metadata 属性，带 LOD 的 splat tileset 承载近景真实感，两者必须在同一 CRS/ECEF 上叠加。
15. 拼接路线从帧间单应升级为「前馈/SfM 位姿 + 双产出」：实时通道只出预览，交付级正射与 splat 走离线通道，中间用 COLMAP 格式解耦。

3.9 动手练习

• 练习 1（坐标链路自检）：输入 3 个 RTK 实测点（WGS84）与对应的 CGCS2000 高斯坐标；用 pyproj 做双向转换、numpy 手写 ENU 正反算，并补 pytest 用例。验收：往返平面残差 < 1 mm、高程 < 2 mm，已知点对偏差平面 < 2 cm、高程 < 5 cm。
• 练习 2（GSD 反算航线）：输入 1 英寸 20 MP 相机（像元 2.4 μm、8.8 mm 定焦）与 2 km × 3 km 任务区；给出 1:1000 成图对应的 GSD、航高、航带间距、航带数与预估航时。验收：GSD ≤ 10 cm、旁向重叠 60%，航时估计与实飞记录偏差 < 15%。
• 练习 3（标记降落闭环）：输入 40 cm 标记板、1080p 下视相机与 PX4 SITL；实现 3.3.4 的五段状态机，含超时与失锁分级。验收：仿真中连续 20 次降落，水平偏差 < 15 cm，触地速度不超过 0.3 m/s。
• 练习 4（最小 AR 投影）：输入一段带 POS 的 30 s 视频与 1 m DEM；按 DEM 网格做反投影，输出每帧地面覆盖多边形与偏差统计。验收：平坦区偏差 < 2 m、起伏区 < 1 × GSD，俯仰角超 70 度时正确拒绝显示。

3.10 延伸阅读

• PROJ 官方文档
• pyproj 文档
• CesiumJS 版本发布记录（1.147）
• OpenCV calib3d 模块文档
• COLMAP 官方仓库
• OpenDroneMap 文档
• GDAL 官方文档
• 3D Tiles 规范与工具（CesiumGS/3d-tiles）