熱帶氣旋 > 基於統計模型之熱帶氣旋信號機率 - 理論篇 [更新] [English]

熱帶氣旋信號機率計算 - 理論篇

本文詳述我們用作擬合的統計模型及計算熱帶氣旋信號機率的過程及方法。注意：此頁載有數學符號，需使用 Javascript 以顯示此等內容。

請按此顯示本產品，或按此了解使用方法。

1. 有關迴歸分析

熱帶氣旋越接近香港，發出熱帶氣旋信號的機會越高。而一般來說，較強的熱帶氣旋破壞力較大，影響範圍普遍較廣，因此發出熱帶氣旋信號的機會亦較高。是次分析基於統計學的迴歸法，此法嘗試建立反應變數 (因變數，response variable) 及解釋變數 (自變數，explanatory variable(s)) 之間的關係。在此應用中，反應變數為信號是否發出的二元變數，而解釋變數可以多於一個，如熱帶氣旋的地理位置 (座標) 及強度。此關係可寫成

\begin{equation} (時間 t 之信號狀態) \sim (時間 t 之氣旋座標) + (時間 t 之氣旋強度)。 \label{eq:regeqn1} \end{equation}

上式中波浪號 ($\sim$) 可解作「相關因素為」，而加號則解作「及」。對於此類反應變數是二元變數的問題，一般是應用邏輯迴歸法。亦即是說，我們以信號發出的機率而不是數值只能為零或一的二元變數為反應變數。因此，基於熱帶氣旋的座標及強度，模型 \eqref{eq:regeqn1} 能協助我們估計信號發出的機會，而這關係亦符合常理。

2. 歷史信號及強度的重要性

模型 \eqref{eq:regeqn1} 假設熱帶氣旋現時狀況以外的因素與信號發出機會無關。由於氣旋過往動態及趨勢往往影響它對香港的威脅，此假設並不實際。於熱帶氣旋來襲時，信號一般會於一段時間內生效，期間不會多次發出及取消，因此過往跟現時信號狀態間有極大關係。與此同時，氣旋的過往路徑包含其移動方向及速度等資料，我們能以這些資料決定氣旋是正在靠近或遠離本港。基於這些因素，以下的關係較為合理：

\begin{eqnarray} (時間 t 之信號狀態) & \sim & (時間 t 之氣旋座標) + (時間 t 之氣旋強度) + (時間 t-k 之氣旋座標) + (時間 t-k 之信號狀態)。 \label{eq:regeqn2} \end{eqnarray}

上式中 $k$ 為時間間隔。在一、三及八號信號的分析中，時間間隔被定為 6 小時，而在九及十號信號的分析則為 2 小時。對於較低信號，6 小時的長度不至於因間隔太長而失去信號轉換的數據，亦不會因間隔太短而導致數據集中有太多信號維持的情況。由於高信號生效時間較短，使用 6 小時的間隔或致一些信號發出數據被忽略。

留意模型 \eqref{eq:regeqn2} 可視為對於信號狀態的時間序列。此模型只假設於給定各自 $t-k$ 信號狀態時，不同時間的信號狀態為獨立事件，理論上較為有效。

3. 利用懲罰樣條令模型更有彈性

模型 \eqref{eq:regeqn1} 及 \eqref{eq:regeqn2} 皆為相加性模型。亦即是說，如其他變數數值維持，某一解釋變數轉變幅度對應反應變數轉變特定幅度。例子有熱帶氣旋於 (20°N, 118°E) 處的信號機率跟 (20°N, 115°E) 處的信號機率之差異，與熱帶氣旋於 (20°N, 115°E) 處的信號機率跟 (20°N, 112°E) 處的信號機率之差異相同1。但現實中信號機率與氣旋座標存在著非線性的關係，而樣條可協助我們處理此情況。

在進行樣條迴歸中，我們需指定一些節點，而節點的數目與所得出模型的彈性有關：節點越多，模型的彈性越大。透過非線性處理，樣條能更準確地吻合數值間的趨勢；我們應用二維及三維樣條以達到所需彈性。二維樣條描述經緯度轉變跟信號機率轉變的關係，可同時兼顧兩座標系中的交互作用 (即緯度改變所引致的機率改變取決於氣旋當時的經度，反之亦然)，而三維樣條則同時描述經緯度及強度對信號機率的影響。

