这篇能做出什么
跑完这篇的步骤，你会得到一条可以反复重跑的流水线：给定一片研究区和一批公开 LiDAR 点云，自动切成上千个瓦片，由 48 个并发智能体在 72 小时内分头扫描，产出带置信度的候选标记（建筑物轮廓异常、地面沉降疑似区、植被高度突变、人工地物新增等），经过多角色交叉验证后，把“有分歧”和“高价值”的瓦片排进人工复核队列。
最终目录大致长这样：
```text
scan_run/
tiles.jsonl # 瓦片清单（核心区 + 带 buffer 的读取范围）
features/ # 每个瓦片的统计特征 JSON
renders/ # 每个瓦片的俯视渲染 PNG（给智能体看）
queue.db # 任务状态机（SQLite）
results/ # 每个智能体的结构化输出
verdicts.jsonl # 交叉验证后的结论
review_queue.csv # 人工复核队列（已排序）
report.html # 一页式汇总报告
```
关于“48 智能体”要先澄清一件事：它不是 48 个不同的大模型，而是 48 个并发执行单元。每个单元是“一个模型 + 一套角色提示 + 一个瓦片”的组合。这样做的好处是隔离、可重试、可审计，出问题能定位到具体是哪个 role 在哪个 tile 上错了。
粗略的容量估算：48 并发 × 72 小时 = 3456 机时。若单次任务平均耗时 5–10 分钟（含重试、限流等待、渲染），总量级在 2 万到 4 万次调用之间；每个瓦片被 3 个角色各看一遍，那么可覆盖数千到一万个瓦片。实际数取决于你的单任务耗时和限流策略，先拿几十个瓦片计时再外推，不要拍脑袋。
前置条件清单
- 一台能连续跑 72 小时的机器：云主机或工作站均可，重点是磁盘。点云切片和渲染图很占空间，按研究区规模预留，建议先算一遍“瓦片数 × 单瓦片渲染图大小”再决定。
- Python 环境：Python 3.9 以上，依赖
laspy、numpy、pandas、pyproj、matplotlib，可选pdal（命令行裁切更快）。版本以各项目官方文档当前版本为准，装最新稳定版即可。 - 一个支持图像输入的模型 API：有速率限制和并发上限，需要看服务商官方文档的配额说明。预算上按“调用次数上限”而非金额做闸门更稳。
- 公开 LiDAR 数据：USGS 3DEP、OpenTopography、国内各省级自然资源部门的开放点云门户、国家青藏高原科学数据中心等都可以。下载方式、坐标系、点云分类字段（是否含地面点标注）各家不同，以各平台官方页面当前版本为准。下载前先读许可条款，确认允许二次加工和成果发布。
- 一个明确的判定问题：这一步最容易被跳过。不要写“分析点云质量”，要写成“这块瓦片里是否存在高于局部地面 2 米以上、面积超过 50 平方米、且边界呈直线的点簇”。智能体不会替你定义什么算异常。
分步骤
第 1 步：把研究区切成带 buffer 的瓦片
切片的两个关键参数是瓦片边长和边界缓冲。缓冲的作用是防止跨瓦片的线性地物（道路、堤坝、围墙）被切断后，两边都判成“无目标”。
```python
01_make_tiles.py
import json, math
from pyproj import Transformer
AOI = {"west": -105.30, "south": 39.90, "east": -105.00, "north": 40.10} # 换成你的研究区
SRC_EPSG = "EPSG:4326"
UTM_EPSG = "EPSG:32613" # 按研究区所在 UTM 带号修改，务必先确认
TILE_M = 1000.0 # 瓦片核心区边长（米）
BUFFER_M = 50.0 # 边界缓冲（米）
OUT = "tiles.jsonl"
fwd = Transformer.from_crs(SRC_EPSG, UTM_EPSG, always_xy=True)
inv = Transformer.from_crs(UTM_EPSG, SRC_EPSG, always_xy=True)
x0, y0 = fwd.transform(AOI["west"], AOI["south"])
x1, y1 = fwd.transform(AOI["east"], AOI["north"])
nx = math.ceil((x1 - x0) / TILE_M)
ny = math.ceil((y1 - y0) / TILE_M)
with open(OUT, "w", encoding="utf-8") as f:
for i in range(nx):
for j in range(ny):
cx0, cy0 = x0 + i * TILE_M, y0 + j * TILE_M
cx1, cy1 = cx0 + TILE_M, cy0 + TILE_M
bx0, by0 = cx0 - BUFFER_M, cy0 - BUFFER_M
bx1, by1 = cx1 + BUFFER_M, cy1 + BUFFER_M
w, s = inv.transform(bx0, by0)
e, n = inv.transform(bx1, by1)
f.write(json.dumps({
"tile_id": f"T{i:04d}_{j:04d}",
"crs": UTM_EPSG,
"core_bbox": [cx0, cy0, cx1, cy1],
"read_bbox": [bx0, by0, bx1, by1],
"read_bbox_wgs84": [w, s, e, n],
}, ensure_ascii=False) + "\n")
print("total tiles:", nx * ny)
```
UTM_EPSG 一定要和研究区匹配。带号错了，切出来的网格会整体偏移，后面所有坐标都对不上。
第 2 步：把点云变成智能体能读的东西
不要指望模型直接读 .las。稳妥做法是每个瓦片产出两样东西：一份统计特征 JSON，一张固定尺寸的俯视渲染图。模型看这两样做判断，判定逻辑才可复现。
```python
02_extract_features.py
import json, os
import numpy as np
import laspy
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
CELL = 1.0 # 栅格分辨率（米）
def summarize(las_path, tile):
las = laspy.read(las_path)
x = np.asarray(las.x, dtype="float64")
y = np.asarray(las.y, dtype="float64")
z = np.asarray(las.z, dtype="float64")
if z.size < 50:
return None, None
area = max((x.max() - x.min()) * (y.max() - y.min()), 1.0)
ground = np.percentile(z, 1) # 粗估局部地面；有分类地面点时改用 classification==2
above = z - ground
feat = {
"tile_id": tile["tile_id"],
"n_points": int(z.size),
"density_per_m2": round(z.size / area, 3),
"z_min": round(float(z.min()), 2),
"z_max": round(float(z.max()), 2),
"z_p50": round(float(np.percentile(z, 50)), 2),
"z_p99": round(float(np.percentile(z, 99)), 2),
"z_std": round(float(z.std()), 2),
"above_ground_2m_ratio": round(float((above > 2.0).mean()), 4),
"above_ground_5m_ratio": round(float((above > 5.0).mean()), 4),
"intensity_p50": round(float(np.percentile(np.asarray(las.intensity), 50)), 2),
}
俯视最大值渲染，固定像素尺寸，便于模型稳定阅读
gx = ((x - x.min()) // CELL).astype(int)
gy = ((y - y.min()) // CELL).astype(int)
grid = np.full((gy.max() + 1, gx.max() + 1), -np.inf)
np.maximum.at(grid, (gy, gx), above)
grid[~np.isfinite(grid)] = np.nan
fig, ax = plt.subplots(figsize=(6, 6), dpi=100)
im = ax.imshow(grid, cmap="inferno", origin="lower")
ax.set_title(tile["tile_id"]); ax.set_xticks([]); ax.set_yticks([])
fig.colorbar(im, ax=ax, shrink=0.7, label="height above ground (m)")
png_path = os.path.join("renders", tile["tile_id"] + ".png")
fig.savefig(png_path, bbox_inches="tight"); plt.close(fig)
return feat, png_path
if __name__ == "__main__":
os.makedirs("features", exist_ok=True)
os.makedirs("renders", exist_ok=True)
for line in open("tiles.jsonl", encoding="utf-8"):
tile = json.loads(line)
las_path = f"las/{tile['tile_id']}.las" # 上一步裁好的瓦片点云
feat, _ = summarize(las_path, tile)
if feat is None:
continue
with open(f"features/{tile['tile_id']}.json", "w", encoding="utf-8") as f:
json.dump(feat, f, ensure_ascii=False)
```
裁切点云用 PDAL 的 readers.las + filters.crop 组合即可，命令按 pdal --help 和官方文档写，不用记。
第 3 步：建任务队列与状态机
72 小时长跑一定会中断。不要用“内存里的循环”，要用能断点续跑的任务表。SQLite 够用且零依赖。
```sql
-- 03_schema.sql
CREATE TABLE IF NOT EXISTS tasks (
tile_id TEXT NOT NULL,
role TEXT NOT NULL, -- structure / terrain / change
priority INTEGER NOT NULL DEFAULT 0,
attempt INTEGER NOT NULL DEFAULT 0,
state TEXT NOT NULL DEFAULT 'pending', -- pending/running/done/failed
worker TEXT,
lease_until TEXT,
PRIMARY KEY (tile_id, role)
);
CREATE TABLE IF NOT EXISTS results (
tile_id TEXT NOT NULL,
role TEXT NOT NULL,
worker TEXT NOT NULL,
payload TEXT NOT NULL, -- 智能体输出的 JSON 原文
created_at TEXT NOT NULL DEFAULT (datetime('now')),
PRIMARY KEY (tile_id, role, worker)
);
```
领取任务用“租约”而不是“标记为运行中”，这样宕机的 worker 留下的任务会自动过期被重新领取：
```python
03_claim.py
def claim(conn, worker_id, lease_seconds=1800):
cur = conn.execute("""
UPDATE tasks
SET state='running',
worker=?,
attempt=attempt+1,
lease_until=datetime('now', ?)
WHERE rowid = (
SELECT rowid FROM tasks
WHERE state='pending'
OR (state='running' AND lease_until < datetime('now'))
ORDER BY priority DESC, attempt ASC, tile_id
LIMIT 1)
RETURNING tile_id, role, attempt
""", (worker_id, f"+{lease_seconds} seconds"))
return cur.fetchone()
```
RETURNING 需要较新的 SQLite；若你的环境不支持，改写成 BEGIN IMMEDIATE 包住“先 SELECT 再 UPDATE”两步，效果一样。
第 4 步：定死智能体的输入输出契约
契约不清楚，交叉验证就没法做。输入固定为“特征 JSON + 渲染图 + 一句判定标准”，输出固定为一个 JSON 对象。
```python
04_agent.py
import json, re
SCHEMA_HINT = """只输出一个 JSON 对象，不要输出任何解释性文字：
{
"tile_id": "字符串",
"has_candidate": true 或 false,
"candidates": [
{"label": "building_edge|ground_subsidence|vegetation_change|unknown",
"confidence": 0.0 到 1.0 之间的小数,
"bbox_px": [x0, y0, x1, y1],
"reason": "不超过 60 字，引用你看到的数值或图斑位置"}
],
"quality_flags": ["low_density", "edge_clipped", "no_ground_class"],
"notes": "不超过 80 字"
}"""
ROLE_PROMPTS = {
"structure": "你只关注几何形状：直线边界、规则多边形、直角转折。忽略颜色与整体亮度。",
"terrain": "你只关注地形连续性：地面是否出现局部凹陷、台阶、坡度突变。忽略单体建筑。",
"change": "你只关注与周边瓦片统计分布明显不一致的地方，先看 above_ground_2m_ratio 与 z_std。",
}
def build_messages(feat, png_b64, role):
text = f"""判定标准：瓦片内是否存在高于局部地面 2 米、水平尺度超过 50 平方米、
且不是孤立噪点的目标簇。
瓦片统计：
{json.dumps(feat, ensure_ascii=False, indent=2)}
角色约束：{ROLE_PROMPTS[role]}
{SCHEMA_HINT}"""
return [{"role": "user", "content": [
{"type": "text", "text": text},
{"type": "image_url", "image_url": {"url": f"data:image/png;base64,{png_b64}"}},
]}]
def parse_json_block(text):
m = re.search(r"\{.*\}", text, re.S)
if not m:
raise ValueError("no json found")
obj = json.loads(m.group(0))
for key in ("tile_id", "has_candidate", "candidates", "quality_flags"):
if key not in obj:
raise ValueError(f"missing key: {key}")
return obj
```
调用部分按你所用服务商的 SDK 写，注意三件事：把 temperature 调低以保证可复现；设置单次请求超时；失败重试最多 3 次，第 3 次仍失败就把任务置为 failed 而不是无限重试。
第 5 步：把 48 个 worker 跑起来
每个 worker 是一个独立进程，循环“领任务 → 取特征和渲染图 → 调模型 → 校验 JSON → 落库 → 标记完成”。三件事必须做：
- 幂等：
results表主键里带worker，重复写入不会互相覆盖；同一个(tile_id, role)已有结果就跳过。 - 限流：全局用信号量控制在服务商允许的并发以内，遇到限流错误按指数退避，不要所有 worker 同时重试。
- 日志：每个 worker 单独一个日志文件，记录
tile_id / role / attempt / 耗时 / token 用量 / 错误类型。72 小时后你要靠它算账。
```bash
启动（示例，进程管理用 systemd 或 supervisor 更稳）
for i in $(seq 1 48); do
nohup python 04_worker.py --worker w$i >> logs/w$i.log 2>&1 &
done
```
启动前先用 20 个瓦片跑一轮全流程，量出单任务耗时，再乘 48 推总时长。这一步能省掉后面 60 小时的心慌。
第 6 步：交叉验证，把分歧挑出来
三个角色对同一瓦片各看一遍，然后比对候选集合。比对前要先把“同一个候选”对齐：同一 label 且 bbox_px 的 IoU 超过阈值，就算同一处。
```python
05_crosscheck.py
def iou(a, b):
ax0, ay0, ax1, ay1 = a; bx0, by0, bx1, by1 = b
ix = max(0, min(ax1, bx1) - max(ax0, bx0))
iy = max(0, min(ay1, by1) - max(ay0, by0))
inter = ix * iy
union = (ax1-ax0)*(ay1-ay0) + (bx1-bx0)*(by1-by0) - inter
return inter / union if union > 0 else 0.0
def agree(ra, rb, iou_thr=0.3):
hits = 0
for ca in ra["candidates"]:
for cb in rb["candidates"]:
if ca["label"] == cb["label"] and iou(ca["bbox_px"], cb["bbox_px"]) >= iou_thr:
hits += 1
break
denom = max(len(ra["candidates"]), len(rb["candidates"]), 1)
return hits / denom
def verdict(results_by_role):
roles = list(results_by_role)
scores = []
for i in range(len(roles)):
for j in range(i + 1, len(roles)):
scores.append(agree(results_by_role[roles[i]], results_by_role[roles[j]]))
consistency = sum(scores) / len(scores) if scores else 0.0
votes = sum(1 for r in results_by_role.values() if r["has_candidate"])
flags = {f for r in results_by_role.values() for f in r["quality_flags"]}
if votes == len(roles) and consistency >= 0.6:
tier = "auto_accept"
elif votes >= 2:
tier = "review_low"
elif votes == 1:
tier = "review_high"
else:
tier = "auto_reject"
if "low_density" in flags or "edge_clipped" in flags:
tier = "review_high" # 数据质量可疑一律上浮
return {"consistency": round(consistency, 3), "votes": votes, "tier": tier, "flags": sorted(flags)}
```
这里有个反直觉的坑：如果三个角色全说“无候选”，一致率会算成 1.0，看着很漂亮，其实什么都没检出。所以一致率不能单独作为质量指标，必须和“投票数分布”“人工抽检的召回率”一起看。稳妥做法是每 200 个瓦片随机抽 10 个人工标注，算一次真实召回。
第 7 步：人工复核队列
复核不是“全部重看一遍”，那会直接把人压垮。按优先级排序，只把该看的送到人眼前：
```python
06_export_review.py
import csv, json
rows = []
for line in open("verdicts.jsonl", encoding="utf-8"):
v = json.loads(line)
if v["tier"] not in ("review_high", "review_low"):
continue
rows.append({
"tile_id": v["tile_id"],
"tier": v["tier"],
"votes": v["votes"],
"consistency": v["consistency"],
"flags": ";".join(v["flags"]),
"render": f"renders/{v['tile_id']}.png",
"candidates": json.dumps(v["candidate_summary"], ensure_ascii=False),
"verdict": "", # 人工填：confirm / reject / unclear
"note": "",
})
rows.sort(key=lambda r: (r["tier"] != "review_high", -r["votes"], r["consistency"]))
with open("review_queue.csv", "w", encoding="utf-8-sig", newline="") as f:
w = csv.DictWriter(f, fieldnames=list(rows[0].keys()))
w.writeheader(); w.writerows(rows)
print("review items:", len(rows))
```
复核界面不需要复杂系统：一个本地 HTTP 服务，左边显示 PNG，右边放三个按钮（确认 / 否决 / 存疑），提交后写回 CSV。关键是复核结果要回流：把人工结论按 tile_id 存下来，下次跑同区域时优先用人工标签校准提示词和阈值。
推荐的分层策略：review_high 全查，review_low 按 30% 抽样查，auto_accept 按 2% 抽查做质控。抽样比例按你的人力和总瓦片数调。
第 8 步：72 小时运行手册
把长跑当成值班，不是当成一次脚本执行。可以按这个节奏安排：
- 每 6 小时看一次：任务表里
pending / running / done / failed四个状态的计数，失败率突增通常意味着限流或密钥问题。 - 每 12 小时看一次：新增
verdicts的 tier 分布。如果auto_accept突然变成 0，多半是某个 role 的提示词或图片渲染出问题了。 - 每天一次人工抽检 10 个瓦片，和智能体结论对照，估算当前召回。
- 告警只设两条：失败率超过阈值、队列剩余时间超过剩余窗口。其余噪音全部忽略。
中断恢复很简单：重启 worker 进程即可，租约过期的任务会被自动重新领取，results 里已有的不会重算。
常见坑与排错
瓦片边界漏检。 表现为线性地物在瓦片接缝处断开。检查 read_bbox 是否真的比 core_bbox 大了一圈，以及裁切时用的是 read_bbox 而不是 core_bbox。
坐标系混用。 特征里的坐标是 UTM，渲染图是像素坐标，复核时的 bbox_px 又是另一套。建议全程只在渲染时做一次坐标换算，其余环节统一用 tile_id + 像素坐标，并在结果里显式记录 crs 字段。
模型输出不是合法 JSON。 三种处理：提示里明确“只输出 JSON 对象”；解析时先剥掉 ``json 围栏；解析失败重试一次并附上错误信息。三次都失败就落 failed`，别在原地打转消耗配额。
一致率虚高。 见第 6 步的说明。补一个动作：统计每个 role 报“有候选”的比例，如果某个 role 常年报 0，它就是在给一致率注水，需要单独检查。
智能体之间互相“抄答案”。 如果把 A 角色的结论放进 B 角色的提示里，交叉验证就失效了。三个角色的输入必须完全隔离，只共享特征和渲染图。
渲染图尺寸不统一。 不同瓦片的图大小不一，会让模型表现飘忽。固定 figsize 和 dpi，bbox_inches="tight" 之外再加一次 pad_inches=0，保证输出像素尺寸一致。
复核变瓶颈。 如果 review_high 占了 40% 以上，说明判定标准太模糊或提示词太松。先回去改判定标准，不要靠加人硬扛。
磁盘写满。 渲染图和中间点云最容易吃满盘。跑之前 df -h 看一眼，并在 worker 里加一条检查：剩余空间低于阈值就停止领新任务。
下一步建议
引入主动学习。 人工复核的结果不要只存着，用它训练一个轻量分类器，对 review_low 做预排序：模型不确定的排前面，人只看最有信息量的那批，同样人力能覆盖更多瓦片。
规则与模型混合。 纯规则（高程阈值 + 形态学）跑得快且完全可复现，适合做第一道粗筛；把规则筛出来的边界样本交给智能体判断，成本和准确率都更划算。
结果版本化。 每次跑完把 tiles.jsonl、提示词、阈值、verdicts.jsonl 打成一个带日期的包。同一区域换参数重跑时能直接对比，也方便回答“这个结论是怎么来的”。
扩到大区域前先压测。 把研究区扩大到 10 倍，瓶颈通常不在模型调用，而在切片、渲染和磁盘 IO。先单独跑一遍“只切片不调用模型”的流程，量出每小时能产出多少瓦片，再决定并发数。
最后提醒一句：这套流水线的产出是“候选标记”，不是“结论”。人工复核那一环不是形式，是整条链路里唯一能给出可信判断的地方，别把它省掉。