第08章 - 空间数据底座与时空基准

项目做到第二个月，总会出现一场所有人都绕不开的会：同一栋楼在两张图上差了 47 米；轨迹回放时车总是提前两秒拐弯；三个月前那份统计报告今天怎么都复算不出来。三个问题分别指向坐标基准、时间基准与数据版本——这就是本章的三条主线。

空间智能的上层全是算法：检测、分割、建图、规划、问答。可算法跑在数据上，底座歪一寸，上层要还的债是一丈。本章不讲模型，只讲怎么把柱子打到岩层里：基准怎么定、矢量与影像怎么建库、三维成果怎么组织、数据怎么治理，以及怎么把这些检查变成 DeepSeek Harness（下称 DSH）能自动执行的技能。凡涉及具体版本、日期的判断，都以本章末尾「延伸阅读」里的官方文档为准。

学习目标：能独立交付一套「基准—组织—治理」方案——写得出一条 WGS84 ↔ CGCS2000 高斯 ↔ 地方独立系 ↔ ENU 的往返代码并自检到厘米级；能判定 UTC、GPS 周秒、本地时间与视频 PTS 的偏差来源，把端到端残差压到 20 ms 以内；能建出带 SRID 约束、分区与索引的 PostGIS 3.7 库，并用 EXPLAIN 证明查询真的走了索引；能把影像组织成 COG + STAC 1.1 的可检索集合；能用 DuckDB 直查 GeoParquet 完成千万级要素统计；能给出三维底座的组织结构与显存预算；最后把数据体检做成 DSH 技能接进 CI。 前置知识：会写 Python 与 SQL；用过 GDAL 命令行；知道 EPSG 编码、GiST 索引、COG 这几个词；能读懂 EXPLAIN 的基本输出。

8.1 为什么先讲数据底座：三条血泪教训

空间智能项目很少败在模型精度上。查到最后，「算法不灵」的现场大都是数据的问题：两张图对不上、两路视频不同步、三个月前的结论今天复现不出来。这三类故障分别对应坐标基准、时间基准与数据版本。它们不上大屏、不进论文、不出 demo，却决定了上层一切结论能不能被信任。

下面三个案例都来自真实项目，数字做过脱敏，但量级与成因没有改。

8.1.1 教训一：坐标系混用，全市管网整体偏移 47 米

某市地下管网普查，作业队把一批 2018 年的「西安80」成果直接当成 CGCS2000 贴进一张图，省掉了七参数转换这一步。入库后图面一切正常：管线接得严丝合缝、阀门正好在管点上。直到现场核验，才发现所有要素相对实测控制点西向偏 31 m、北向偏 35 m，合成 46.7 m——管线横穿河道，阀门落在楼里。

这类错误最阴险的地方是「自洽」：整批数据一起偏，图面拓扑关系仍然正确，任何只在本数据集内部做的检查都发现不了。只有与外部基准（实测控制点、影像、其它部门成果）套合才会暴露。

常见的基准混用形态与量级：

8.1.2 教训二：时间不同步，轨迹整体错位 30 米

某园区巡检项目，POS 是 10 Hz、时间戳为 UTC；视频用 MP4 封装，PTS 从 0 开始。开发图省事，直接按「第 n 帧 = n/25 秒」把视频映射到 POS 首包。实际录制比 POS 首包晚 1.96 s，于是每一帧都取到了 1.96 s 之后的位置：60 km/h 的车速对应错位 32.7 m，15 m/s 的无人机对应 29.4 m。AR 标注框在画面上往回「爬」，越到后段越偏。

定位这个问题只需要一个动作：把两路时间戳按同一个物理事件（相机触发电平、画面里的 LED 闪光、场景里的闸机横杆）对齐，偏移量立刻现形。时间基座的要求不是「差不多」，是「有残差统计」。

常见时间源和它们各自的脾气：

8.1.3 教训三：版本不清，结论不可复现

同一片区域，2024Q3 的影像抽检出 1200 个图斑，2025Q1 的影像抽检出 1380 个。任务书只写了「使用最新影像」，报告里没写用了哪一份、哪天的、什么处理的。三个月后上级复核，谁也算不出 1200 这个数，整份成果只能作废重做。

版本问题的解法不是「约定一下」，而是让每一条数据自带身份证：采集时间、处理版本、内容哈希、上游来源，全部落进可查询的元数据表。复现时按 manifest 拉回同一批哈希，结论才可复算。这一节的具体做法放在 8.8。

8.1.4 「基准—组织—治理」三件事

三个案例对应三层工作，顺序不能颠倒：

              ┌──────────────────────────────────────────┐
   第 3 层    │ 治理：版本 / 血缘 / 权限 / 合规            │  可审计、可复现
              └────────────────────▲─────────────────────┘
              ┌────────────────────┴─────────────────────┐
   第 2 层    │ 组织：矢量库 / 影像集合 / 分析表 / 三维瓦片 │  可检索、可扩展
              └────────────────────▲─────────────────────┘
              ┌────────────────────┴─────────────────────┐
   第 1 层    │ 基准：坐标基准 / 时间基准 / 高程基准        │  可对齐、可量化
              └──────────────────────────────────────────┘
                       先钉死基准，再谈组织，最后谈治理

8.1.5 底座的五个硬指标

底座做得好不好，不能靠「感觉对了」，要能量化。下面五个指标建议写进项目验收文档，每条都配一条可执行的命令：

体检的最小命令集如下，每条都能给出一个数而不是一句「看着没问题」：

# 1) 矢量：要素数、几何类型、范围、CRS 一次看全
ogrinfo -so -al data/pipe.gpkg | grep -E "Layer name|Geometry:|Feature Count|Extent"

# 2) 影像：尺寸、波段数、CRS、是否 COG
gdalinfo -json data/ortho_2025q1.tif | jq '{size: .size, bands: (.bands|length), crs: .coordinateSystem.wkt[0:40]}'

# 3) 数据库：SRID 缺失数、无效几何数、空几何数
psql -d gis -tAc "SELECT count(*) FILTER (WHERE ST_SRID(geom) = 0) AS no_srid,
                         count(*) FILTER (WHERE NOT ST_IsValid(geom)) AS invalid,
                         count(*) FILTER (WHERE ST_IsEmpty(geom)) AS empty
                  FROM parcel;"

坑 1：先建库、后补基准。表建好了、数据进了、索引建了，才发现坐标系是错的，此时改一个字段要重导几十亿条记录。基准永远是第一步。

坑 2：迷信「万能中间坐标系」。每多跳一次转换就多叠一次误差，长链路还会把椭球高和正常高搅在一起。能一次转到位就不要分三段转。

坑 3：把「坐标看起来差不多」当验收。人眼在 1:500 图上分辨不出 5 m，但管线会穿墙、图斑会越界。验收必须写数值阈值，不能写「基本吻合」。

8.2 坐标基准：先把「在哪」钉死

坐标基准不是一个数，是「参考椭球 + 原点与定向 + 投影方式 + 单位与轴序」的组合。同样一个点，换一套组合就得到一组新数值，而所有数值都能自称「经纬度」。判断一份数据能不能用，先看它声明的 CRS 完整到什么程度：只有「WGS84」四个字的，等于没说；写到 EPSG 编码、带号、是否含带号、高程基准的，才算可用的声明。

8.2.1 三套基准与 EPSG 编码

国内项目绕不开三套基准，它们不是「精度不同」，而是「定义不同」：

CGCS2000 的投影 EPSG 编码有四个连续段，记住它们的对应关系，就能在拿到一个编码时立刻判断是几度带、带不带带号：

CGCS2000 与 WGS84 在 ITRF 框架下的差异在厘米级，工程上按无转换处理通常够用；但报批成果必须使用官方 CGCS2000 定义与官方转换参数，不要拿「差不多」去顶。

8.2.2 高斯-克吕格：中央子午线、带号与东坐标

高斯-克吕格是横轴墨卡托在椭球上的实现：中央子午线投影后是一条直线且长度不变形，离中央子午线越远变形越大，长度变形近似为

Δs / s ≈ (Δλ² · cos²φ) / 2

按这个式子代进去：3 度带边缘（Δλ = 1.5°）在纬度 30 度附近最大长度变形约 1/4000，也就是 0.25 m/km；6 度带边缘（Δλ = 3°）约 1/1000，即 1 m/km。城市大比例尺测图因此用 3 度带，或者干脆建地方独立系，用抵偿高程面把变形压下去。

带号和中央子午线的关系可以直接算，不用查表：