雖然使用樣條可令模型更有彈性，但由於這些模型有較多參數，有機會出現過度擬合的情況，即擬合得出的模型對於數據中的雜訊太敏感，隱沒了大趨勢。要應付這個問題，我們可以使用懲罰樣條；於模型擬合時，樣條越曲折其懲罰越大，這誘使所得的模型，其樣條不至於過度擬合。使用懲罰樣條所得的模型一般擁有較平滑的表面，對於間中一兩個偏離值不會太敏感。

1 準確來說，這應該是經轉換後的機率差異相同。

4. 數據格式及候選模型

是次使用的數據集包含自 1961 至 2014 年2的每 6 小時熱帶氣旋位置，強度 (根據 10 分鐘平均風速標準) 及相應熱帶氣旋信號資料，共有 6,389 項觀測，並已經排除一長方區外、發出信號機率極低的觀測點。下表節錄 2013 年超級颱風天兔所對應的資料列：

名稱 年 月 日 時間
(UTC)		 現時 6 小時前
緯度 經度 強度 (節) 信號 緯度 經度 強度 (節) 信號
USAGI 2013 9 21 0 20.7 121.7 105 0 20.4 122.5 110 0
USAGI 2013 9 21 6 20.8 120.7 105 1 20.7 121.7 105 0
USAGI 2013 9 21 12 21.0 119.7 95 1 20.8 120.7 105 1
USAGI 2013 9 21 18 21.4 118.8 90 3 21.0 119.7 95 1
USAGI 2013 9 22 0 21.7 118.0 90 3 21.4 118.8 90 3
USAGI 2013 9 22 6 22.4 116.8 90 3 21.7 118.0 90 3
USAGI 2013 9 22 12 22.8 115.4 85 8 22.4 116.8 90 3
USAGI 2013 9 22 18 23.1 113.9 70 8 22.8 115.4 85 8
USAGI 2013 9 23 0 23.7 112.6 40 8 23.1 113.9 70 8
USAGI 2013 9 23 6 24.3 111.2 25 0 23.7 112.6 40 8

用作擬合九及十號信號資料的數據集格式相似，但時間間隔為 2 小時，該數據集共有 1,234 項觀測。

我們考慮的候選模型如下：

模型 關係方程式 參數數量
M0 \begin{eqnarray*} (時間 t 之信號狀態) & \sim & (時間 t 之經、緯、強度樣條) + (時間 t-6 之信號狀態) \end{eqnarray*} 76
M1 \begin{eqnarray*} (時間 t 之信號狀態) & \sim & (時間 t 之經、緯、強度樣條) + (時間 t 之徑向速率樣條) + (時間 t-6 之信號狀態) \end{eqnarray*} 85
M2 \begin{eqnarray*} (時間 t 之信號狀態) & \sim & (時間 t 之經、緯、強度樣條) + (時間 t-6 之經、緯、強度樣條) \\ & & + (時間 t-6 之信號狀態) \end{eqnarray*} 150
M3 \begin{eqnarray*} (時間 t 之信號狀態) & \sim & (時間 t-6 每個信號狀態之時間 t 經、緯、強度樣條) \\ & & + (時間 t-6 之經、緯、強度樣條) \end{eqnarray*} 224
M4 \begin{eqnarray*} (時間 t 之信號狀態) & \sim & (時間 t-6 每個信號狀態之時間 t 經、緯、強度樣條) \\ & & + (時間 t-6 每個信號狀態之時間 t-6 經、緯、強度樣條) \end{eqnarray*} 298
M5 \begin{eqnarray*} (時間 t 之信號狀態) & \sim & (時間 t-6 每個信號狀態之時間 t 經、緯度樣條) \\ & & + (時間 t-6 之經、緯度樣條) + (時間 t 之強度樣條) \end{eqnarray*} 82

