蒲公英空間擴散與入侵物種影響評分

封存文章

本文是 HiMCM 2023 Problem A 的教學重構。舊版中的氣候倍率、物種分數和權重都是示範假設，不能當作蒲公英、葛藤或鳳眼藍的實證風險排名。本版本保留演算法概念，同時指出校準與地區定義所需資料。

蒲公英（Taraxacum officinale）的瘦果連着冠毛，可以由風搬運。題目要求由一株位於一公頃空地旁的蒲公英出發，預測 1、2、3、6、12 個月的擴散，比較溫帶、乾旱與熱帶條件，再建立可套用於其他入侵植物的 impact factor。

這裏其實有兩個不同問題：

動態問題： 種子何時產生、落在何處、能否萌發並繁殖？
決策問題： 一個物種在指定地區造成多大生態、經濟與管理負擔？

第一個問題需要空間與時間；第二個問題需要資料來源、地區語境與可解釋的權重。

一、由 Logistic 到空間模型

總量 Logistic 方程

\frac{dP}{dt}=rP\left(1-\frac{P}{K}\right)

可描述密度制約，但看不到風向、邊界、空地形狀或散播距離。它適合作為總量 benchmark，而不是題目主模型。

把一公頃正方形表示成

\Omega=[0,100]\times[0,100]\ {\rm m^2},

並切成 cell。每個 cell $i$ 可有成年株數 $A_i(t)$ 、seed bank $B_i(t)$ ，以及資源狀態 $R_i(t)$ 。一個月的更新可拆成五步。

1. 種子產生

對每株成熟植物，種子數可用 over-dispersed distribution：

S_i(t)\sim\operatorname{NegBin} \bigl(A_i(t)\,q(T_t,M_t),\ \theta\bigr),

其中 $q$ 連接溫度、土壤水分與 phenology， $\theta$ 控制個體差異。這比固定「spread probability」更容易對接每株花序、每花序瘦果數和成熟月份資料。

2. 風媒散播核

由 cell $i$ 到 $j$ 的概率可寫成方向性 kernel：

K_{ij}(t) \propto \exp\left[ -\frac12 (\Delta x_{ij}-\mu_t)^\mathsf T \Sigma_t^{-1} (\Delta x_{ij}-\mu_t) \right],

$\mu_t$ 沿平均風向， $\Sigma_t$ 表示順風與側風擴散。若要捕捉罕見長距離事件，可用 Gaussian core 與 fat-tailed component 的 mixture：

K=(1-\omega)K_{\rm local}+\omega K_{\rm long}.

每批種子依 multinomial distribution 分配至目的地。邊界外的種子應記錄為離開研究地塊，而不是從另一邊重新出現。

3. 萌發與建立

落地不等於成株。對 cell $j$ ：

G_j(t)\sim \operatorname{Binomial} \left[ B_j(t)+D_j(t),\ g(T_t,M_t,L_j,N_j) \right],

其中 $D_j$ 是新到達種子， $L_j,N_j$ 是光與養分指標。建立概率可再乘上密度制約：

p_{{\rm est},j} =p_0 f_T(T_t)f_M(M_t) \left(1-\frac{A_j}{K_j}\right)_+.

4. 存活與成熟

新苗、juvenile 與 adult 最好分開，形成 stage-structured model：

\mathbf n_j(t+1) =\mathbf M_j(t)\mathbf n_j(t) +\mathbf I_j(t),

而不是讓新佔據 cell 下一個 time step 立刻產生同等數量種子。若省略生命史延遲，擴散速度通常會被嚴重高估。

5. 競爭

若要研究另一物種 $Q$ ，可在 cell level 使用競爭係數：

\begin{aligned} P_{t+1} &=P_t+r_PP_t \left(1-\frac{P_t+\alpha_{PQ}Q_t}{K_P}\right),\\ Q_{t+1} &=Q_t+r_QQ_t \left(1-\frac{Q_t+\alpha_{QP}P_t}{K_Q}\right). \end{aligned}

但 $\alpha_{PQ}$ 不能靠物種名稱猜測；它需要共同園圃實驗、observational co-occurrence 或文獻先驗。陽光、水和養分亦未必可由單一競爭係數充分表示。

二、氣候不是三個任意倍率

舊程式將基本散播概率在乾旱氣候乘 0.5、熱帶乘 1.5。這種分類方便畫圖，卻混合了不同機制：

風速與湍流影響散播距離；
降雨與土壤水分影響萌發及存活；
溫度與日長影響物候及繁殖時點；
地被與擾動影響建立機會。