def gk_zone(lon_deg: float, width_deg: int = 3):
    """返回 (带号, 中央子午线经度)。lon_deg 为东经正、西经负。"""
    if width_deg == 3:
        n = round(lon_deg / 3.0)          # 深圳 114.31E → 38 带，中央子午线 114E
        cm = 3.0 * n
    elif width_deg == 6:
        n = int(lon_deg // 6) + 1         # 114.31E → 20 带，中央子午线 117E
        cm = 6.0 * n - 3.0
    else:
        raise ValueError("width_deg 只支持 3 或 6")
    return n, cm

for lon in (114.31, 116.40, 121.47):
    print(lon, "3 度带:", gk_zone(lon, 3), " 6 度带:", gk_zone(lon, 6))

同一个中央子午线，带号变体与不带带号变体的区别只有东坐标前面那两三位：

这两种变体在 EPSG 里是两组编码，混用一次就是三万八千公里的偏移——数据直接落到地球外面去。看到 easting 大于 1e6 时，先别急着画图，先确认它是不是带了带号。

8.2.3 ECEF、ENU 与地方独立坐标

• ECEF（Earth-Centered Earth-Fixed，EPSG:4978）：原点在地心，Z 轴指向协议北极，X 轴指向本初子午线与赤道交点。卫星轨道、惯导解算、三维引擎内部都用它，因为它是全局无奇点的直角坐标。
• ENU（East-North-Up）：以某个原点切出的局部平面，东-北-天三轴。飞控的本地坐标、点云局部拼接、机器人导航都用它。它是近似：切平面与地球表面在距离 d 处的偏差约为 d²/(2R)，10 km 处约 8 m，所以 ENU 只在十几公里内可信。
• 地方独立坐标系：测区自定义的中央子午线（常为过测区中心的子午线）加上抵偿高程面，再与 CGCS2000 之间用四参数（平面）或七参数（三维）联系。它的好处是长度变形接近零，坏处是每换一个测区就要换一套参数。

8.2.4 一条往返链的完整代码

工程实现的关键不在公式，在三个开关：轴序（always_xy）、带号（easting 是否含带号前缀）、高程基准（椭球高还是正常高）。下面这段代码把 WGS84 → CGCS2000 高斯 → 地方独立系 → ENU → 回到 WGS84 全部串起来，依赖只有 pyproj 与 numpy：

# crs_roundtrip.py  |  pip install "pyproj>=3.6" "numpy>=1.26"
import numpy as np
from pyproj import CRS, Geod, Transformer

CRS_WGS84 = CRS.from_epsg(4326)     # WGS84 地理
CRS_GK    = CRS.from_epsg(4547)     # CGCS2000 / 3 度带 CM 114E，不带带号
CRS_ECEF  = CRS.from_epsg(4978)     # ECEF

# 地方独立系四参数（示例值；真实项目必须用测绘部门下发的参数）
P = {"dx": -12584.319, "dy": 33871.207, "theta": np.deg2rad(0.00137), "k": 1.0000217}
ORG_E, ORG_N = 500000.0, 2490000.0  # 地方系原点在 114E 带内的平面位置

t_wgs_gk   = Transformer.from_crs(CRS_WGS84, CRS_GK, always_xy=True)
t_gk_wgs   = Transformer.from_crs(CRS_GK, CRS_WGS84, always_xy=True)
t_wgs_ecef = Transformer.from_crs(CRS_WGS84, CRS_ECEF, always_xy=True)
geod = Geod(ellps="WGS84")

def gk_to_local(e, n):
    c, s = np.cos(P["theta"]), np.sin(P["theta"])
    de, dn = e - ORG_E, n - ORG_N
    return (P["dx"] + P["k"] * (c * de - s * dn),
            P["dy"] + P["k"] * (s * de + c * dn))

def local_to_gk(x, y):
    c, s = np.cos(P["theta"]), np.sin(P["theta"])
    x0, y0 = (x - P["dx"]) / P["k"], (y - P["dy"]) / P["k"]
    return c * x0 + s * y0 + ORG_E, -s * x0 + c * y0 + ORG_N

def ecef_to_enu(xyz, org_xyz, org_ll):
    lon, lat = np.deg2rad(org_ll[0]), np.deg2rad(org_ll[1])
    R = np.array([[-np.sin(lon),                np.cos(lon),                0.0],
                  [-np.sin(lat) * np.cos(lon), -np.sin(lat) * np.sin(lon), np.cos(lat)],
                  [ np.cos(lat) * np.cos(lon),  np.cos(lat) * np.sin(lon), np.sin(lat)]])
    return R @ (np.asarray(xyz, float) - np.asarray(org_xyz, float))

if __name__ == "__main__":
    org = (114.305392, 22.546117, 41.35)        # 测区原点
    pt  = (114.308177, 22.547290, 38.62)        # 待转换点

    e, n = t_wgs_gk.transform(pt[0], pt[1])     # ① WGS84 → 高斯平面
    x, y = gk_to_local(e, n)                    # ② 高斯平面 → 地方系
    enu  = ecef_to_enu(t_wgs_ecef.transform(*pt),
                       t_wgs_ecef.transform(*org), org[:2])   # ③ → 本地 ENU

    e2, n2 = local_to_gk(x, y)                  # ④ 地方系 → 高斯平面
    lon2, lat2 = t_gk_wgs.transform(e2, n2)     # ⑤ 高斯平面 → WGS84
    _, _, res = geod.inv(pt[0], pt[1], lon2, lat2)

    print(f"WGS84 : {pt[0]:.6f}, {pt[1]:.6f}, {pt[2]:.2f}")
    print(f"高斯  : E={e:.3f} N={n:.3f}")
    print(f"地方系: X={x:.3f} Y={y:.3f}")
    print(f"ENU   : E={enu[0]:.3f} N={enu[1]:.3f} U={enu[2]:.3f}")
    print(f"回程  : {lon2:.9f}, {lat2:.9f}")
    print(f"平面往返残差 {res * 1000:.4f} mm（应 < 0.1 mm）")

注意高程没有参与这条链：EPSG:4547 是二维投影，高度原样穿过。椭球高到 1985 国家高程基准正常高的换算要单独做（用高程异常或似大地水准面模型），投影库不会顺手帮你换——这正是 8.1.1 里那个「高程差 40 米」的来源。

8.2.5 自检：三条不变量加一组控制点

没有外部控制点时，坐标链路也能自检，因为投影定义本身蕴含了可验证的数学性质。三条不变量不需要任何外部数据，跑通即说明实现与 EPSG 定义一致：

# crs_selfcheck.py —— 不需要外部控制点，30 秒跑完
import math
from pyproj import CRS, Transformer

A, F = 6378137.0, 1 / 298.257223563
E2 = F * (2 - F)                 # 第一偏心率平方
CM, FE = 114.0, 500000.0         # EPSG:4547 的中央子午线与东偏常数

t_fwd = Transformer.from_crs(CRS.from_epsg(4490), CRS.from_epsg(4547), always_xy=True)
t_inv = Transformer.from_crs(CRS.from_epsg(4547), CRS.from_epsg(4490), always_xy=True)

# 不变量 1：中央子午线上 easting 恒等于东偏常数
for lat in (18.0, 22.5, 30.0, 45.0):
    e, _ = t_fwd.transform(CM, lat)
    assert abs(e - FE) < 1e-6, (lat, e)

# 不变量 2：往返一致，1e-9 度约等于 0.1 mm
for lon, lat in ((114.308177, 22.547290), (113.9, 23.1), (114.9, 21.8)):
    lon2, lat2 = t_inv.transform(*t_fwd.transform(lon, lat))
    assert abs(lon2 - lon) < 1e-9 and abs(lat2 - lat) < 1e-9, (lon, lat, lon2, lat2)

# 不变量 3：离 CM 一毫度处的 easting 增量等于 N·cosφ·Δλ（N 为卯酉圈曲率半径）
D = math.radians(0.001)
for lat in (18.0, 22.5, 30.0):
    s = math.sin(math.radians(lat))
    expect = A / math.sqrt(1 - E2 * s * s) * math.cos(math.radians(lat)) * D
    e1, _ = t_fwd.transform(CM, lat)
    e2, _ = t_fwd.transform(CM + 0.001, lat)
    assert abs((e2 - e1) - expect) < 1e-3, (lat, e2 - e1, expect)

print("三条不变量全部通过：EPSG:4547 的实参与官方定义一致")

再由项目实测控制点兜底，把「实现对不对」升级为「参数对不对」：

# control_points.csv 表头：name,lon,lat,e_expected,n_expected（平面 m，控制点由测绘部门提供）
import csv
TOL_E, TOL_N = 0.010, 0.010      # 平面 1 cm；高程单独用 0.03 m 判
with open("control_points.csv", newline="", encoding="utf-8") as fh:
    for row in csv.DictReader(fh):
        e, n = t_fwd.transform(float(row["lon"]), float(row["lat"]))
        assert abs(e - float(row["e_expected"])) <= TOL_E, row["name"]
        assert abs(n - float(row["n_expected"])) <= TOL_N, row["name"]
print("控制点比对通过")

四层校验按成本从低到高排列，每次都从第一层开始跑：

坑 4：EPSG:4326 的官方轴序是（纬度, 经度）。用 pyproj 时忘写 always_xy=True，传入 (lon, lat) 会被当成 (lat, lon)，输出一个「看着还挺合理」的坐标，直到套合时才发现南北互换。

坑 5：把带带号变体（EPSG:4513–4533）当成不带带号变体（EPSG:4534–4554）。easting 多出三万八千米量级，投影后数据落在地球外面，还以为是渲染问题。

坑 6：拿 ENU 当区域坐标用。ENU 是切平面近似，10 km 外与地球表面的偏差约 8 m；跨区拼接请回到高斯平面或 ECEF。

8.3 时间基准：把「何时」对齐到毫秒

一次采集的时间戳至少经过五跳：传感器时钟 → 采集板 → 编码器 → 封装容器 → 存储介质。每一跳都可能引入固定偏移或随机抖动，你在终端看到的那串时间，只是最后一跳的结果。做时间基准，本质是把这五跳的偏移逐一标定，并留下可复算的残差统计。

8.3.1 四种时间源与它们的偏移

一行算术就能说明为什么不能只靠软件打点：1 ppm 的晶振漂移等于 86.4 ms/天，20 ms 的容差最多撑 5 小时。任何跨设备、跨天的对齐，都必须有外部授时源参与。

8.3.2 授时：NTP、PTP 与 GPS 授时的选择

8.3.3 偏差检测脚本：把「同步了没」变成两个数

判据先定死：端到端残差 RMS < 20 ms、P95 < 50 ms。20 ms 不是拍脑袋，它是「空间误差 30 cm 以内」的等价条件：

检测方法是拟合一条直线，而不是逐帧做差。因为两路时钟之间存在固定偏移 b 和相对漂移 (a-1)，逐帧相减会把漂移混进随机误差里：

# time_check.py  —— 判据：RMS < 20 ms 且 P95 < 50 ms
import numpy as np
import pandas as pd

TOL_RMS, TOL_P95 = 0.020, 0.050          # 秒

def check(path, col_a="trigger_utc", col_b="pts_utc"):
    df = pd.read_csv(path, parse_dates=[col_a, col_b]).dropna()
    x = df[col_a].astype("int64").to_numpy() / 1e9   # 参考时间（日志/授时打标）
    y = df[col_b].astype("int64").to_numpy() / 1e9   # 被测时间（容器 PTS 换算的 UTC）
    A = np.vstack([x - x[0], np.ones_like(x)]).T
    (a, b), *_ = np.linalg.lstsq(A, y - y[0], rcond=None)
    res = (y - y[0]) - (a * (x - x[0]) + b)          # 去趋势后的残差
    rms = float(np.sqrt(np.mean(res ** 2)))
    p95 = float(np.percentile(np.abs(res), 95))
    ok = rms < TOL_RMS and p95 < TOL_P95
    print(f"{path}: n={len(df)} offset={b * 1000:+.1f} ms "
          f"drift={(a - 1) * 1e6:+.1f} ppm rms={rms * 1000:.2f} ms "
          f"p95={p95 * 1000:.2f} ms -> {'PASS' if ok else 'FAIL'}")
    return ok

if __name__ == "__main__":
    assert all(check(p) for p in ("logs/cam01.csv", "logs/lidar_front.csv")), "时间基准未达标"

输出里两个数各有用途：offset 是固定偏移，说明起始时刻没对齐，改时间戳即可；drift 是相对漂移，单位 ppm，说明两路时钟不同源，必须引入授时，改代码没用。残差出现「台阶」则说明中途换了时钟域（闰秒、夏令时、设备重启），台阶位置就是排查入口。

8.3.4 存储约定：只存 UTC，只存有类型的时刻

数据库里时间字段的规矩只有一条：全部存 UTC，全部用带时区类型。PostgreSQL 的 timestamptz 存的是 UTC 瞬时值、显示时才按时区转换；timestamp 什么都不做，你把本地时间写进去，它原样存下来，跨时区查询时静默错 8 小时。

CREATE TABLE base.frame (
  id          bigserial PRIMARY KEY,
  stream_id   text        NOT NULL,
  pts         bigint      NOT NULL,            -- 容器时间基计数
  time_base   text        NOT NULL,            -- 例如 1/90000
  t_start_utc timestamptz NOT NULL,            -- 该流首帧对应的绝对 UTC
  acquired_at timestamptz NOT NULL,
  ingested_at timestamptz NOT NULL DEFAULT now(),
  CONSTRAINT frame_pts_nonneg CHECK (pts >= 0)
);
COMMENT ON COLUMN base.frame.pts IS '绝对时间 = t_start_utc + pts * time_base';

从视频里把 PTS 和时间基取出来，一条命令就够：

ffprobe -v error -select_streams v:0 -show_entries stream=time_base,start_time -show_entries frame=pts_time -of csv=p=0 clip.mp4 | head -5

坑 7：把 timestamp（无时区）当 UTC 存。开发机在东八区、服务器在 UTC，同一个值在两台机器上代表两个时刻，测试全过、上线全错。

坑 8：把视频 PTS 直接当 UTC 用。PTS 是相对量，原点是录制起点或首帧；MP4 的 edit list 还能让 start_time 为正、为负或不为零。

坑 9：用「帧号 ÷ 帧率」推算时刻。安防摄像头普遍是可变帧率（VFR），帧号与时间不是线性关系，必须读每帧 PTS。

8.4 矢量底座：PostGIS 3.7 建库、索引与查询

矢量底座要同时满足三个约束：存得进（CRS 明确、几何有效）、查得快（索引与分区匹配查询形态）、查得准（同一张表里只有一个 SRID）。PostGIS 把能力都给了，但规矩得自己立——它不会阻止你在一个 geometry 列里混放三种坐标系。

8.4.1 版本、依赖与自检

本章以 PostGIS 3.7.0（2026-10-29 发布）为准；需要更保守的维护版本时用 3.6.4（2026-06 稳定维护版）。3.7 系列的要求是 GEOS 3.10+、PostgreSQL 14–19、PROJ 6.1+，装之前先确认这三个组件都在射程内。

-- 扩展是否可装、装的是哪个版本
SELECT name, default_version, installed_version
FROM pg_available_extensions
WHERE name IN ('postgis', 'postgis_raster', 'postgis_topology');

CREATE EXTENSION IF NOT EXISTS postgis;
CREATE EXTENSION IF NOT EXISTS postgis_raster;

-- 一次看全：PostGIS / GEOS / PROJ 各自版本
SELECT postgis_full_version();

8.4.2 geometry 还是 geography

选错类型不会报错，只会让你的查询慢十倍或者算错单位：

判据很直接：测区在一个 3 度带内、面积不超过一个省，用 geometry + CGCS2000 投影坐标（如 EPSG:4547）；数据跨带或跨半球，用 geography，或者存 4326 并按需投影到局部坐标系再算。两者都想要的做法是「主表存 geometry，另建 geography 生成列用于距离查询」，代价是多一倍存储。

8.4.3 建库 DDL：把规矩写进类型里

下面这份 DDL 把三条规矩做成了数据库层面的硬约束：SRID 写进类型修饰符（写错 SRID 直接失败）、几何有效性写成 CHECK、按时间分区（让索引重建和 VACUUM 可控）：

CREATE SCHEMA IF NOT EXISTS base;

CREATE TABLE base.parcel (
  id          bigserial,
  code        text        NOT NULL,
  geom        geometry(Polygon, 4490) NOT NULL,   -- SRID 进类型，别的 SRID 写不进来
  quality     smallint    NOT NULL DEFAULT 1,
  acquired_at timestamptz NOT NULL,
  ingested_at timestamptz NOT NULL DEFAULT now(),
  props       jsonb       NOT NULL DEFAULT '{}'::jsonb,
  PRIMARY KEY (id, acquired_at),                  -- 分区表主键必须含分区键
  CONSTRAINT parcel_geom_valid    CHECK (ST_IsValid(geom)),
  CONSTRAINT parcel_geom_not_empty CHECK (NOT ST_IsEmpty(geom)),
  CONSTRAINT parcel_quality_range CHECK (quality BETWEEN 0 AND 5)
) PARTITION BY RANGE (acquired_at);

CREATE TABLE base.parcel_2026q1 PARTITION OF base.parcel
  FOR VALUES FROM ('2026-01-01') TO ('2026-04-01');
CREATE TABLE base.parcel_2026q2 PARTITION OF base.parcel
  FOR VALUES FROM ('2026-04-01') TO ('2026-07-01');

-- 主力：几何检索走 GiST
CREATE INDEX parcel_geom_gix  ON base.parcel USING GIST (geom);
-- 时间顺序写入的大表用 BRIN，体积只有 B-tree 的百分之一量级
CREATE INDEX parcel_time_brin ON base.parcel USING BRIN (acquired_at) WITH (pages_per_range = 64);
-- 属性等值查询走普通 B-tree
CREATE INDEX parcel_code_btree ON base.parcel (code);

分区粒度按写入量倒推：单个分区控制在 2000 万行以内，季度分区在千万级要素的区县项目里刚好；如果一个月就写 5000 万行，改按月分区。分区键必须是查询里最常出现的时间范围条件，否则分区裁剪（partition pruning）不会生效，等于白分。

8.4.4 三类查询与 EXPLAIN 判读

张量底座的查询形态只有三类：范围检索、最近邻、聚合统计。三类各配一种索引，判读 EXPLAIN 时只看三个数：节点类型（有没有 Index Scan）、Rows Removed by Filter（索引筛得干不干净）、Buffers（读了多少页）。

第一类，范围检索。关键是先用 && 这个纯几何包围盒运算走 GiST，再用精确谓词过滤：

-- 检索某季度内落在给定矩形中的地块，并给出平方米面积
EXPLAIN (ANALYZE, BUFFERS)
WITH win AS (
  SELECT ST_MakeEnvelope(113.90, 22.50, 113.95, 22.55, 4490) AS g
)
SELECT p.id, p.code, round(ST_Area(p.geom::geography)::numeric, 2) AS area_m2
FROM base.parcel p, win
WHERE p.acquired_at >= '2026-01-01' AND p.acquired_at < '2026-04-01'   -- 分区键
  AND p.geom && win.g                                                  -- 走 GiST
  AND ST_Intersects(p.geom, win.g)                                     -- 精确过滤
ORDER BY p.id
LIMIT 100;

Limit  (cost=0.55..412.37 rows=100 width=40) (actual time=0.041..2.918 rows=100 loops=1)
  Buffers: shared hit=186
  ->  Nested Loop  (cost=0.55..41822.11 rows=10124 width=40) (actual time=0.039..2.881 rows=100 loops=1)
        ->  Index Scan using parcel_geom_gix on parcel_2026q1 p
              Index Cond: (geom && '0103...'::geometry)
              Filter: st_intersects(geom, '0103...'::geometry)
              Rows Removed by Filter: 231
Execution Time: 2.9 ms

判读要点：出现 Index Scan using parcel_geom_gix 说明 GiST 生效；Rows Removed by Filter 只有 231，说明包围盒筛选已经很接近精确结果，索引质量好；如果这个数涨到几万，说明数据里混了很多「大而空」的面，先用 ST_Subdivide 切小再入库。实际数字随数据量浮动，看量级与节点类型，不要照抄数字。

第二类，最近邻。KNN 靠 GiST 的 <-> 距离算子排序，LIMIT 会触发索引的「边扫边排序」：

-- 500 m 内的最近 10 个地块，单位米，用 geography 做球面距离
EXPLAIN (ANALYZE, BUFFERS)
SELECT p.id, p.code,
       round(ST_Distance(p.geom::geography, pt.g::geography)::numeric, 1) AS dist_m
FROM base.parcel p,
     (SELECT ST_SetSRID(ST_MakePoint(114.305392, 22.546117), 4490) AS g) pt
WHERE ST_DWithin(p.geom::geography, pt.g::geography, 500)
ORDER BY p.geom::geography <-> pt.g::geography
LIMIT 10;

第三类，聚合统计。分组聚合的瓶颈通常在 IO 而不是索引，用分区裁剪把范围压到 1–2 个分区最有效：

-- 按质量等级统计图斑数与面积（只扫 2026Q1 分区）
SELECT quality,
       count(*)                                              AS n,
       round((sum(ST_Area(geom::geography)) / 10000)::numeric, 2) AS area_ha
FROM base.parcel
WHERE acquired_at >= '2026-01-01' AND acquired_at < '2026-04-01'
GROUP BY quality
ORDER BY quality;

8.4.5 索引维护与批量导入

批量导入的顺序错了，后面所有查询都要还债。推荐流程：建无索引的临时表 → COPY 灌数 → 校验几何 → 一次建索引 → ANALYZE → 按几何列 CLUSTER 一次。

-- 1) 灌数（比逐条 INSERT 快一个数量级以上）
\copy base.parcel (code, geom, quality, acquired_at) FROM 'parcel_2026q1.csv' WITH (FORMAT csv, HEADER true)

