# 实验 16：Z 序曲线（Z-order / Morton code） > 把二维的经纬度压成一维的整数，让"空间上靠近的点"在一维顺序里也尽量靠近，从而用普通的 B 树索引就能做范围查询。 ## 一、实验目标 1. 理解**空间填充曲线**（space-filling curve）的作用：把多维数据线性化（linearization）成一维。 2. 掌握 Z 序码（z-code / z-value / Morton code）的构造：位交织（bit interleaving），以及它的逆变换。 3. 会用字符串形式的 z 码做教学演示，也会用 numpy 写向量化的整数 Morton 码。 4. 理解"用 z 码做范围查询"的两阶段流程：先用 z 区间取候选集（candidate），再精确过滤（refine）。 5. 用仓库自带的真实经纬度数据实测：顺序扫描、Z 序索引、Hilbert 序索引三种方案的候选率。 6. 知道 Z 序在实际系统中的落地形态（Geohash、S2、H3）与相关争议。 ## 二、预备知识 ### 2.1 为什么要把多维压成一维 数据库的 B+ 树索引只能对**一个**排序键做范围查找。经纬度是二维的，"找出某个矩形框内的点"没法直接用一个 B+ 树区间表达。解决办法是找一个映射： ``` f: (x, y) → z（一个整数） ``` 使得空间上邻近的点，其 z 值也尽量邻近。这样"矩形框查询"就近似变成"若干个 z 区间查询"，可以复用现成的 B+ 树。 ### 2.2 空间填充曲线家族 常见的空间填充曲线（space-filling curve，SFC）： - **Z 序**（Z-order / Morton order）：本实验的主角，构造最简单，只需位交织。 - **Gray 码变体**（Gray-coded）：把 Z 序的码字再做一次格雷码变换。 - **Hilbert 曲线**：递归的"U"形拼接，通常被认为局部性更好（但有争议，见 2.5）。 - **Peano 曲线 / Sierpinski 曲线**：其他递归构造。 ### 2.3 位交织规则 设递归层级（recursion level）为 L，则 x、y 各用 L 个二进制位表示，z 码共 2L 位。**本实验采用仓库 notebook 的位序约定**： ``` z 的位（从高到低）= y_{L−1} x_{L−1} … y_1 x_1 y_0 x_0 ``` 即：y 的最高位在前，x 的最低位在最后。举例（L = 3，x = 2，y = 1）： ``` x = 2 = 010 y = 1 = 001 y2 x2 y1 x1 y0 x0 z = 0 0 0 1 1 0 = "000110" ``` 网格大小与码长的关系： | 递归层级 L | 网格 | 单元数 | z 码位宽 | |---|---|---|---| | 1 | 2 × 2 | 4 | 2 | | 2 | 4 × 4 | 16 | 4 | | 3 | 8 × 8 | 64 | 6 | | L | 2^L × 2^L | 2^(2L) | 2L | ### 2.4 量化（quantization） 真实经纬度是浮点数，做位交织前必须先**量化**成整数网格坐标。本实验用线性映射（L = 16）： ``` x = round((lon + 180) / 360 × (2^L − 1)) y = round((lat + 90) / 180 × (2^L − 1)) ``` 量化层级越高，精度越高，但 z 码越长；量化本身就是一次精度损失，这也是空间索引共同的取舍。 ### 2.5 关于"Hilbert 一定比 Z 好"的说明 必须说清楚：**Hilbert 与 Z 序的局部性优劣在文献中存在不一致结论**。Dai & Su（DMTCS 2003）在"平均跨簇距离"与"平均总跨簇距离"两项统计上得到的是 **Z 序家族优于 Hilbert 家族**的结论。因此"Hilbert 一定更好"是有争议的说法，取决于用什么度量口径。工程上应当**基于自己的负载实测**来选择，而不是照搬结论。 ### 2.6 范围查询的两阶段 ``` ① 候选阶段（candidate）：把查询矩形翻译成一段（或几段）z 区间，用 B+ 树取出区间内所有点 ② 精化阶段（refine） ：对候选点逐个用真实坐标判断是否在矩形内 ``` 第 ① 步不能保证取到的点都在框内（z 区间是矩形的**超集**），所以第 ② 步不可省略。本实验用最简单的做法：把矩形左下角与右上角的 z 值作为一个区间 `[z_lo, z_hi]`。工业实现会用 `BIGMIN`/`LITMAX` 把大区间切成若干段，跳过中间大量不属于矩形的空洞。 衡量指标：**候选率** = 候选点数 / 真实命中点数。候选率越接近 1 越好；顺序扫描的候选率 = 总点数 / 命中点数。 ## 三、实验原理 ### 3.1 Z 曲线的形状 下图左半是 L = 3（8×8 网格）的 Z 曲线，右半是同网格下的 Hilbert 曲线。Z 曲线的每一级递归都在四个象限里重复一个"Z"字。 ![Z 曲线与 Hilbert 曲线](../figures/fig_zorder_curve.png) Z 曲线的关键性质（也是它的缺点）：相邻的两个 z 值在空间上几乎总是相邻的（局部性好），但在象限边界上会出现**大跳跃**（例如从 (3,0) 跳到 (0,1)）。 ### 3.2 用真实数据测局部性 实验数据来自仓库自带的 `notebooks/latitude_longitude.bin`： - 197,229 个点 - 纬度范围 47.28° – 55.07°，经度范围 −2.01° – 50.51° - 存储格式：float32 交替存放（纬度、经度、纬度、经度……），小端序 实验做法：随机生成 400 个小窗口矩形查询框，只统计真实命中数 ≥ 5 的查询（命中太少的查询不足以体现差异），记录每种方案的"候选数 / 真实命中数"的中位数。 ![真实经纬度上的局部性对比](../figures/fig_geo_locality.png) 本机实测结果（中位数）： | 方案 | 候选数 / 真实命中数 | |---|---| | 顺序扫描 | 12.64× | | Z 序索引 | 3.19× | | Hilbert 序索引 | 1.77× | （89 个有效查询窗口；Z 序与 Hilbert 序均建立在**同一个** 2^16 网格上。） **为什么要强调"网格一致"。** 本案例库早期版本的图件里，Z 码用 2^16 网格、Hilbert 码用 2^10 网格。网格越粗，同一个查询框在码空间里覆盖的区间越小，候选集自然更少——那时测出的 1.29× 有相当一部分来自量化层级差异，是**假优势**，不是曲线本身的性质。现在两者统一到 2^16，Hilbert 的领先幅度从"2.5 倍"回落到 **1.8 倍**（3.19 → 1.77）。 这是一个值得单独讲的实验方法学教训：**比较两条空间填充曲线之前，先确认它们站在同一把尺子上。** 任何对比实验，只要两方的"其他条件"不一致，读出来的差距就不可信。 ### 3.3 落地形态 - **Geohash**：基于 Z 序的字符串编码，把经纬度编成 base32 字符串，相邻前缀表示相邻区域。 - **Google S2**：用六条 Hilbert 曲线首尾相连覆盖球面，层级 level 0–30。 - **Uber H3**：六边形层次网格，分辨率 0–15，总单元数 c(r) = 2 + 120·7^r，且**每层恰有 12 个五边形**（这是六边形铺满球面不可避免的拓扑代价）。 - **PostGIS**：`ST_DWithin` 是索引感知（index-aware）的距离阈值谓词，能走 GiST 索引；`ST_Distance` 不使用索引（PostGIS 官方文档 Chapter 5）。写查询时用 `ST_DWithin(a.geom, b.geom, d)` 而不是 `ST_Distance(...) < d`。 ## 四、实验环境与准备 - Python ≥ 3.9；需要 `numpy`、`matplotlib` - 可选：`ipycanvas`（仓库 notebook 用它画交互式画布，本指导书改用 matplotlib，避免额外依赖） - 数据文件：`notebooks/latitude_longitude.bin`（仓库自带） - 对应 notebook：`notebooks/Z-Order-Curve.ipynb` ```bash python -m pip install numpy matplotlib ``` ## 五、实验步骤 ### 步骤 1：教学版 z 码（字符串形式） 先用最直白的字符串形式把位交织写出来，方便观察每一位从哪来。 ```python def scalar_to_z_code(recursion_level: int, x: int, y: int) -> str: """把 (x, y) 转成 z 码字符串。位序：y_{L−1} x_{L−1} … y_0 x_0""" if recursion_level == 0: return "" x_binary = bin(x)[2:].zfill(recursion_level) y_binary = bin(y)[2:].zfill(recursion_level) assert len(x_binary) == len(y_binary) == recursion_level return "".join(a + b for a, b in zip(y_binary, x_binary)) def z_code_to_scalar(z_code: str) -> tuple[int, int, int]: """z 码字符串 → (递归层级, x, y)，scalar_to_z_code 的逆变换""" assert len(z_code) % 2 == 0 recursion_level = len(z_code) // 2 x = y = 0 exponent = recursion_level - 1 for i in range(0, len(z_code), 2): y += int(z_code[i]) * 2 ** exponent x += int(z_code[i + 1]) * 2 ** exponent exponent -= 1 return recursion_level, x, y print(scalar_to_z_code(3, 2, 1)) # 000110 print(z_code_to_scalar("000110")) # (3, 2, 1) print(z_code_to_scalar(scalar_to_z_code(2, 3, 1))) ``` 请手工核对第一个结果：x = 2 = `010`，y = 1 = `001`，按 `y2 x2 y1 x1 y0 x0` 交织得 `0 0 0 1 1 0`。 ### 步骤 2：打印整张 z 码表 ```python def print_z_codes(recursion_level: int): """按行优先打印整个网格的 z 码（y 从 0 到 2^L−1，x 从 0 到 2^L−1）""" n = 2 ** recursion_level for y in range(n): print(" ".join(scalar_to_z_code(recursion_level, x, y) for x in range(n))) print_z_codes(3) ``` 观察：沿每一行看，z 码的低位在变；跨行时高位在变。整体读起来是一条"Z"字形的路径。 ### 步骤 3：整数版 Morton 码（numpy 向量化） 字符串版只适合教学。真实系统用整数，且要一次算几百万个点，所以写成向量化形式。 ```python import numpy as np def morton(x, y, levels=16): """向量化的 Morton 码。位序与步骤 1 一致：y 位在高，x 位在低。""" x = np.asarray(x, dtype=np.int64) y = np.asarray(y, dtype=np.int64) z = np.zeros_like(x) for i in range(levels): z |= ((y >> i) & 1) << (2 * i + 1) # y 的位放在奇数位 z |= ((x >> i) & 1) << (2 * i) # x 的位放在偶数位 return z print(morton(np.array([2]), np.array([1]), levels=3)) # [6]，对应 "000110" ``` ### 步骤 4：把 Z 曲线画出来（matplotlib，不依赖 ipycanvas） ```python import matplotlib.pyplot as plt def draw_z_curve(L=3): n = 1 << L xs, ys = np.meshgrid(np.arange(n), np.arange(n)) xs, ys = xs.ravel(), ys.ravel() z = morton(xs.astype(np.int64), ys.astype(np.int64), levels=L) order = np.argsort(z, kind="stable") plt.rcParams["font.sans-serif"] = ["Microsoft YaHei", "SimHei", "DejaVu Sans"] fig, ax = plt.subplots(figsize=(6, 6)) ax.plot(xs[order], ys[order], "-o", color="#1F4E79", ms=8, lw=2) ax.set_aspect("equal") ax.set_xlim(-0.8, n - 0.2) ax.set_ylim(-0.8, n - 0.2) ax.set_xlabel("x") ax.set_ylabel("y") ax.set_title(f"Z 曲线（递归层级 {L}，{n}×{n} 网格）") ax.grid(ls=":", alpha=0.4) plt.tight_layout() plt.show() draw_z_curve(3) ``` 把 L 改成 1、2、4 各画一次，看递归结构怎么长出来。 ### 步骤 5：读入真实经纬度数据 数据文件里 float32 交替存放：纬度、经度、纬度、经度…… ```python import numpy as np raw = np.fromfile("notebooks/latitude_longitude.bin", dtype=" 0: rx = 1 if (x & s) > 0 else 0 ry = 1 if (y & s) > 0 else 0 d += s * s * ((3 * rx) ^ ry) if ry == 0: # 旋转象限 if rx == 1: x = s - 1 - x y = s - 1 - y x, y = y, x s //= 2 return d Lh = L # 与 Z 序用同一个网格，比较才公平 xh, yh = x >> (L - Lh), y >> (L - Lh) # Lh == L 时右移 0 位，即原值 h = np.array([xy2d(1 << Lh, int(a), int(b)) for a, b in zip(xh, yh)]) h_order = np.argsort(h, kind="stable") ``` 注意这里**没有**给 Hilbert 单独设一个更粗的网格。两条曲线必须用同一套 $2^{16}$ 网格、同一批查询窗口、同一份点集，否则比出来的倍数是量化层级差出来的，不是曲线本身的性质。代价是这一段要跑几秒（`xy2d` 是逐点的位循环，19.7 万个点各跑 16 轮）。 ### 步骤 8：范围查询候选率实测 ```python rng = np.random.default_rng(0) queries = [] for _ in range(400): cx = rng.integers(int(x.min()), int(x.max())) cy = rng.integers(int(y.min()), int(y.max())) half = rng.integers(maxv // 500, maxv // 100) # 小窗口 queries.append((int(cx - half), int(cy - half), int(cx + half), int(cy + half))) results = {"顺序扫描": [], "Z 序索引": [], "Hilbert 序索引": []} zz = z[z_order] hh = h[h_order] for x0, y0, x1, y1 in queries: in_box = (x >= x0) & (x <= x1) & (y >= y0) & (y <= y1) true_hits = int(in_box.sum()) if true_hits < 5: # 命中太少的查询不参与统计 continue results["顺序扫描"].append(len(x) / true_hits) # Z 序：把 [z(x0,y0), z(x1,y1)] 这段区间作为候选 zlo = morton(np.array([x0]), np.array([y0]), levels=L)[0] zhi = morton(np.array([x1]), np.array([y1]), levels=L)[0] lo = np.searchsorted(zz, zlo) hi = np.searchsorted(zz, zhi, side="right") results["Z 序索引"].append(max(1, hi - lo) / true_hits) # Hilbert 序：同理 hlo = xy2d(1 << Lh, int(x0 >> (L - Lh)), int(y0 >> (L - Lh))) hhi = xy2d(1 << Lh, int(x1 >> (L - Lh)), int(y1 >> (L - Lh))) if hhi < hlo: # Hilbert 码对 x/y 不单调，端点必须兜底 hlo, hhi = hhi, hlo lo = np.searchsorted(hh, hlo) hi2 = np.searchsorted(hh, hhi, side="right") results["Hilbert 序索引"].append(max(1, hi2 - lo) / true_hits) for k, v in results.items(): print(f"{k:14s} 候选率中位数 {np.median(v):.2f}× （参与统计的查询数 {len(v)}）") ``` 本机输出： ``` 顺序扫描 候选率中位数 12.64× Z 序索引 候选率中位数 3.19× Hilbert 序索引 候选率中位数 1.77× ``` ### 步骤 9：画出局部性 ```python fig, axes = plt.subplots(1, 2, figsize=(12.6, 5.4)) axes[0].scatter(lon[::3], lat[::3], s=2, color="#BFBFBF", label="全部点（真实经纬度）") idx = z_order[:4000] axes[0].scatter(lon[idx], lat[idx], s=4, color="#1F4E79", label="Z 序最前的 4000 个点") axes[0].set_xlabel("经度") axes[0].set_ylabel("纬度") axes[0].set_title("Z 序把空间邻近的点聚到了一起") axes[0].legend(markerscale=2) axes[0].grid(ls=":", alpha=0.3) labels = list(results.keys()) vals = [float(np.median(results[k])) for k in labels] axes[1].bar(labels, vals, color=["#B0B0B0", "#1F4E79", "#2E7D52"], width=0.55) axes[1].set_yscale("log") for i, v in enumerate(vals): axes[1].text(i, v, f"{v:.2f}×", ha="center", va="bottom") axes[1].set_ylabel("候选数 / 真实命中数（对数轴，越小越好）") axes[1].set_title("小窗口范围查询的候选率") axes[1].grid(axis="y", ls=":", alpha=0.4) plt.tight_layout() plt.show() ``` ### 步骤 10（选做）：用 notebook 的 ipycanvas 版本 如果装了 `ipycanvas`，可以直接看仓库 notebook 里的交互式画布版本： ```python from ipycanvas import Canvas, hold_canvas def visualize_z_codes(recursion_level: int): n = 2 ** recursion_level results = [(scalar_to_z_code(recursion_level, x, y), x, y) for y in range(n) for x in range(n)] results.sort() size = 300 canvas = Canvas(size=(size, size), sync_image_data=True) mult = size / n _, sx, sy = results.pop(0) with hold_canvas(canvas): for _, nx, ny in results: canvas.begin_path() canvas.move_to((sx + 0.5) * mult, (sy + 0.5) * mult) canvas.line_to((nx + 0.5) * mult, (ny + 0.5) * mult) canvas.stroke() sx, sy = nx, ny return canvas for L_ in range(1, 6): display(visualize_z_codes(L_)) ``` ## 六、观测记录表 z 码手算练习（用步骤 1 的函数验证）： | L | x | y | x 的二进制 | y 的二进制 | z 码 | 逆变换是否还原 | |---|---|---|---|---|---|---| | 3 | 2 | 1 | 010 | 001 | 000110 | | | 3 | 5 | 2 | | | | | | 4 | 7 | 3 | | | | | | （自定） | | | | | | | 网格规模： | 递归层级 L | 网格边长 2^L | 单元数 | z 码位宽 | |---|---|---|---| | 3 | 8 | 64 | 6 | | 8 | | | | | 16 | | | | 真实数据信息（把本机输出填入）： | 项目 | 参照值（本机实测） | 你的实测 | |---|---|---| | 点数 | 197,229 | | | 纬度范围 | 47.28 … 55.07 | | | 经度范围 | −2.01 … 50.51 | | 范围查询候选率（随机小窗口，命中数 ≥ 5 才计入，本机有效窗口 89 个；两种曲线同为 2^16 网格）： | 方案 | 参照值（中位数） | 你的实测 | |---|---|---| | 顺序扫描 | 12.64× | | | Z 序索引（2^16 网格） | 3.19× | | | Hilbert 序索引（2^16 网格） | 1.77× | | 改变条件后再测一次： | 实验条件 | 顺序扫描 | Z 序索引 | Hilbert 序索引 | |---|---|---|---| | 查询框放大 10 倍 | | | | | L 从 16 降到 10 | | | | | 两者都用同一层级 L=16（现在的默认口径） | 12.64× | 3.19× | 1.77× | ## 七、思考题 1. 手算 `scalar_to_z_code(3, 5, 2)`，写出 x、y 的二进制与交织过程，再用步骤 1 的代码验证。 2. 求 `z_code_to_scalar("011001")` 的结果，并说明逆变换为什么能还原出 x 与 y。 3. 递归层级 L = 16 时，网格边长是多少？单元总数是多少？z 码需要几个二进制位、能否放进一个 64 位整数？ 4. 为什么范围查询必须"先取候选再精确过滤"？请举一个具体例子：某个点的 z 值落在 `[z(x0,y0), z(x1,y1)]` 区间内，但点本身不在查询矩形内。 5. 关于 Hilbert 与 Z 序的优劣，文献结论一致吗？请说明 Dai & Su（DMTCS 2003）得到的是什么结论，以及工程上应该怎么处理。 6. 比较两条空间填充曲线时，如果一方用 2¹⁶ 网格、另一方用 2¹⁰ 网格，会出现什么问题？本实验是怎么保证公平的？请再设计一个**可能证伪**"Hilbert 一定优于 Z 序"的实验方案。 7. Z 序在真实系统里有哪些落地形态？请说明 Geohash、Google S2、Uber H3 各自的编码方式，以及 H3 每层有多少个五边形。 8. 在 PostGIS 里想查"距离某点 500 米以内的餐厅"，为什么应该写 `ST_DWithin` 而不是 `ST_Distance < 500`？ ## 八、思考题答案与解析 **第 1 题** ``` x = 5 = 101 y = 2 = 010 位序 y2 x2 y1 x1 y0 x0 0 1 1 0 0 1 z 码 = "011001" ``` 验证：`print(scalar_to_z_code(3, 5, 2))` 输出 `011001`。 **第 2 题** `"011001"` 长度 6 → L = 3。按 2 位一组拆：`01|10|01`，指数从 2 递减到 0。 ``` 第 1 组 "01"（指数 2）：y += 0×4 = 0，x += 1×4 = 4 第 2 组 "10"（指数 1）：y += 1×2 = 2，x += 0×2 = 0 第 3 组 "01"（指数 0）：y += 0×1 = 0，x += 1×1 = 1 结果：x = 5，y = 2 → 返回 (3, 5, 2) ``` 能还原的原因是位交织是**一一映射**：z 码的偶数位（从低位数的第 0、2、4…位）恰好是 x 的各位，奇数位是 y 的各位，把它们按原来的权重 2^i 累加回去即可。 **第 3 题** - 网格边长 = 2^16 = 65,536 - 单元总数 = 2^16 × 2^16 = 2^32 = 4,294,967,296（约 43 亿） - z 码位宽 = 2L = 32 位，**能**放进 64 位整数（也能放进无符号 32 位整数） 如果 L 到 32，z 码就是 64 位，必须注意有符号整数的溢出问题（要用 `uint64`）。 **第 4 题** 取 L = 2（4×4 网格），查询矩形 x ∈ [1, 2]、y ∈ [0, 1]。 ``` z(1, 0)：x=01, y=00 → y1 x1 y0 x0 = 0 0 0 1 = 1 z(2, 1)：x=10, y=01 → y1 x1 y0 x0 = 0 1 1 0 = 6 候选区间 = [1, 6] ``` 区间里的 z = 2 对应二进制 `0010`：`y1 x1 y0 x0 = 0 0 1 0` → y = `01` = 1，x = `00` = 0，即点 **(0, 1)**。它的 x = 0 不在 [1, 2] 内——**落在 z 区间里却不在矩形里**。 原因：z 区间是矩形的一个超集，一维区间无法精确表达二维矩形。所以精化阶段（逐个用真实坐标判断）不可省略。工业实现用 `BIGMIN` 算法把大区间切成若干"真正落在矩形内"的小段，减少这种无用候选。 **第 5 题** 不一致。Dai & Su（DMTCS 2003）在"平均跨簇距离"与"平均总跨簇距离"两项统计上得到的结论是 **Z 序家族优于 Hilbert 家族**；而很多教材和工程经验又说 Hilbert 局部性更好。分歧来自**度量口径**不同：用哪种距离、按哪种查询分布来统计，结论会变。 工程做法：不要照搬结论，用自己的真实数据 + 真实查询分布跑一遍本实验的步骤 8，按候选率（或直接按端到端查询耗时）来选。本实验在**同一 2¹⁶ 网格、同一批 89 个有效窗口**下测得 Hilbert 1.77× 优于 Z 序 3.19×，但这只是"这批数据、这批窗口"上的结论，换数据可能翻转（见第 6 题）。 **第 6 题** **会得出假优势。** 网格越粗，同一个查询框落到码空间里覆盖的区间越短，候选集天然更小。粗网格那一方"赢"下来的部分，来自量化层级差异，不是曲线本身的性质。 本案例库早期版本就犯过这个错：Z 码用 2¹⁶ 网格、Hilbert 码用 2¹⁰ 网格，测出 Hilbert 1.29×，看上去优势很大。把两者统一到 **同一 2¹⁶ 网格**、同一批 197,229 个点、同一组 89 个有效查询窗口之后，Hilbert 的候选率变成 **1.77×**（Z 序 3.19×、顺序扫描 12.64×）——领先幅度明显收窄，但**这个数才是可复现的**。1.29× 已作废，不要再引用。 **保证公平的三条**（缺一条，比较就没有意义）： 1. 网格层级相同（本实验 `Lh = L = 16`）； 2. 查询窗口集合相同（同一批 400 个随机窗口，同样按"真实命中 ≥ 5"筛出 89 个）； 3. 点集相同且排序方式一致（都用 `argsort(kind="stable")`）。 **可能证伪"Hilbert 一定优于 Z 序"的实验设计**（任选其一或组合）： ```python # 思路：把"窗口大小"和"点集分布"作为两个变量扫一遍，看结论是否翻转 for half in (maxv // 1000, maxv // 500, maxv // 100, maxv // 20): # 窗口从很小到很大 queries = make_queries(half) # 同一 seed，保证四条曲线面对同一批窗口 print(half, med_ratio(z_order), med_ratio(h_order)) ``` - **扫窗口大小**：窗口极小时两者候选率都接近 1，差距被抹平；窗口大到覆盖多个"簇"时，Z 序的大跳变反而可能更少。若在某个窗口区间出现 Z 序优于 Hilbert，就证伪了。 - **换点集分布**：把真实经纬度换成均匀分布 / 强聚簇分布 / 沿一条线分布，观察排序是否改变。曲线局部性的好坏本来就依赖数据分布。 - **换度量**：除了候选率，再统计"候选区间的段数"（Z 序一个大区间往往要切成多段）与端到端耗时。候选率低但段数多，实际不一定更快。 - **换维度**：把 2 维换成 3 维（见思考题里 Morton 的三维扩展），高维下 Hilbert 的优势通常会缩小。 做完这四组，你大概率会得到一个比教材结论更谨慎的说法：**在真实数据、小窗口范围查询这个具体场景下，Hilbert 优于 Z 序；但这个优势不是普遍规律。** **第 7 题** - **Geohash**：基于 Z 序，把经纬度交替二分后编成 base32 字符串；前缀相同表示区域相邻（常用作前缀匹配的粗过滤）。 - **Google S2**：把球面（立方体投影后的六个面）用 **六条 Hilbert 曲线首尾相连** 覆盖，层级 level 0–30，单元 ID 是 64 位整数。 - **Uber H3**：六边形层次网格，分辨率 0–15，总单元数 c(r) = 2 + 120·7^r；由于六边形无法完美铺满球面，**每层恰有 12 个五边形**（这是拓扑上不可避免的）。 三者的取舍：Geohash/S2 的四边形单元便于位运算与前缀查询；H3 的六边形单元到邻居的距离更均匀（每个单元有 6 个等距邻居），适合做邻域聚合，但代价是那 12 个五边形。 **第 8 题** `ST_DWithin(geom1, geom2, d)` 是**索引感知**（index-aware）的谓词：PostGIS 会先用 GiST 索引的包围盒做粗过滤，再算精确距离，因此能走索引。 而 `ST_Distance(geom1, geom2) < d` 是一个普通的表达式：优化器不知道它等价于"距离阈值"，只能对每一行都算一次距离再做比较，**不使用索引**，在大表上会退化成全表扫描。 写法对照： ```sql -- 推荐：能走 GiST 索引 SELECT * FROM restaurants WHERE ST_DWithin(geom, ST_MakePoint(:lon, :lat)::geography, 500); -- 不推荐：全表扫描 SELECT * FROM restaurants WHERE ST_Distance(geom, ST_MakePoint(:lon, :lat)::geography) < 500; ``` （依据：PostGIS 官方文档 Chapter 5。） ## 九、扩展任务（选做） 1. **位并行的 Hilbert 码**：用查表法（magic bits）把 `xy2d` 改成 O(1) 的位运算版本，再把两者都放到 2^16 网格上公平对比。 2. **BIGMIN / LITMAX**：实现真正的 Z 序范围查询算法——给定当前 z 值与查询矩形，直接算出"下一个落在矩形内的最小 z 值"，跳过大段空洞。统计它比"单个区间"少取多少候选。 3. **多段区间切分**：不实现 BIGMIN 也行，可以把矩形按 z 码前缀切成若干子矩形（decomposition），每个子矩形对应一个连续区间，统计总候选数。 4. **三维扩展**：把位交织推广到三维（x, y, t），思考时间维度加入后局部性如何变化。 5. **Geohash 编解码**：写一个 base32 的 Geohash 编解码器，验证"前缀相同 = 区域相邻"，并找出反例（边界处前缀不同却相邻）。 6. **真实引擎对照**：如果有 PostGIS 环境，把这 197,229 个点导入，建 GiST 索引，比较 `ST_DWithin` 与 `ST_Distance < d` 的执行计划与耗时。 ## 十、常见错误与排查 | 现象 | 原因 | 处理 | |---|---|---| | `z_code_to_scalar` 还原出的 x、y 与输入不符 | 位序搞反了（写成 x 在高位、y 在低位） | 统一约定：y 位在奇数位、x 位在偶数位，与 `scalar_to_z_code` 配套 | | `morton()` 结果与字符串版不一致 | 循环里的移位写成 `(2*i)` 与 `(2*i+1)` 时放反了 | y 用 `<< (2*i+1)`，x 用 `<< (2*i)` | | 读入的经纬度范围离谱（如几百上千） | 数据读取时 lat/lon 取反了，或字节序不对 | 用 `dtype="