更合理的情景表會直接指定月度 forcing：

情景	月溫度	土壤水分	風速／風向	萌發季節
溫帶	有明顯季節	中等	觀測分布	春、秋較高
乾旱	日夜差較大	低且間歇	可有強陣風	降雨脈衝後
濕熱	長期偏暖	高	季風／對流	受熱與淹水限制

參數要由指定地區資料估計，因為「熱帶」不等於所有生命階段都更適合蒲公英。

三、可重現的模擬骨架

以下程式只展示邊界與方向散播的基本設計：

import numpy as np

def disperse(seed_counts, wind_vector, rng, sigma_parallel=8, sigma_cross=3):
    arrivals = np.zeros_like(seed_counts, dtype=int)
    escaped = 0
    rows, cols = seed_counts.shape

    direction = np.asarray(wind_vector, dtype=float)
    speed = np.linalg.norm(direction)
    unit = direction / speed if speed > 0 else np.array([1.0, 0.0])
    cross = np.array([-unit[1], unit[0]])

    for i, j in np.argwhere(seed_counts > 0):
        count = int(seed_counts[i, j])
        downwind = rng.normal(loc=2 * speed, scale=sigma_parallel, size=count)
        lateral = rng.normal(loc=0, scale=sigma_cross, size=count)
        offsets = downwind[:, None] * unit + lateral[:, None] * cross
        destinations = np.rint(offsets + [i, j]).astype(int)

        inside = (
            (destinations[:, 0] >= 0)
            & (destinations[:, 0] < rows)
            & (destinations[:, 1] >= 0)
            & (destinations[:, 1] < cols)
        )
        escaped += np.count_nonzero(~inside)
        for di, dj in destinations[inside]:
            arrivals[di, dj] += 1

    return arrivals, escaped

這裏沒有使用 modulo 邊界。舊版

(i + di + wind_x) % size

令右邊飛走的種子從左邊回來，等於把一公頃地塊變成 torus，會低估外流並扭曲前緣速度。

完整 monthly update 應明確處理 seed production、dispersal、germination、 survival 和 maturation：

def one_month(state, climate, rng):
    seeds = produce_seeds(state.adults, climate, rng)
    arrivals, escaped = disperse(seeds, climate.wind, rng)
    germinated = germinate(state.seed_bank + arrivals, climate, rng)
    seedlings = establish(germinated, state.resources, rng)
    next_state = survive_and_age(state, seedlings, climate, rng)
    return next_state, escaped

每個函數都應有單位、參數來源與測試。

四、輸出與不確定性

題目要求 1、2、3、6、12 個月，所以 simulation step 必須清楚等於一個月。舊程式跑 30 steps，圖軸卻以 steps / 12 標成年，動畫又叫它月份；這使時間尺度不一致。

每個 checkpoint 應至少報告：

成株總數；
佔據面積與 occupancy probability map；
由初始株到前緣的距離；
地塊外流種子數；
種群滅絕概率；
不同氣候情景間的差。

checkpoints = np.array([1, 2, 3, 6, 12])
counts = np.stack([
    simulate(seed=seed, months=12).adult_counts
    for seed in range(500)
])

median = np.median(counts[:, checkpoints], axis=0)
lower = np.quantile(counts[:, checkpoints], 0.025, axis=0)
upper = np.quantile(counts[:, checkpoints], 0.975, axis=0)

lower 至 upper 是 Monte Carlo outcome 的 95% prediction interval，不是估計平均值的 confidence interval。舊文只跑 10 次便稱作 95% CI，樣本太少而且術語不準確。

原模型報告 30 steps 後平均 2,856 個佔據 cell、覆蓋 28.56%。由於一個 cell 被當成一株、繁殖沒有成熟延遲、參數未校準且採 toroidal boundary，該數字不能轉換成真實一公頃內的蒲公英數目。

五、Impact factor 的設計

舊版選七個因素：

growth rate；
reproductive capacity；
dispersal ability；
competitive ability；
environmental tolerance；
ecosystem impact；
economic impact。

加權分數為

I=\sum_{k=1}^{7}w_kz_k,\qquad \sum_kw_k=1,\quad0\le z_k\le1.

形式簡單，但要先處理三個問題。

1. 地區與受體

一個物種是否「入侵」取決於地區、引入歷史與受影響生態系。葛藤在美國東南部與鳳眼藍在某些淡水系統的風險，不能與某地原生族群直接混成全球分數。每次評估必須列明：