-- 2) 校验：无效几何、空几何、SRID 错误各有多少
SELECT count(*) FILTER (WHERE NOT ST_IsValid(geom))       AS invalid,
       count(*) FILTER (WHERE ST_IsEmpty(geom))           AS empty,
       count(*) FILTER (WHERE ST_SRID(geom) <> 4490)      AS wrong_srid
FROM base.parcel WHERE acquired_at >= '2026-01-01' AND acquired_at < '2026-04-01';

-- 3) 修掉自相交（只对无效行做，别全表跑）
UPDATE base.parcel SET geom = ST_MakeValid(geom)
WHERE NOT ST_IsValid(geom)
  AND acquired_at >= '2026-01-01' AND acquired_at < '2026-04-01';

-- 4) 建索引 → 统计信息 → 物理聚簇
CREATE INDEX IF NOT EXISTS parcel_geom_gix ON base.parcel USING GIST (geom);
ANALYZE base.parcel;
CLUSTER base.parcel USING parcel_geom_gix;   -- 大表耗时，放在业务低峰

坑 10：在 EPSG:4490 的 geometry 上直接调 ST_Area，得到的是平方度而不是平方米。要么转成 geography（ST_Area(geom::geography)），要么先投影到 4547 再算。距离同理，ST_Distance 在 4490 上返回度。