模型 M0 可視作基準模型，用以對比其他模式的改善程度。留意模型 M0 至 M4 使用時間 $t$ (或 $t-6$) 之經、緯、強度三維樣條，而模型 M5 則只使用經、緯度二維樣條。模型 M1 利用徑向速率來代替氣旋先前位置。模型 M3 至 M5 有每個信號狀態之樣條，即時間 $t-6$ 信號並未生效時有一組樣條，而時間 $t-6$ 信號生效時有另一組，這令參數數量加倍。由於此組合已包含信號狀態之 (線性) 作用，該三個模型並沒有時間 $t-6$ 信號狀態之獨立項。模型 M0 至 M4 的參數數量逐步增加，而模型 M5 的參數數目比 M0 多，但比其餘少。

利用統計軟件 R，我們為每個信號擬合上述模型。每組數據中的反應變量為代表該信號 (或更高) 是否發出之二元變量。由於我們使用懲罰樣條，參數數值大小受到限制，模型之實質彈性比所列出的參數數量為低。我們使用交叉驗證來決定樣條懲罰程度，此步驟於下一節說明。

九及十號信號的模型擬合過程大致相同，但上述所有時間 $t-6$ 的變量皆轉為時間 $t-2$。

2 為提高準確度，數據集已更新至 2020 年的資料，並使用下節所選取的模型進行擬合。

5. 交叉驗證及模型驗證

使用懲罰樣條時，如果對擬合結果曲折度的懲罰太大會令擬合表面過於平滑，但如懲罰太小則有機會令擬合表面捕捉過多雜訊。我們可利用交叉驗證來選擇最佳懲罰程度，其中$k$-折交叉驗證須將數據集隨機分割為$k$等份，首次模型擬合使用除第一份外的所有數據，而擬合得出的模型則用來預測第一份數據的反應變數 (即估計機率)。重複此步驟$k$次，每次剔除第$i$份數據 ($i=2,3,\ldots,k$)，從而得出$k$組預測。利用一決定法則，即如預測機率高於某數便將反應變數視為 TRUE 值，將預測機率轉化為 TRUE/FALSE 的二元變數，並比對數據集中的反應變數，計算出誤辨率。經交叉驗證得出的誤辨率較能反映擬合模型的預測誤差；如將模型擬合至整組數據再作預測，數據將被重複使用 (模型擬合及預測過程皆使用)，得出的誤辨率無法反映模型套用至將來的數據時的實際誤差。

是次模型擬合期間，程式自動作交叉驗證運算，懲罰程度為得出最小誤辨率的數值。但是，由於數據集的反應變數中的 TRUE 值 (即信號正在生效) 比例偏低，我們可以使用其他指標來比較各模型的表現，而這些指標較為著重模型能否辨認 TRUE 值。設觀測及預測信號狀況的矩陣如下：

數目 觀測值為 FALSE 觀測值為 TRUE 總數
預測值為 FALSE $X_{FF}$ $X_{FT}$ $X_{F+}$
預測值為 TRUE $X_{TF}$ $X_{TT}$ $X_{T+}$
總數 $X_{+F}$ $X_{+T}$ $X_{++}$

這些指標分別為：

表現較佳模型的 CSI 或 HIT 值較高，而 FAR 值則較低。針對每個信號及每個擬合模型，我們使用「訓練/驗證集」的方法來計算此三指標的數值：從數據集中隨機選出三份之二為訓練集，其餘三份之一為驗證集。候選模型被擬合至訓練集，得出模型用來預測驗證集中的資料。我們重複此步驟三次，並計算每個信號/模型組合的平均指標值。是次使用的決定法則是將預測機率大於 0.5 的資料點歸類為 TRUE。驗證得出結果如下；每欄最佳的數字以粗體顯示：