assessment region；
habitat／ecosystem；
time horizon；
ecological receptors；
management objective。

2. 證據表

每個 $z_k$ 應由觀測或文獻指標映射，例如：

構面	可觀測指標	方向
擴散壓力	年擴張距離、propagule pressure	越高越嚴重
建立能力	不同 habitat 的 establishment probability	越高越嚴重
生物多樣性	native richness／cover 改變	損失越大越嚴重
生態功能	水流、養分循環、火風險改變	偏離越大越嚴重
經濟損失	管理成本、產量或基建損失	越高越嚴重
可逆性	清除後恢復時間	越慢越嚴重

缺資料不應用一個看似精確的 0–10 主觀均值掩蓋，而應標記 evidence gap。

3. 權重與非補償性

純加權和允許極高生態破壞被低經濟成本抵消。若某些 criterion 是不可接受門檻，可用 veto rule 或分開報告 ecological severity 與 management burden。權重應由 stakeholder elicitation、政策文件或 sensitivity analysis 支持：

I(\mathbf w,\mathbf z) =\mathbf w^\mathsf T\mathbf z.

檢查在合理權重集合 $\mathcal W$ 中，物種排序是否穩定：

\Pr_{\mathbf w\in\mathcal W} \left[I_a>I_b\right].

六、Monte Carlo 不會自動令分數客觀

舊版把主觀設定的均值與標準差放入 truncated normal，再抽樣 10,000 次，得到：

物種	舊版平均分數
蒲公英	0.658
葛藤	0.788
鳳眼藍	0.797

大量抽樣只會精確傳播輸入假設；它不能驗證均值 7、8、9 或權重 0.15 是否正確。此外，mean ± 1.96 × sample SD 描述的是假設分布的近似範圍，不是 Monte Carlo mean 的 confidence interval。

因此這三個數字應標成illustrative score under assumed inputs，不能說明鳳眼藍在所有地區必然比葛藤更有害。

七、模型驗證

一個更可信的驗證流程包括：

用已知釋放點與後續調查地圖校準散播 kernel；
以不同年份或不同地塊作 out-of-sample validation；
比較 occupancy、front speed 和 density，而非只比較總數；
對 germination、survival、fecundity 做 global sensitivity；
將 impact score 與獨立 expert assessment 或既有風險框架比較；
報告錯誤分類的管理後果。

若沒有這些資料，模型仍可作 scenario generator，答案應寫成「在假設 X 下會出現 Y」，而不是「模型預測真實地塊必然有 Y 株」。

八、管理用途

空間模型與 impact score 可組合成一個優先次序工具：

\text{Priority}_s = \Pr(\text{establishment}_s) \times \mathbb E[\text{impact}_s] \times \text{exposure}_s.

但實際決策還要加入偵測概率、清除成本與延遲造成的損失。例如 early detection 模型可比較：

邊界巡查與全面網格巡查；
每月與每季監測；
優先移除 seed-producing adults；
在 prevailing downwind direction 增加樣帶。

這比單獨製作一個沒有地區與時間的「入侵分數」更能支持行動。

結論

蒲公英問題展示了由簡單總量模型走向空間、生態與決策模型的過程：

Logistic growth 只能提供總量基準；
風媒散播需要方向性 kernel、開放邊界和長距離尾部；
種子到成株之間要有萌發、建立、存活與成熟階段；
氣候應分解為溫度、水分、風與物候，而非任意倍率；
impact factor 必須指定地區、證據、權重與不可補償門檻；
Monte Carlo 傳播不確定性，但不會令主觀輸入變成實證。

因此，原有 2,856 株和三個 0–1 分數的真正作用，是提醒我們下一步要收集甚麼資料，而不是提供完成的生態風險結論。

參考資料

COMAP, Inc. (2023). HiMCM 2023 Problem A and accompanying problem materials.
Nathan, R., Katul, G. G., Bohrer, G., et al. (2011). Mechanistic models of seed dispersal by wind. Theoretical Ecology, 4, 113–132.
Soons, M. B., & Ozinga, W. A. (2005). How important is long-distance seed dispersal for the regional survival of plant species? Diversity and Distributions, 11(2), 165–172.
Pyšek, P., Hulme, P. E., Simberloff, D., et al. (2020). Scientists’ warning on invasive alien species. Biological Reviews, 95(6), 1511–1534.
Blackburn, T. M., Pyšek, P., Bacher, S., et al. (2011). A proposed unified framework for biological invasions. Trends in Ecology & Evolution, 26(7), 333–339.