坑 11：在 WHERE 里给列套函数，例如 WHERE ST_Transform(geom, 4547) && …。索引是对列本身建的，列被函数包住就用不上，执行计划退化为全表扫。要用投影坐标就建表达式索引，或直接存两列。

坑 12：批量导入后不 ANALYZE。规划器按默认选择率瞎估，rows 估成一万实际一千万，于是选了 Nested Loop，查询从毫秒变成分钟。

8.5 影像与栅格：COG、STAC 与对象存储

影像底座的矛盾很集中：单景动辄几 GB 到几十 GB，用户却只想看一个 256×256 的瓦片。解法是把「一个文件」变成「一个可按字节区间读取的容器」，再把「一堆文件」变成「一个可检索的目录」。

8.5.1 COG：为什么 range 请求能救你的带宽

COG（Cloud Optimized GeoTIFF）不是新格式，是对 GeoTIFF 的布局约束：文件头与索引放在最前面，金字塔（overview）按「由粗到细」排列，瓦片数据分块存储并带内部索引。客户端先取文件头，读到各层级的瓦片偏移表，然后只请求需要的那几百 KB。

COG 文件内部布局
┌────────────┬─────────────┬───────────────────┬────────────────────┐
│ TIFF 头    │ 主 IFD      │ 概览 IFD 1..n     │ 瓦片数据（分块）    │
│ (前 16 KB) │ 指向各层    │ 由粗到细           │ 每块 256/512 像素   │
└────────────┴─────────────┴───────────────────┴────────────────────┘
      ▲              ▲                                    ▲
      │ 第 1 次 range│ 第 2 次 range（读瓦片偏移表）        │ 第 3 次 range
      └──────────────┴────────────────────────────────────┘
              三次请求即可显示一屏，与整景大小无关