模型 信號 1+ 信號 3+ 信號 8+ 信號 9+ 信號 10
CSI HIT FAR CSI HIT FAR CSI HIT FAR CSI HIT FAR CSI HIT FAR
M0 0.79 0.84 0.06 0.70 0.75 0.10 0.63 0.68 0.10 0.63 0.72 0.17 0.69 0.81 0.17
M1 0.80 0.84 0.06 0.72 0.77 0.08 0.61 0.66 0.10 0.63 0.71 0.16 0.69 0.80 0.17
M2 0.80 0.85 0.07 0.72 0.78 0.09 0.62 0.67 0.11 0.67 0.75 0.13 0.68 0.78 0.15
M3 0.78 0.82 0.06 0.70 0.75 0.09 0.53 0.57 0.12 N/A N/A
M4 0.78 0.82 0.05 0.70 0.75 0.07 0.50 0.53 0.11 N/A 0.70 0.85 0.20
M5 0.79 0.83 0.06 0.70 0.76 0.10 0.60 0.62 0.06 0.70 0.75 0.09 0.78 0.86 0.10

我們從一、三、八號信號中選取一個模型 (M2)；模型 M2 的 CSI 及 HIT 普遍較高，而 FAR 亦與其他模型相若。至於九及十號信號，模型 M5 的表現明顯較佳。留意由於模型 M3 及 M4 所含參數較多，它們的表現較不穩定，間中更未能成功擬合。

6. 擬合所得模型的意義

下文以一、三、八號信號的分析中所選取的 M2 模型來解說其用途。這個模型中，我們利用以下解釋變數來預測某時間 ($t$) 時一號或更高信號發出的機率：

對於時間 $t$ 時一網格中的每個經緯度，我們預測發出信號的機率，集合成擬合得出的機率表面。在以下第一個例子中，假設 6 小時前熱帶氣旋的位置為 18°N, 119°E (即圖一中各圖的粉紅圈位置)，當時強度為 55 節 (每小時 102 公里，強烈熱帶風暴級)。圖一中間假設 6 小時前信號並未發出、氣旋的強度於時間 $t-6$ 至 $t$ 間維持不變時的機率表面。由於氣旋集結於香港東南方，圖中西北方的位置表示氣旋於過去 6 小時移近本港，發出信號的機率理所當然較高。但是詮釋此圖需要小心；圖一中的箭頭顯示數據集中所有由 18°N, 119°E 周邊 0.5 經緯度內開始的 6 小時路徑，清楚可見熱帶氣旋於 6 小時內移動距離不多於 300 公里。因此圖中的機率表面只於此範圍內有代表性，範圍外的機率有外推成分，誤差或會較大，而預測機率的準繩度與該地歷史資料多寡有密切關係。

圖一 - 上文及下文所述之機率表面，圖中內弧與外弧分別標示香港 400 及 800 公里半徑範圍 (香港的位置為 22.3°N, 114.2°E)，而右下角可見菲律賓北部的地形輪廓。按圖可以原大小顯示。

圖一左圖與中間類似，但我們假設 6 小時前信號已經發出。可以見到此圖中紅色 (即高機率) 區域明顯較大。由於假設信號已經發出，加上氣旋正在靠近，信號繼續生效的機率自然較高。注意圖中 800 公里圈南端處預測機率仍較高，但由於該區缺乏觀測資料 (即很少熱帶氣旋以高速向西南移動)，此機率沒有代表性。最後，圖一右圖與中間唯一的分別為右圖假設氣旋於 6 小時內由 55 節增強至 70 節 (每小時 130 公里，颱風級)。此情況下機率較中間為高，但沒有像左圖、信號已經發出時般高。

於第二個例子中，我們假設一颱風 6 小時前位置為 20°N, 114°E，強度 70 節 (每小時 130 公里)，並預測發出八號或更高信號的機率。同樣，圖二中間假設 6 小時前並未發出八號或更高信號，且氣旋強度維持不變。由此圖可見，除非氣旋北移靠近本港使其威脅增加，否則於此 6 小時間改發八號或更高信號的機會不大。圖二左圖假設信號 6 小時前已經發出，明顯可見此 6 小時間八號或更高信號很大機會維持；右圖則假設 6 小時前信號並未生效，且氣旋由 70 節減弱至 50 節 (每小時 93 公里，強烈熱帶風暴級)，可見信號生效機率比中間低。圖二中的箭頭顯示數據集中所有由 20°N, 114°E 周邊 0.5 經緯度內開始的 6 小時路徑，對於預測機率有效範圍有參考作用。

