這個項目的構思從一開始是如此的，譬如我獲取到城市的 44 個地點 × 7 個指標 × 8 個年份的城市數據，我想要判斷「高建築密度的點是不是傾向跟其他高密度點聚在一起」，後來我使用了空間自相關概念 —— Moran's I。

一、空間自相關

假設你手上有城市 44 個採樣點的建築密度數據。如果你只看「平均值 + 標準差」，你只知道整體分佈的形狀，卻無法觀察到空間結構信息：

  • 隨機分佈：高密度點到處散落，沒有規律
  • 聚集分佈：高密度的點都擠在某幾個區塊（譬如市中心 + 副都心）
  • 棋盤分佈：高低密度交錯，鄰居總是相反

這三種情況的「均值 + 標準差」可以完全一樣，但對都市規劃的意義天差地遠。而 Moran's I 是把「空間結構」壓成一個數字的方法。

二、Moran's I 的判斷依據來源

Moran's I 在每一個數據點上判斷 「鄰居的值跟自己的值，方向上一致嗎？」，如果是，就是正相關。

具體地如下：

  • 我的建築密度比平均高 → 鄰居也比平均高 → 正貢獻
  • 我比平均高 → 鄰居比平均低 → 負貢獻
  • 把所有點對都這樣算，加起來標準化 → Moran's I

結果區間 ≈ [-1, 1]：

  • I > 0 顯著 → 正空間自相關（聚集）
  • I < 0 顯著 → 負空間自相關（棋盤）
  • I ≈ 0 → 隨機分佈

期望值不是 0 而是 E[I] = -1/(n-1)，n=44 時約等於 -0.023，誤差很小可以忽略。

三、公式

        n     Σᵢⱼ wᵢⱼ (xᵢ - x̄)(xⱼ - x̄)
I  =  ─── × ──────────────────────────
       W           Σᵢ (xᵢ - x̄)²
  • n：地點數
  • wᵢⱼ：地點 i 到地點 j 的「空間權重」
  • W = Σᵢⱼ wᵢⱼ：所有權重之和
  • xᵢ：地點 i 的觀測值
  • x̄：所有觀測值的平均

兩個關鍵零件：

  1. 空間權重矩陣 W —— 怎麼定義「鄰居」？
  2. 協方差項 (xᵢ - x̄)(xⱼ - x̄) —— 鄰居 i 和 j 的值是不是同方向偏離平均？

四、空間權重矩陣怎麼建

city-ana 的做法：逆距離 + 距離帶截斷 + 行標準化。

def build_weight_matrix(lons, lats, threshold=0.015, row_standardize=True):
    coords = np.column_stack([lons, lats])
    n = len(coords)

    # 廣播計算成對距離
    diff = coords[:, np.newaxis, :] - coords[np.newaxis, :, :]
    dist = np.sqrt((diff ** 2).sum(axis=2))
    np.fill_diagonal(dist, np.inf)  # 對角線設 inf，避免自己當自己鄰居

    # 距離帶 + 逆距離
    with np.errstate(divide="ignore"):
        W = np.where(dist <= threshold, 1.0 / dist, 0.0)

    if row_standardize:
        row_sums = W.sum(axis=1, keepdims=True)
        row_sums[row_sums == 0] = 1.0
        W = W / row_sums

    return W, dist

四個設計：

1. threshold = 0.015 (≈1500m)

採樣點間距大約 500m。閾值 1500m 涵蓋每個點的「1~3 階鄰居」（直接鄰居 + 隔一個的鄰居）。

太小（譬如 500m）：很多點變成孤島，W 大量為 0，結果不穩。
太大（譬如 5000m）：「鄰居」太多，等於把整個城市當鄰居，自相關訊號被稀釋。

建議經驗法則：閾值 ≈ 採樣間距的 2~3 倍。

2. 逆距離 1/d

常用的權重設計會考慮近的鄰居權重大，遠的權重值小，除此之外還有k-最近鄰「Queen contiguity」等）。

3. 距離帶截斷

當超過閾值的則權重設為 0，有點類似 Relu 在 0 值則直接砍掉。這讓 W 變稀疏，計算快也避免遠端噪音。

4. 行標準化

W[i, :] / sum(W[i, :])，讓每一行加起來等於 1。這意味著「每個點的鄰居加權平均」是個真正的平均，而不是「鄰居越多分數越高」。

行標準化後 Moran's I 的解讀變得很直觀：I 接近 1 = 鄰居加權平均跟自己幾乎一致。

五、根据 I 值判斷顯著

計算出 I = 0.3 之後我們需要思考這是真的聚集特性還是只是隨機波動，因此要做統計顯著性檢定。下面有兩種方法：

方法 A：常態近似 + Z 檢定

在「值是隨機分佈」的虛無假設 H₀ 下，計算 I 的期望值 E[I] 和方差 Var[I]，算 Z 分數和 p 值：

z = (I - E[I]) / sqrt(Var[I])
p = 2 × (1 - Φ(|z|))   # 雙尾

E[I] 的公式很簡單：-1/(n-1)
Var[I] 的公式很長：