量级对照（以一景 20 cm 分辨率、覆盖一个区县的正射影像为例，整景压缩后约 8 GB）：

生成与校验 COG 用 rio-cogeo，两个命令就够：

# 生成：内部 512 分块、DEFLATE 压缩、自动建概览、保留 6 级金字塔
rio cogeo create ortho_2025q1.tif ortho_2025q1_cog.tif \
  --cog-compress DEFLATE --cog-blocksize 512 --overview-level 6 \
  --overview-resampling average

# 校验：不通过会列出具体原因（非瓦片、概览缺失、IFD 顺序错等）
rio cogeo validate ortho_2025q1_cog.tif

校验通过只说明布局合规，不说明服务端支持 range。上线前做一次「半开区间」测试：用 curl 只取前 16 KB 和末尾 16 KB，两者都应该在百毫秒内返回，且 Content-Range 头正确。

curl -s -r 0-16383 -o /dev/null -w "%{http_code} %{size_download} %{time_total}\n" \
  https://oss.example.com/ortho/2025q1/ortho_2025q1_cog.tif

8.5.2 STAC：把「一堆文件」变成可检索目录

STAC（SpatioTemporal Asset Catalog）1.1.0 是目前稳定版本，它把影像目录抽象成三层：Catalog（目录）、Collection（集合）、Item（单景资产）。一个行政区一年的影像就是一个 Collection，每景是一个 Item，栅格文件挂在 Item 的 assets 里。

最小可用的 Item 长这样，字段一个都不能少：

{
  "type": "Feature",
  "stac_version": "1.1.0",
  "id": "GD-SZ-2025Q1-ortho-0p2m",
  "geometry": {"type": "Polygon", "coordinates": [[[113.75, 22.40], [114.62, 22.40],
                [114.62, 22.86], [113.75, 22.86], [113.75, 22.40]]]},
  "bbox": [113.75, 22.40, 114.62, 22.86],
  "properties": {
    "datetime": "2025-02-18T02:35:00Z",
    "gsd": 0.2,
    "proj:epsg": 4547,
    "eo:cloud_cover": 1.8,
    "platform": "aerial-photo"
  },
  "assets": {
    "data": {"href": "ortho_2025q1_cog.tif", "type": "image/tiff; application=geotiff; profile=cloud-optimized",
             "roles": ["data"]},
    "thumbnail": {"href": "thumb.jpg", "type": "image/jpeg", "roles": ["thumbnail"]},
    "metadata": {"href": "metadata.xml", "type": "application/xml", "roles": ["metadata"]}
  },
  "links": [
    {"rel": "root", "href": "../../catalog.json"},
    {"rel": "collection", "href": "./collection.json"},
    {"rel": "self", "href": "./GD-SZ-2025Q1-ortho-0p2m.json"}
  ]
}

STAC 1.1 的几条实践要点，直接用就能少踩坑：id 一旦发布就不要改（下游用它做外键）；datetime 必填，只有区间时用 start_datetime 加 end_datetime；links 里的 href 用相对路径，部署换域名不用改数据；proj:epsg 与 gsd 写在 properties 里，客户端才能在检索阶段就过滤掉分辨率不合格的景。

8.5.3 对象存储前缀设计：先想清楚怎么删

对象存储没有目录，前缀就是一切。前缀设计错了，最直接的代价是「想删一批数据时删不掉、想统计一批数据时数不清」。推荐按「数据类型 / 区域 / 时间 / 集合」四级分：

oss://gis-data/
└── ortho/                          # 数据类型：影像正射
    └── GD-SZ/                      # 区域：广东-深圳
        └── 2025Q1/                 # 时间：采集季度（用于生命周期与批量删除）
            └── ortho-0p2m/         # 集合：分辨率规格
                ├── collection.json
                ├── GD-SZ-2025Q1-ortho-0p2m.json     # Item
                ├── GD-SZ-2025Q1-ortho-0p2m.tif      # COG
                └── thumb.jpg

三条实践约定。一，时间用「年 + 季度」而不是完整日期，前缀数量可控、扫描成本低。二，不可变数据用版本前缀（v1/ v2/），需要回滚时换前缀而不是覆盖文件。三，清理规则挂在第 3 级：临时成果 90 天过期，正式成果永不自动删除，靠第 1 级的权限区分。

8.5.4 一个行政区影像集合的最小生成代码

围绕「一个行政区、一年、一档分辨率」这一组约束，pystac 三十行就能把 Collection 与 Item 写全：

# build_stac.py  |  pip install "pystac>=1.10" "shapely>=2.0"
import pystac
from shapely.geometry import box, mapping

COLL_ID = "gd-sz-ortho-0p2m"
coll = pystac.Collection(
    id=COLL_ID,
    description="深圳 0.2 m 正射影像，季度更新",
    extent=pystac.Extent(
        spatial=pystac.SpatialExtent([[113.75, 22.40, 114.62, 22.86]]),
        temporal=pystac.TemporalExtent([[None, None]]),
    ),
    license="proprietary",
)
coll.add_link(pystac.Link.root("../../catalog.json"))

def add_item(coll, item_id, dt, href, bbox):
    item = pystac.Item(id=item_id, geometry=mapping(box(*bbox)), bbox=bbox,
                       datetime=dt, properties={"gsd": 0.2, "proj:epsg": 4547})
    item.add_asset("data", pystac.Asset(
        href=href, media_type="image/tiff; application=geotiff; profile=cloud-optimized",
        roles=["data"], extra_fields={"gsd": 0.2}))
    item.add_asset("thumbnail", pystac.Asset(href="thumb.jpg", media_type="image/jpeg",
                                             roles=["thumbnail"]))
    item.set_self_href(f"./{item_id}.json")
    coll.add_item(item)          # 会同时写入 rel=item 链接

add_item(coll, "GD-SZ-2025Q1-ortho-0p2m", "2025-02-18T02:35:00Z",
         "GD-SZ-2025Q1-ortho-0p2m.tif", [113.75, 22.40, 114.62, 22.86])
coll.normalize_hrefs("./stac")   # 统一成相对路径
coll.save(catalog_type=pystac.CatalogType.SELF_CONTAINED)
print(coll.get_self_href())

生成的目录是静态的，直接扔对象存储就能被任何支持 STAC 的客户端检索，不需要部署服务。要素规模上千以后再考虑换 STAC API（stac-fastapi 一类），静态 catalog 的启动成本低，回滚也简单。

坑 13：把非 COG 的 GeoTIFF 直接挂到对象存储当影像服务。普通 GeoTIFF 的索引在文件末尾，客户端必须整景下载才能读到索引，range 请求等于没用。

坑 14：Item 的 id 随生成时间变化（如带时间戳后缀）。下游用它做外键，一改就断链；id 必须稳定且唯一。

坑 15：STAC 的 datetime 用本地时间且不带时区。所有客户端按 UTC 解析，检索区间会整体偏移 8 小时，季度边界上的景会漏检。

8.6 分析型底座：GeoParquet 与 DuckDB

PostGIS 擅长「事务 + 并发 + 空间谓词」，不擅长「扫一遍全表做聚合」。凡是「对一亿条记录做一次分组统计」这类任务，把数据导出成列式文件、用向量化引擎扫，比在数据库里硬跑快一个数量级。

8.6.1 GeoParquet 的版本选择

GeoParquet 在 Parquet 之上约定了几何列的存法与元数据：几何以 WKB 存在 BINARY 列里，文件元数据写明哪一列是几何、CRS 是什么、包围盒范围。这带来三个直接好处：列裁剪（只读需要的列）、谓词下推（按统计信息跳过整个 row group）、跨语言可读（Python、DuckDB、Spark 都认）。

选 1.1 而不是 2.0 的原因很实际：2.0 依赖 Parquet 侧的原生地理类型支持，读取端（各语言库、引擎、BI 工具）还没铺开，一旦上游改了类型编码，历史文件要重写。生产环境「能用 1.1 说清楚的事，不要用 rc 版本赌未来」。

写出 GeoParquet 用 GeoPandas 一行，写入时注意三件事：CRS 要显式设置、几何列只保留一列、按查询常用的字段排序：

import geopandas as gpd

gdf = gpd.read_file("data/parcel_2026q1.gpkg", layer="parcel")
gdf = gdf.to_crs(4490)                                  # 统一到 CGCS2000 地理
gdf = gdf.sort_values(["quality", "code"]).reset_index(drop=True)
gdf.to_parquet("warehouse/parcel/dt=2026-01-01/part-0.parquet",
               compression="zstd", compression_level=9, row_group_size=200_000)