圖二 - 上文所述之機率表面，圖中內弧與外弧分別標示香港 200 及 400 公里半徑範圍 (香港的位置為 22.3°N, 114.2°E)。按圖可以原大小顯示。

九及十號信號方面，擬合得出的模型 M5 詮釋方法類似，模型所需解釋變數相同，但時間間隔為 2 小時。

7. 從單週期的機率計算出多週期的機率

上述討論集中於單週期 (一個時間間隔) 的機率，即基於 6 或 2 小時前的資訊預測現在的信號機率。但如有一整條路徑預測，我們可以計算每一週期的機率，再整合為多週期的機率。要了解此計算方法，我們以本站於 2016 年 10 月 19 日為超級颱風海馬發出的路徑預測作一例子：

為保持一致，所有的強度預測將被換算為 10 分鐘平均，而每 6 小時的路徑預測以線性插值法計算。基於擬合所得的模型 M2，每 6 小時之三號或更高信號生效預測機率如下 (上列為假設 6 小時前並未發出有關信號時的機率，下列為假設 6 小時已發出有關信號時的機率)：

時間 $t-6$ 及 $t$ 間之單週期發出信號機率
時間 $t$ (小時) 6 12 18 24 30 36 42 48 54 60 66 72
三號或
更高		 時間 $t-6$ 並未發出 0.000 0.004 0.014 0.062 0.206 0.648 0.889 0.981 0.523 0.152 0.013 0.000
時間 $t-6$ 已發出 0.000 0.334 0.658 0.898 0.972 0.996 0.999 1.000 0.993 0.960 0.634 0.052

設 $Y_t$ 為時間 $t$ 時三號或更高信號是否生效的二元變數 (並未生效為 0，正在生效為 1)。利用全概率公式，我們可用以下方法計算多週期生效機率：

\begin{eqnarray*} \mathbb{P}(Y_6 = 1) & = & \mathbb{P}(Y_6 = 1 | Y_0 = y_0); \\ \mathbb{P}(Y_6 = 0) & = & 1-\mathbb{P}(Y_6 = 1); \\ \mathbb{P}(Y_t = 1) & = & \mathbb{P}(Y_t = 1 | Y_{t-6} = 0)\mathbb{P}(Y_{t-6} = 0) + \mathbb{P}(Y_t = 1 | Y_{t-6} = 1)\mathbb{P}(Y_{t-6} = 1); \\ \mathbb{P}(Y_t = 0) & = & 1-\mathbb{P}(Y_t = 1), \end{eqnarray*}

式中 $t$ 的值可以是 $6,12, \ldots ,72$，$y_0$ 為時間 0 的信號狀態 (已觀測)，而直線 ($|$) 意思為「給定」(given that)。因此，這些多週期機率與路徑預測發佈時信號是否生效有直接關係。

另一方面，給定某信號尚未發出 (取消)，我們或想得知信號首次於 $t-6$ 至 $t$ 時之間發出的機率。此機率可由下式得出：

\begin{eqnarray} \mathbb{P}(Y_{t} = 1, Y_{t-6} = 0, \ldots , Y_6 = 0 | Y_0 = 0) & = & \mathbb{P}(Y_{t} = 1 | Y_{t-6} = 0) \prod_{i=6}^{t-6} \mathbb{P}(Y_{i} = 0 | Y_{i-6} = 0); \label{eq:firstissue}\\ \mathbb{P}(Y_{t} = 0, Y_{t-6} = 1, \ldots , Y_6 = 1 | Y_0 = 1) & = & \mathbb{P}(Y_{t} = 0 | Y_{t-6} = 1) \prod_{i=6}^{t-6} \mathbb{P}(Y_{i} = 1 | Y_{i-6} = 1), \label{eq:firstcancel}\\ \end{eqnarray}