S1 = 0.5 * ((W + W.T) ** 2).sum()
S2 = ((W.sum(axis=1) + W.sum(axis=0)) ** 2).sum()

m2 = float((z ** 2).mean())
m4 = float((z ** 4).mean())
b2 = m4 / (m2 ** 2)  # 樣本峰度

A = n * ((n**2 - 3*n + 3) * S1 - n * S2 + 3 * W_sum**2)
B = b2 * ((n**2 - n) * S1 - 2*n * S2 + 6 * W_sum**2)
C = (n - 1) * (n - 2) * (n - 3) * (W_sum**2)

Var_I = (A - B) / C - E_I**2

這是 Moran (1950) 原始論文裡推導好的封閉解。

優點：快，O(n²) 一次算完。
缺點：依賴常態近似，當 n 很小（<20）時不可靠。

方法 B：置換檢定（Permutation test）

把觀測值隨機洗牌 N 次（n_perm=999 通常足夠），每次算一遍 I，看實際 I 落在分佈的哪個尾巴：

if permutation:
    rng = np.random.default_rng(42)
    perm_I = np.zeros(n_perm)
    for i in range(n_perm):
        xp = rng.permutation(x)
        zp = xp - xp.mean()
        num_p = float((W * np.outer(zp, zp)).sum())
        den_p = float((zp ** 2).sum())
        perm_I[i] = (n / W_sum) * (num_p / den_p) if den_p > 0 else 0.0
    p_value = float(((perm_I >= I).sum() + (perm_I <= -I).sum() + 1) / (n_perm + 1))

優點：不依賴常態假設，小樣本更可靠。
缺點：慢，要算 n_perm 次。對 n=44 還好，n=10000 就要等。

規則：n > 30 用常態近似；n < 30 或結果在邊緣（p ≈ 0.05）時用 permutation 復核。

六、批次跑全部年份 × 全部指標

def run_spatial_autocorr(panel, indicators=None, threshold=0.015, permutation=False):
    if indicators is None:
        exclude = {"year", "lon", "lat"}
        indicators = [c for c in panel.select_dtypes(include="number").columns
                      if c not in exclude and not c.startswith("d_")]

    years = sorted(panel["year"].unique())
    results = []

    for year in years:
        df_year = panel[panel["year"] == year].dropna(subset=["lon", "lat"])
        if len(df_year) < 4:
            continue

        lons, lats = df_year["lon"].values, df_year["lat"].values
        W, _ = build_weight_matrix(lons, lats, threshold=threshold)

        for ind in indicators:
            vals = df_year[ind].dropna().values
            if len(vals) < 4 or vals.std() == 0:
                continue
            mi = morans_i(vals, W, permutation=permutation)
            results.append({"year": year, "indicator": ind, "n_points": len(vals), **mi})

    return pd.DataFrame(results)

執行：

python analysis/spatial_autocorr.py

輸出 output/spatial_autocorr.csv：

year,indicator,n_points,I,E_I,Var_I,z_score,p_value,interpretation
2018,building_count,44,0.412,-0.023,0.004,6.875,0.0000,clustered (high-high / low-low)
2018,edge_density,44,0.198,-0.023,0.004,3.495,0.0005,clustered (high-high / low-low)
2018,poi_diversity,44,0.087,-0.023,0.004,1.738,0.0822,random (not significant)
...
2025,building_count,44,0.523,-0.023,0.005,7.728,0.0000,clustered (high-high / low-low)

七、結果解讀

解讀 1：建築密度有強空間聚集，且隨年份增強

  • 2018: I = 0.412, p < 0.001
  • 2025: I = 0.523, p < 0.001

台中市的建築密度高的地點確實傾向擠在一起，而且 8 年下來這個聚集程度還在增強。市中心愈來愈密集，郊區跟得慢，差距拉大。

這個結論可以作爲強聚落特性的證據，因為：

  • p < 0.001 是強顯著
  • I = 0.5+ 是強烈空間相關
  • 趨勢（隨年份增強）跨年份都成立

解讀 2：POI 多樣性早期不顯著，後期顯著

  • 2018: I = 0.087, p = 0.08（不顯著）
  • 2025: I = 0.245, p = 0.002（強顯著）

早期台中的「功能混合」是隨機分佈的，後期出現空間聚集。意思是：早期市中心和郊區都還算單一功能，後期市中心發展出明顯的「複合機能商圈」、郊區仍以單一住宅為主。

解讀 3：路網長度不顯著

當大部分點路網長度為 0，方差為 0，I 算不出來。這就是 README TODO 提到的「路網查詢標籤需要調整」的副作用。

八、Moran's I 之後可以做什麼

在使用全局 Moran's I 得知整體空間結構特性後，但我們仍舊無法得知在具體的點是熱點/冷點的差異，因此下一步是局部 Moran's I (LISA, Local Indicators of Spatial Association)：

I_i = (xᵢ - x̄) × Σⱼ wᵢⱼ (xⱼ - x̄) / m₂

每個地點都有自己的 I_i，可以畫出「熱點 (HH)、冷點 (LL)、異常 (HL/LH)」的空間分佈圖。再次我們就不做延伸探討。

無標籤

關注作者：

新增評論