分区目录用 Hive 风格（字段=值），DuckDB 能自动把分区键当列读进来，不需要额外声明。row_group_size 按 10 万–50 万行给，太小则元数据膨胀，太大则谓词下推粒度变粗。

8.6.2 DuckDB 1.5 + spatial：直查文件，不建库

DuckDB 1.5.x 的 spatial 扩展提供了 ST_Read、ST_AsMVT、ST_Area_Spheroid 这些函数，可以直接对文件做空间查询。千万级要素的统计任务，不需要导入、不需要建索引：

-- duckdb 命令行：duckdb -c ".read stats.sql"
INSTALL spatial;
LOAD spatial;

-- 1) 直查 Hive 分区：dt 从目录名自动成为列
CREATE OR REPLACE VIEW parcel AS
SELECT * FROM read_parquet('warehouse/parcel/dt=*/part-*.parquet',
                           hive_partitioning = true);

-- 2) 千万级统计：按质量等级求要素数与球面面积（公顷）
SELECT quality,
       count(*)                                                   AS n,
       round(sum(ST_Area_Spheroid(geom)) / 10000.0, 2)            AS area_ha,
       round(avg(ST_Perimeter(geom)), 1)                          AS avg_perim
FROM parcel
WHERE dt >= DATE '2026-01-01' AND dt < DATE '2026-04-01'
GROUP BY quality
ORDER BY quality;

-- 3) 只读需要的列：几何列不参与时，IO 直接省一半
SELECT count(*) FROM read_parquet('warehouse/parcel/dt=2026-01-01/part-0.parquet');

-- 4) 出矢量瓦片：按 z/x/y 切 MVT，直接交给前端
COPY (
  SELECT ST_AsMVT(t, 'parcel', 4096, 'geom') AS mvt
  FROM (SELECT code, quality, ST_AsMVTGeom(geom,
               ST_Extent(ST_MakeEnvelope(113.90, 22.50, 113.95, 22.55, 4490))) AS geom
        FROM parcel) t
) TO 'tiles/parcel.pbf' (FORMAT BLOB);

ST_Area_Spheroid 与 PostGIS 的 ST_Area(geom::geography) 口径一致，都是椭球面积，可以直接对账。这一步很重要：两个引擎算出来的面积必须能对齐到小数点后两位，否则口径差异会污染统计结论。

8.6.3 性能对照：什么情况下该换引擎

下面是参考量级，口径统一为「1000 万条多边形、8 个属性列，按质量等级分组求面积和」，同一台 16 核 / 64 GB 机器，冷缓存：

这些数字会随机器、数据分布、压缩率浮动，结论只看量级关系：扫描型聚合用列式文件快一个数量级，而带空间谓词的交互式查询（范围检索、最近邻）PostGIS 的 GiST 仍然不可替代。两者不是替代关系，是分工关系——热路径走 PostGIS，批量分析走 GeoParquet + DuckDB。

这个分工决定了底座的组织方式：同一份要素，入库 PostGIS 提供在线服务，同时定期导出 GeoParquet 到对象存储供分析。两份数据必须来自同一批源文件、同一套 CRS，用 8.8 的 manifest 保证「两处一致」，否则统计结果和服务结果会打架。

坑 16：把 GeoParquet 当「备份格式」整体导出，所有字段全量落盘。列裁剪的前提是文件里真的分了列，把几何、属性、原始 JSON 一股脑塞进一列，DuckDB 只能整行读，列存优势全没了。

坑 17：两份数据（PostGIS 与 GeoParquet）用不同的 CRS 或不同的面积口径。服务端用 4547 平面算面积，分析端用 4490 球面算，同一块地差出千分之几，统计结论直接对不上。

坑 18：分区字段选错。按采集日期分区却总按行政代码查询，分区裁剪永远不生效，DuckDB 要把所有文件头读一遍，第一次查询慢到让人以为数据库挂了。

8.7 三维底座：3D Tiles 1.1 与 splat 瓦片

三维底座要装两种东西，它们的评价标准完全不同：一是测绘成果（倾斜摄影、白模、DEM 加影像），要求几何准、可量测、带属性；二是近景真实感（3DGS），要求看着真、加载快，但不承担量测任务。混为一谈的后果是拿渲染成果去画竣工图。

8.7.1 两类内容，一种承载方式

3D Tiles 1.1 相对 1.0 有三个对底座有实际影响的变化：瓦片内容可以直接是 glTF（不再强制包一层 b3dm）、支持隐式瓦片（implicit tiling，用子树文件描述空间划分，替代巨大的 tileset.json 树）、支持元数据（metadata，把每栋楼的属性随瓦片一起走）。对一个区县级的倾斜摄影，隐式瓦片能把 tileset.json 从几百 MB 压到几 KB，客户端启动解析时间从分钟级降到毫秒级。

8.7.2 目录结构与 tileset.json

标准布局如下，四类文件各司其职，content 与 subtrees 的数量由隐式划分决定，不需要手写：

tiles/oblique-3dtiles/
├── tileset.json                 # 只有一个根瓦片 + 隐式划分参数
├── content/{level}/{x}/{y}.glb   # 每个瓦片的 glTF 内容
├── subtrees/{level}/{x}/{y}.subtree  # 可用性位图 + 子节点索引
└── metadata/schema.json         # 每瓦片 / 每要素的属性表结构

{
  "asset": {"version": "1.1"},
  "geometricError": 512.0,
  "root": {
    "boundingVolume": {"region": [1.9853, 0.3910, 2.0005, 0.3990, 0.0, 120.0]},
    "geometricError": 256.0,
    "refine": "ADD",
    "content": {"uri": "content/0/0/0.glb"},
    "implicitTiling": {
      "subdivisionScheme": "QUADTREE",
      "subtreeLevels": 6,
      "availableLevels": 18,
      "subtrees": {"uri": "subtrees/{level}/{x}/{y}.subtree"}
    }
  }
}

boundingVolume 用 region 时，前四个数是弧度：西经、南纬、东经、北纬，后两个是高程范围米。这里 [1.9853, 0.3910, 2.0005, 0.3990] 对应东经 113.75–114.62 度、北纬 22.40–22.86 度——弧度与度搞混是三维底座最常见的低级错误，症状是「地球上看不到模型」，因为包围盒被放到了地球内部。

属性随瓦片走，靠的是 metadata：把楼栋数、采集日期、数据版本挂在瓦片级，点击时不用回查数据库：

"metadata": {
  "class": "tile",
  "properties": {"bldgCount": {"minimum": 0, "maximum": 4200}, "surveyYear": {"minimum": 2022}}
}

8.7.3 splat tileset：真实感归真实感

Cesium 在 2026-04 发布了「3D Tiles + 层次 LOD」的 3D Gaussian Splats 方案，CesiumJS 1.147 可以直接加载带 LOD 的 splat tileset。它在底座里的位置很明确：作为「看得清」的补充图层，与「量得准」的测绘瓦片叠在同一个地球上，用图层开关切换。

加载代码与普通三维瓦片一致，关键参数是屏幕空间误差：

import * as Cesium from "cesium";

const viewer = new Cesium.Viewer("cesiumContainer", {
  terrain: Cesium.Terrain.fromWorldTerrain(),
});

const splat = await Cesium.Cesium3DTileset.fromUrl("/tiles/splat-clip/tileset.json", {
  maximumScreenSpaceError: 12,   // 越小越清晰、显存越高；首屏可先给 32 再降到 12
  dynamicScreenSpaceError: true, // 视野边缘降质，省 10–30% 显存
});
viewer.scene.primitives.add(splat);
viewer.zoomTo(splat, new Cesium.HeadingPitchRange(0.0, -0.35, 180.0));

8.7.4 首屏预算与显存估算

三维底座的失败多数不是「打不开」，是「在甲方那台笔记本上转不动」。预算要在建模阶段就定，而不是上线后优化：

估算口径：100 万三角面（含位置、法线、UV、索引）约占 40 MB 显存，纹理按压缩后尺寸的 1.33 倍计；3DGS 按每 splat 只存 DC 项约 56 字节、完整三阶球谐约 236 字节计，2000 万个 splat 在 DC 模式下约 1.1 GB。首屏策略统一为「先粗后细」：初始 maximumScreenSpaceError 给 32，视野稳定后降到 12–16，用户感知的是「秒开」，而不是「等最高精度」。

坑 19：把倾斜摄影的原始 OSGB 目录直接挂服务。几十万个小文件、没有统一的几何误差层级，客户端一次请求几百个文件，首屏半分钟起。

