实验 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 曲线的关键性质（也是它的缺点）：相邻的两个 z 值在空间上几乎总是相邻的（局部性好），但在象限边界上会出现**大跳跃**（例如从 (3,0) 跳到 (0,1)）。

3.2 用真实数据测局部性

实验数据来自仓库自带的 `notebooks/latitude_longitude.bin`：

• 197,229 个点

• 纬度范围 47.28° – 55.07°，经度范围 −2.01° – 50.51°

• 存储格式：float32 交替存放（纬度、经度、纬度、经度……），小端序

实验做法：随机生成 400 个小窗口矩形查询框，只统计真实命中数 ≥ 5 的查询（命中太少的查询不足以体现差异），记录每种方案的"候选数 / 真实命中数"的中位数。


本机实测结果（中位数）：

| 方案 | 候选数 / 真实命中数 |
|---|---|
| 顺序扫描 | 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="<f4")
lat = raw[0::2].astype(np.float64) # 偶数下标 = 纬度
lon = raw[1::2].astype(np.float64) # 奇数下标 = 经度

print(f"点数: {len(lat):,}")
print(f"纬度: {lat.min():.2f} … {lat.max():.2f}")
print(f"经度: {lon.min():.2f} … {lon.max():.2f}")
```

本机输出：

```
点数: 197,229
纬度: 47.28 … 55.07
经度: -2.01 … 50.51
```

步骤 6：量化并排序

```python
L = 16
maxv = (1 << L) - 1
x = np.clip(((lon + 180.0) / 360.0 * maxv).astype(np.int64), 0, maxv)
y = np.clip(((lat + 90.0) / 180.0 * maxv).astype(np.int64), 0, maxv)

z = morton(x, y, levels=L)
z_order = np.argsort(z, kind="stable") # 按 z 码排序后的下标
print("前 5 个 z 码:", z[z_order][:5])
```

步骤 7：Hilbert 码（用于对照）

```python
def xy2d(n, x, y):
"""Hilbert 码：坐标 → d（n 为网格边长）。

与 Wikipedia 标准一致：rot 收到的是当前轮边长 s（xy2d 内调用 rot(s,...)）。
翻转后坐标可能为负，靠补码语义保证后续轮的位判断仍正确。"""
d = 0
s = n // 2
while s > 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="<f4"`（小端 float32），`lat = raw[0::2]`、`lon = raw[1::2]` |
| 读入的点数只有一半 | 把整个数组当成 (lat, lon) 对，没做奇偶拆分 | 总长度应是点数的 2 倍；打印 `len(raw)` 与 `len(lat)` 检查 |
| `np.fromfile` 报文件不存在 | 运行目录不在仓库根 | 用绝对路径，或先 `os.chdir` 到包含 `notebooks/` 的目录 |
| 量化后所有点挤在一个小角落 | 量化公式用了数据的实际范围但没有归一化到 [0, 2^L−1] | 用 `clip` 保证结果落在 [0, maxv]，或用实际 min/max 做归一化 |
| Hilbert 那一步特别慢 | `xy2d` 是 Python 逐点循环，L=16 时 20 万点很慢 | 先用较粗层级（L=10）或抽样；想快就做扩展任务 1 |
| Hilbert 候选率异常偏低（接近 1） | 忘了做 `hhi < hlo` 的端点交换，负区间被 `max(1, …)` 兜成 1 | Hilbert 码对 x/y 不是单调的，`searchsorted` 前先做端点交换 |
| 候选率出现 1.0 或 0 | 查询框太大/太小，或 `searchsorted` 的 `side` 参数用错 | 区间右端用 `side="right"`；候选数用 `max(1, hi-lo)` 兜底 |
| 顺序扫描的候选率看起来"太好了" | 查询框很大，命中数本来就多 | 缩小窗口（`half` 取 `maxv//500 … maxv//100`），并只统计命中数 ≥ 5 的查询 |
| 图中中文变方框 | 缺中文字体 | 设置 `plt.rcParams["font.sans-serif"]`，或改用英文标签 |
| 导入 `ipycanvas` 失败 | 未安装或不在 Jupyter 环境 | 用步骤 4 的 matplotlib 版本替代 |