式中 $t$ 的值可以是 $6,12, \ldots ,72$。由於 $Y_0$ 已觀測，算式 \eqref{eq:firstissue} 及 \eqref{eq:firstcancel} 中只有一式有用。同時留意此機率並不等於信號於 $t-6$ 至 $t$ 時發出的機率，即 $\mathbb{P}(Y_{t} = 1 | Y_{t-6} = 0) \mathbb{P}(Y_{t-6} = 0)$；這是由於信號或於 $t-6$ 時之前曾經生效，但除非氣旋移動路徑迂迴，否則此兩機率差異不大。

利用這些算式，我們可以將上表的單週期機率轉換為下表的多週期機率。由於信號於時間 0 並未生效，$y_0=0$。信號生效及首次發出的機率如下：

多週期機率
時間 $t$ (小時) 6 12 18 24 30 36 42 48 54 60 66 72
三號或
更高		 時間 $t$ 正在生效 0.000 0.004 0.017 0.076 0.264 0.740 0.970 0.999 0.993 0.954 0.605 0.032
時間 $t-6$ 及 $t$ 之間首次發出 0.000 0.004 0.014 0.061 0.189 0.474 0.229 0.028 0.000 0.000 0.000 0.000

8. 對路徑的不確定性作出機率調整

上一節的機率之基於預測路徑為正確的假設得出。由於熱帶氣旋不是依循固定路徑移動，我們必須考慮路徑不確定性對信號發出機率的影響。利用本站歷史路徑誤差的經驗分布，我們模擬一些跟預測路徑有相近形狀但方向稍為更改的假設路徑。以上述的超級颱風海馬為例，模擬路徑圖如下：

除預測路徑外，上圖共有 80 條模擬路徑。於每一時間點的模擬氣旋位置形成 5 個同心圓，分別代表路徑距離誤差分布的第 10、30、50、70 及 90 百分位數。我們以第 7 節所述的方法為每條模擬路徑計算出其對應的多週期機率，然後取平均以反映路徑不確定性所帶來的加權機率。由於模擬點與預測點的距離為間距相等的百分位數，我們只需取簡單平均以得出此加權機率：設 $W$ 為一隨機變量 (此例子中，$W$ 為與預測點之距離)，密度函數及分布函數分別為 $f_W(w)$ 及 $F_W(w)$。設 $g$ 為一函數 (此例子中，$g$ 將距離轉換為根據擬合模型得出的預測機率)，那 $g(W)$ 的期望值為

\begin{equation} \mathbb{E}[g(W)]=\int_0^\infty g(w) f_W(w) \, \textrm{d}w = \int_0^\infty g(w) \, \textrm{d}F_W(w) = \int_0^1 g(F^{-1}(u)) \,\textrm{d}u, \label{eq:expectquant} \end{equation}

上式中，$F$ 的反函數為分位函數 (quantile function)，我們以樣本分位函數 (sample quantile) 取代之。由於地球上有多於一點與預測點距離相等，我們於上述路徑模擬中考慮所有方向，對每個分位 $u$ 的所有路徑所得出的機率取平均數，而最後結果就是將所有路徑得出的多週期機率取簡單平均。由於只取 5 個分位點，上述計算只能得到 \eqref{eq:expectquant} 的近似值，但我們經試驗確認此近似值表現已經不錯，即使取更多分位點及模擬路徑，所得出的機率變化不大。

針對上面的例子，其平均機率詳列於下表，這些數字亦曾於本站發佈中顯示：

對路徑不確定性作出調整後的平均多週期機率
時間 $t$ (小時) 6 12 18 24 30 36 42 48 54 60 66 72
三號或
更高		 時間 $t$ 正在生效 0.000 0.004 0.020 0.110 0.318 0.632 0.853 0.915 0.869 0.741 0.448 0.135
時間 $t-6$ 及 $t$ 之間首次發出 0.000 0.004 0.017 0.092 0.210 0.314 0.223 0.072 0.005 0.001 0.000 0.000

有關詮釋本產品的方法，請參閱此頁。

最近訪問日期: Sun Sep 08 2024 13:03:22 HKT
最近修訂日期: Wed Jun 29 2022