坑 20：geometricError 拍脑袋填。它与瓦片包围盒、屏幕空间误差是联动关系，填小了客户端死等最高层级，填大了永远加载不到细节；用转换工具生成，别手改。

坑 21：拿 3DGS 成果去量测。3DGS 是渲染表示，没有稳定的几何面，两点之间的「距离」随视角与训练质量变。量测只能用 mesh、点云或正射影像。

8.8 数据治理：版本、血缘、权限、合规

治理听起来抽象，落到操作上只有四个问题：这份数据是哪一版、从哪来、谁能看、能不能对外发。每个问题都要有可查询的答案，而不是靠某个人的记忆。

8.8.1 版本：manifest 加内容哈希

manifest 是一张表，一行描述一个成果文件；哈希用来判断「是不是同一份内容」，而不是「文件名一样不一样」。文件名会重复、会被复制粘贴，哈希不会：

CREATE TABLE gov.manifest (
  id           bigserial PRIMARY KEY,
  dataset_id   text        NOT NULL,            -- 逻辑数据集，例如 ortho-0p2m
  version      text        NOT NULL,            -- 语义版本，例如 2026.1
  rel_path     text        NOT NULL,            -- 对象存储内的相对路径
  sha256       char(64)    NOT NULL,            -- 内容哈希
  bytes        bigint      NOT NULL,
  crs          text        NOT NULL,            -- 例如 EPSG:4547
  producer     text        NOT NULL,            -- 生产工具与版本
  produced_at  timestamptz NOT NULL,
  source_ids   text[]      NOT NULL DEFAULT '{}',  -- 上游 manifest.id
  UNIQUE (dataset_id, version, rel_path)
);
CREATE INDEX manifest_ds_ver ON gov.manifest (dataset_id, version);
CREATE INDEX manifest_src    ON gov.manifest USING GIN (source_ids);

生成哈希与登记台账用一条命令串起来，写进流水线而不是手工执行：

# 1) 算哈希并追加到清单（find -print0 处理中文与空格文件名）
find output/ -type f -name "*.tif" -print0 \
  | xargs -0 sha256sum > manifest.sha256

# 2) 与上一版比对，只看差异，避免整批重传
sha256sum -c manifest.sha256 --quiet 2>&1 | head

# 3) 登记到台账（psql 的 \copy 支持从 stdin 读 CSV）
psql -d gis -c "\copy gov.manifest (dataset_id, version, rel_path, sha256, bytes, crs, producer, produced_at) FROM 'manifest.csv' WITH (FORMAT csv, HEADER true)"

复现一份三个月前的报告，流程是：从报告里读出 dataset_id 与 version → 查 manifest 拿到全部 rel_path 与 sha256 → 按清单拉回文件 → 校验哈希 → 重跑。任何一步对不上就中止，不猜。

8.8.2 血缘：从原始影像到服务图层

血缘回答「这个图层是怎么来的」。它是一张有向图，每条边记录一次加工，节点是 manifest 里的一行：

原始影像(2025-02-18 航摄)
   │ 空三 + 密集匹配          [tool: 空三软件 x.y]
   ▼
点云 / DSM
   │ 正射纠正 + 匀色           [tool: 正射流程 z]
   ▼
DOM(0.2 m) ────────────────┐
   │ COG 转换 + 概览           │ 入库
   ▼                          ▼
COG(2026.1) ──► PostGIS 矢量层 ──► 瓦片服务
   │                                   │
   └──────► GeoParquet 分析表 ─────────┘
                    │
                    ▼
              统计报告(2026-04)

落地方式就两条约束：每个产物必须记录 source_ids（上游 id 数组），每个加工步骤必须记录工具名与版本。有了这两样，任意一份结论都能沿图回溯到原始影像，中间缺一环就说明血缘断了。

8.8.3 权限：行级、列级与服务账号

空间数据天然带行政区划属性，权限也按区划切。PostgreSQL 的行级安全（RLS）能把这层约束放进数据库，而不是指望每个应用都记得加 WHERE：

-- 1) 行级：每个区县只能看到本区县的地块
ALTER TABLE base.parcel ENABLE ROW LEVEL SECURITY;

CREATE POLICY parcel_by_district ON base.parcel
  FOR SELECT
  USING (props->>'district' = current_setting('app.district', true));

-- 2) 会话建立时设置当前区县（应用连接池里 set_config 一次）
SET app.district = '440305';

-- 3) 列级：对外视图隐藏敏感字段
CREATE VIEW pub.parcel_public AS
SELECT id, code, geom, quality, acquired_at
FROM base.parcel;                 -- props 里的权属人信息不外露
REVOKE ALL ON base.parcel FROM PUBLIC;
GRANT SELECT ON pub.parcel_public TO svc_tiles;

原则只有一条：服务账号按「最小可用的 schema 加最小的动词」授予，绝不给 superuser。空间数据的破坏往往来自一次全表 UPDATE，而不是外部攻击。

8.8.4 合规：分级、分辨率红线与脱敏

测绘成果在国内是有级别的：涉及国家秘密的成果按密级管理，公开发布另有分辨率与要素的限制。这一块不能靠经验，必须按现行的《测绘法》《测绘成果管理条例》以及自然资源主管部门的最新规定执行；本章不展开条文，只给工程上要落的三件事。

第一，给每份成果打密级标签，写进 manifest 与对象存储的元数据，并在发布链路上做强制校验。第二，公开底图与内部底图分开存放，公开版按现行分辨率规定处理，处理动作本身留记录。第三，涉密要素（军事设施、重要管线、精确坐标的控制点）在出图、出瓦片、出 API 三个阶段分别脱敏，脱敏规则进代码库接受评审。具体阈值与要求逐年调整，上线前请以主管部门的现行文件为准。

8.8.5 治理检查表

坑 22：manifest 只登记「最终成果」，中间产物不登记。三个月后要复现，发现缺了正射纠正那一步的输入，血缘在这里断掉，只能整条重跑。

坑 23：用文件名当版本。加了「最终版」「真的最终版」后缀的文件名在两次交付后彻底失控；版本号进表，文件名保持稳定。

坑 24：公开底图与内部底图共用一份对象存储前缀，只靠应用层判断。一次配错权限就是整批影像对外可读；前缀与桶都要物理隔离。

8.9 用 DeepSeek Harness 维护底座

底座的问题有两个特点：种类固定、出现频繁。固定到可以写成脚本，频繁到必须自动化。这一节把前面八节的检查项收拢成一个「数据体检」技能，接进 CI，让底座在每次入库前自己证明干净。

8.9.1 体检脚本：五类问题一次查完

五类问题是长期统计下来最高频的：CRS 缺失、CRS 不在白名单（等同混合坐标系）、几何自相交、空几何、时间戳异常（未来时间、全空、重复）。脚本输出机器可读的 JSON，退出码用来说「过不过」：

# tools/data_health.py —— 底座体检，退出码 0 = 通过，1 = 有阻断项
import json
import sys

import geopandas as gpd
import pandas as pd
import pyogrio

CRS_WHITELIST = {"EPSG:4490", "EPSG:4547"}     # 项目允许的 CRS，按需收窄
TS_COLUMNS = ("acquired_at", "survey_time", "created_at")
NOW = pd.Timestamp.now(tz="UTC").tz_localize(None)   # 统一用 UTC 比较，避免本地时区偏移

def check_layer(path: str, layer: str) -> list[dict]:
    info = pyogrio.read_info(path, layer=layer)
    gdf = gpd.read_file(path, layer=layer, engine="pyogrio")
    issues: list[dict] = []

    crs = str(info.get("crs")) if info.get("crs") else None
    if crs is None:
        issues.append({"kind": "no_crs", "layer": layer, "detail": "未声明 CRS"})
    elif crs not in CRS_WHITELIST:
        issues.append({"kind": "crs_not_allowed", "layer": layer, "detail": crs})

    geom = gdf.geometry
    invalid = int((~geom.is_valid).sum())
    empty = int((geom.is_empty | geom.isna()).sum())
    if invalid:
        issues.append({"kind": "invalid_geometry", "layer": layer, "detail": invalid})
    if empty:
        issues.append({"kind": "empty_geometry", "layer": layer, "detail": empty})

    for col in TS_COLUMNS:
        if col not in gdf.columns:
            continue
        ts = pd.to_datetime(gdf[col], errors="coerce", utc=True).dt.tz_localize(None)
        future = int((ts > NOW).sum())
        missing = int(ts.isna().sum())
        if future or missing:
            issues.append({"kind": "bad_timestamp", "layer": layer,
                           "detail": {"column": col, "future": future, "missing": missing}})
    return issues

def main(argv: list[str]) -> int:
    path = argv[1]
    report = {"file": path, "layers": 0, "issues": []}
    for layer in pyogrio.list_layers(path)[:, 0]:
        report["layers"] += 1
        report["issues"].extend(check_layer(path, layer))
    blocking = [i for i in report["issues"]
                if i["kind"] in {"no_crs", "crs_not_allowed", "invalid_geometry"}]
    report["blocking"] = len(blocking)
    print(json.dumps(report, ensure_ascii=False, indent=2))
    return 1 if blocking else 0

if __name__ == "__main__":
    sys.exit(main(sys.argv))

8.9.2 接进 CI：入库前的门禁

把体检脚本挂成 CI 的一个 job，用 DSH 做「判断与摘要」，让脚本保持纯粹的规则执行：

#!/usr/bin/env bash
set -euo pipefail

# 1) 规则检查：失败即阻断流水线
python tools/data_health.py build/parcel_2026q1.gpkg > build/health.json

# 2) 只读模式跑诊断，产出可读摘要与修复建议；--json 便于后续消费
dsh --profile headless \
  "读取 build/health.json，按严重程度分组，输出三段：阻断项与修复顺序、告警项与业务影响、本次体检的通过判据。不要修改任何文件。" \
  --permission-mode read-only --json > build/health_summary.json

门禁的判据要写清楚：阻断项为 0 才允许进入建索引与发布；告警项不为 0 时必须在工单里留一条处置记录。这样「数据脏了」不再是口头提醒，而是流水线上一个红叉。

8.9.3 沉淀成技能：下一次不用重新讲

把体检流程写进 SKILL.md，任何人接手都能一次跑对：

---
name: data-health
description: 空间数据底座体检——CRS、几何、时间戳、manifest 四项，输出阻断/告警分组
---

1. 运行 python tools/data_health.py <数据文件>，把 JSON 结果保存到 build/health.json
2. 检查 gov.manifest：dataset_id 与 version 是否齐全，sha256 是否为空
3. 对阻断项输出：问题类型、涉及图层、最小复现命令（不超过 10 行）
4. 对告警项输出：业务影响、是否可带病发布、建议的工单标题
5. 不改任何数据文件；需要修复时只给 SQL 或命令，由人确认后执行

坑 25：把体检写成一次性脚本，跑完就扔。底座的问题会随数据量线性增长，规则必须跟着数据走，沉淀成技能与 CI job 才有复利。

坑 26：体检只查几何、不查时间。时间戳的「未来时间」与「全空」在入库当时看不出来，等到做时序分析时才发现整批记录都不可用。

坑 27：让模型直接改数据来「修好」检查。修复动作必须由人确认、由脚本执行、留执行记录；模型的职责是判断与汇总，不是拿 UPDATE 去碰生产表。

8.10 本章小结

本章的九节可以压缩成三句话：

第一，基准先于一切。坐标基准决定「在哪」，时间基准决定「何时」，两者都要用数值判据验收——平面 1 cm、时间 20 ms。基准没钉死之前，建库和展示都是白做。

第二，组织要按访问方式设计。事务型空间查询交给 PostGIS（SRID 进类型、GiST 加分区）；栅格交给 COG 加 STAC（range 请求加可检索目录）；批量分析交给 GeoParquet 加 DuckDB（列裁剪加向量化）；三维交给 3D Tiles 1.1（隐式瓦片加 metadata），近景真实感单独用 splat 瓦片承载，不承担量测。

第三，治理要落到表和流水线。manifest 加哈希解决「哪一版」，血缘加 source_ids 解决「从哪来」，RLS 加最小权限解决「谁能看」，密级标签加发布校验解决「能不能发」。全部做成检查项，才有可执行的约束力。

8.11 动手练习

五个练习按「基准 → 时间 → 矢量 → 影像 → 分析」排，每个都能在半天内做完，验收全部是数值而非常见评语。

练习 1：坐标往返链与自检

输入：一个测区的实测控制点 CSV（至少 5 个点，含 CGCS2000 经纬度与高斯平面坐标），以及该测区的地方独立坐标四参数。

要求：实现 WGS84 → CGCS2000 高斯 → 地方独立系 → ENU → 回到 WGS84 的完整链路；实现 8.2.5 的三条不变量检查；把往返残差与控制点残差都打出来。

验收：三条不变量全部通过；往返残差小于 0.1 mm；控制点平面残差小于 1 cm、高程残差小于 3 cm；把容差改成 0.5 cm 后脚本必须失败——证明它真的在判。

练习 2：时间偏差检测与判定

输入：两路带时间戳的数据（一份视频的逐帧 PTS 加起录 UTC，一份 10 Hz 的 POS 日志）。没有实测数据时，用脚本合成：给 POS 加 1.96 s 固定偏移与 8 ppm 漂移。

要求：写出检测脚本，做线性回归分离固定偏移与相对漂移，输出 RMS、P95、offset、drift 四个数；写清楚 offset 与 drift 各对应什么工程问题。

验收：合成数据上测得 offset 与 1.96 s 相差小于 1 ms，drift 与 8 ppm 相差小于 0.2 ppm；判据在 20 ms 处能正确 PASS / FAIL；残差图上人为插入一个 100 ms 台阶后，脚本能指出台阶位于第几帧。

练习 3：PostGIS 建库与查询优化

输入：一份 100 万条以上的面要素数据（可用公开行政区划数据），PostgreSQL 17 加 PostGIS 3.7.0。

要求：按 8.4.3 建分区表与三类索引；导入并 ANALYZE；写出范围检索、最近邻、分组聚合三条 SQL；对每条查询给出 EXPLAIN (ANALYZE, BUFFERS) 输出与判读。

验收：范围检索走 GiST（计划里出现 Index Scan 或 Bitmap Index Scan）；最近邻走 <-> 排序且 LIMIT 在 10 时耗时小于 50 ms；分区裁剪生效（计划里只出现 1 个分区）；把 WHERE 里的 && 换成只写 ST_Intersects 后耗时明显上升，能解释原因。

练习 4：把一景影像做成可检索集合

输入：一景正射影像（GeoTIFF，含 CRS）。

要求：用 rio-cogeo 转成 COG 并通过 validate；用 pystac 生成 Collection 与一个 Item，字段按 8.5.2 写全；在本地起一个静态文件服务，用 curl 只取前 16 KB 验证 range 请求可用。

验收：rio cogeo validate 输出通过；Item 的 datetime 是带 Z 的 UTC；proj:epsg 与文件真实 CRS 一致；三次 range 请求的响应时间均在 300 ms 内，且 Content-Range 头正确。

练习 5：千万级要素统计与两引擎对账

输入：把练习 3 的数据导出为 GeoParquet 1.1（Hive 分区），以及同一个库里的 PostGIS 表。

要求：用 DuckDB 1.5 加 spatial 扩展完成按质量等级的分组统计（要素数、球面面积、平均周长）；在同一台机器上跑 PostGIS 的等价 SQL；把两份结果逐组比对。

验收：两个引擎的面积差值相对误差小于 1e-6（口径对齐的证明）；DuckDB 在 1000 万行上的耗时显著低于 PostGIS 且能给出自己的量级；写出「什么查询该走哪个引擎」的一页结论。

8.12 延伸阅读

1. PROJ 官方文档：坐标转换引擎与参数定义，遇到 EPSG 编码拿不准时的第一手资料。 proj.org
2. GDAL 官方文档：矢量与栅格的读写、投影转换、COG 相关的命令行与 API 参考。 gdal.org
3. PostGIS 官方手册：几何类型、索引、空间谓词与性能调优建议的权威来源。 postgis.net/docs
4. COG 官方说明：Cloud Optimized GeoTIFF 的布局约束与校验方法。 cogeo.org
5. STAC 规范与生态：Catalog / Collection / Item 三层模型与最佳实践。 stacspec.org
6. GeoParquet 规范：几何列的存储约定、元数据字段与版本演进。 geoparquet.org
7. DuckDB 与 spatial 扩展：列式分析引擎的 SQL 方言、扩展安装与空间函数。 duckdb.org
8. 3D Tiles 规范（CesiumGS）：隐式瓦片、metadata 与瓦片内容格式。 github.com/CesiumGS/3d-tiles