報(bào)對比到排行榜的Python實(shí)現(xiàn))
天氣模型排名指的是把不同數(shù)值天氣預(yù)報(bào)模型放到同一批歷史天氣場景里用實(shí)際觀測結(jié)果檢驗(yàn)它們的預(yù)報(bào)誤差最后按誤差大小排出先后。這個(gè)做法在氣象服務(wù)、能源調(diào)度、農(nóng)業(yè)決策、物流規(guī)劃和量化交易場景里非常有用業(yè)務(wù)方不需要關(guān)心模型內(nèi)部有多少物理過程只需要知道過去一段時(shí)間里哪個(gè)模型的預(yù)報(bào)更接近真實(shí)天氣從而決定后續(xù)引用哪條預(yù)報(bào)源。標(biāo)題里的 Show HN: Ranking weather models by how their forecasts turned out 對應(yīng)的正是這樣一個(gè)系統(tǒng)把“事后驗(yàn)證”做成排行榜讓模型質(zhì)量可比較、可追蹤、可復(fù)盤。下面會(huì)圍繞如何搭建一套最小可用的天氣模型排名系統(tǒng)展開。主線是從數(shù)據(jù)準(zhǔn)備、對齊邏輯、評分指標(biāo)到 Python 實(shí)現(xiàn)、異常排查和工程化擴(kuò)展。所有代碼和數(shù)據(jù)格式都是演示結(jié)構(gòu)用于說明思路真實(shí)項(xiàng)目接入任何模型輸出、觀測站點(diǎn)或數(shù)據(jù)文件時(shí)需要按實(shí)際字段和授權(quán)調(diào)整。1. 天氣預(yù)報(bào)模型排名要解決什么問題1.1 什么是天氣模型排名為什么不能只信模型說明數(shù)值天氣預(yù)報(bào)模型會(huì)用物理方程描述大氣運(yùn)動(dòng)再通過超級計(jì)算機(jī)求解得到未來若干小時(shí)甚至十幾天的氣象要素預(yù)報(bào)。不同模型在物理方案、分辨率、資料同化方式上差別很大所以同一天、同一個(gè)地點(diǎn)、同一個(gè)變量的預(yù)報(bào)結(jié)果經(jīng)常不一致。使用者最關(guān)心的問題是哪個(gè)模型更可信天氣模型排名就是給這個(gè)問題提供一個(gè)可量化的答案。它不評價(jià)模型內(nèi)部的物理設(shè)計(jì)只評價(jià)“歷史輸出和真實(shí)觀測的差距”。通俗地說一次預(yù)報(bào)就是一次考試觀測就是標(biāo)準(zhǔn)答案。模型在大量歷史樣本上考出來的平均得分就是它的排名依據(jù)。技術(shù)定義可以寫成給定一組模型、一組起報(bào)時(shí)間、一組站點(diǎn)或格點(diǎn)、一組氣象變量和一組驗(yàn)證時(shí)間逐個(gè)比較預(yù)報(bào)值與觀測值計(jì)算統(tǒng)計(jì)算法指標(biāo)再按指標(biāo)排序。不能只信模型說明原因主要有三點(diǎn)。第一模型宣傳材料通常強(qiáng)調(diào)分辨率、同化系統(tǒng)等能力卻不一定提供與業(yè)務(wù)場景匹配的獨(dú)立驗(yàn)證。第二模型在不同地區(qū)、季節(jié)和天氣類型下的表現(xiàn)可能差異很大一個(gè)總體分?jǐn)?shù)不能覆蓋所有場景。第三模型版本會(huì)更新排名會(huì)變化只有持續(xù)驗(yàn)證才能跟蹤真實(shí)變化。一個(gè)榜單系統(tǒng)的價(jià)值不是給出一次結(jié)論而是讓結(jié)論可以被持續(xù)更新和質(zhì)疑。1.2 排名系統(tǒng)的核心鏈路預(yù)報(bào)、觀測、對齊、評分一套最小可用的天氣模型排名系統(tǒng)核心鏈路只有五步。第一步收集多個(gè)模型在過去一段時(shí)間的預(yù)報(bào)。注意“預(yù)報(bào)”不是一句話而是結(jié)構(gòu)化數(shù)據(jù)至少要包含起報(bào)時(shí)間、有效時(shí)間、站點(diǎn)或經(jīng)緯度、氣象變量和預(yù)報(bào)值。第二步收集同時(shí)段的觀測數(shù)據(jù)。觀測可以是氣象站、網(wǎng)格分析場或衛(wèi)星反演產(chǎn)品但要和預(yù)報(bào)使用同一套空間和時(shí)間口徑。第三步對齊。這一步最容易被忽略卻最影響結(jié)果。兩個(gè)模型必須比較同一個(gè)有效時(shí)間、同一個(gè)地點(diǎn)、同一個(gè)變量的預(yù)報(bào)才公平。第四步計(jì)算誤差指標(biāo)。比如絕對誤差、均方根誤差、偏差、降水命中率等。第五步按模型聚合指標(biāo)并生成排行榜再展示給業(yè)務(wù)方或下游系統(tǒng)。后面所有代碼和配置都圍繞這條鏈路展開。很多人一開始就把精力放在“把模型跑起來”或“把畫圖做漂亮”上其實(shí)排名系統(tǒng)最核心的工作量在前三步數(shù)據(jù)怎么來、字段怎么對應(yīng)、口徑怎么統(tǒng)一。數(shù)據(jù)和字段錯(cuò)了后面的指標(biāo)再專業(yè)也沒有意義。1.3 適用場景與讀者定位這套方案適合以下幾類讀者。第一類是氣象數(shù)據(jù)產(chǎn)品開發(fā)人員需要為業(yè)務(wù)方提供模型質(zhì)量看板。第二類是數(shù)據(jù)工程師正在對接氣象預(yù)報(bào)和觀測數(shù)據(jù)需要設(shè)計(jì)統(tǒng)一驗(yàn)證流程。第三類是研究人員想快速驗(yàn)證新模型或新參數(shù)方案是否優(yōu)于既有模型。第四類是后端開發(fā)人員需要把歷史驗(yàn)證結(jié)果做成排行榜頁面或 API。本文的例子使用 Python 和 pandas不依賴專業(yè)氣象軟件。這樣做的目的是先把排名邏輯講清楚再讓讀者根據(jù)自己的數(shù)據(jù)格式替換加載層。如果你已經(jīng)在使用 xarray、cfgrib 或 GRIB2 文件只需要替換數(shù)據(jù)加載部分評分和排名邏輯可以原樣復(fù)用。2. 數(shù)據(jù)準(zhǔn)備歷史預(yù)報(bào)和觀測數(shù)據(jù)如何配對2.1 先統(tǒng)一時(shí)間口徑預(yù)報(bào)發(fā)布時(shí)間與有效時(shí)間氣象數(shù)據(jù)里至少有三個(gè)時(shí)間概念起報(bào)時(shí)間issue_time、預(yù)報(bào)時(shí)效lead_hour和有效時(shí)間valid_time。起報(bào)時(shí)間是模型開始計(jì)算的時(shí)刻通常是 00 時(shí)或 12 時(shí)預(yù)報(bào)時(shí)效是預(yù)測未來多少個(gè)小時(shí)后的大氣狀態(tài)有效時(shí)間是預(yù)報(bào)值對應(yīng)的真實(shí)時(shí)刻它等于起報(bào)時(shí)間加上預(yù)報(bào)時(shí)效。排名系統(tǒng)比較的是“對同一個(gè)真實(shí)時(shí)刻的預(yù)報(bào)誰更準(zhǔn)”所以最終必須以 valid_time 作為對齊鍵。如果兩個(gè)模型起報(bào)時(shí)間相同但一個(gè)預(yù)報(bào) 24 小時(shí)另一個(gè)預(yù)報(bào) 48 小時(shí)它們對應(yīng)的 valid_time 完全不同不能放在一起比較。同樣如果兩個(gè)模型 valid_time 相同但起報(bào)時(shí)間不同它們利用的初始資料可能不同可以比較但要注意時(shí)效差異。下面是一個(gè)典型字段示例model,issue_time,valid_time,lead_hour GFS,2025-01-01T00:00:00Z,2025-01-02T00:00:00Z,24 ECMWF,2025-01-01T00:00:00Z,2025-01-02T00:00:00Z,24在讀取數(shù)據(jù)時(shí)所有時(shí)間字段都應(yīng)該轉(zhuǎn)成帶時(shí)區(qū)的 UTC 時(shí)間。不要使用本地時(shí)間字符串做 join因?yàn)椴煌瑪?shù)據(jù)源可能使用不同時(shí)區(qū)容易出現(xiàn)“看起來相同實(shí)際差 8 小時(shí)”的問題。2.2 先統(tǒng)一空間口徑格點(diǎn)與站點(diǎn)如何對齊天氣預(yù)報(bào)模型通常輸出網(wǎng)格數(shù)據(jù)比如 0.25 度、0.5 度分辨率的經(jīng)緯度格點(diǎn)。氣象觀測站是離散的點(diǎn)可能落在兩個(gè)格點(diǎn)之間。為了比較需要把模型格點(diǎn)插值到觀測站點(diǎn)位置或者把觀測值匹配到網(wǎng)格上。最簡單可靠的做法是使用現(xiàn)成插值工具將模型值插值到站點(diǎn)然后把站點(diǎn)的經(jīng)緯度轉(zhuǎn)換成唯一的 station_id。在演示項(xiàng)目中為了聚焦排名邏輯可以先把數(shù)據(jù)整理成站點(diǎn)的 long 表。每條記錄都包含 model、valid_time、station_id、variable、value。對齊時(shí)forecast 和 observation 都通過 valid_time、station_id、variable 三個(gè)字段連接。這樣能避免插值實(shí)現(xiàn)干擾主流程。如果數(shù)據(jù)源來自網(wǎng)格加載時(shí)要注意經(jīng)緯度表示方式和單位。部分?jǐn)?shù)據(jù)用整數(shù)經(jīng)緯度部分用浮點(diǎn)經(jīng)緯度有的用 lon 表示 -180 到 180有的用 0 到 360。做空間對齊前需要先統(tǒng)一經(jīng)緯度范圍。2.3 用最小 CSV 結(jié)構(gòu)承載評分?jǐn)?shù)據(jù)為了能直接運(yùn)行后面的 Python 代碼這里定義兩個(gè)最小 CSV 文件。第一個(gè)是歷史預(yù)報(bào)文件 forecasts.csvmodel,issue_time,valid_time,station_id,variable,value GFS,2025-01-01T00:00:00Z,2025-01-02T00:00:00Z,S001,temperature_2m,4.2 ECMWF,2025-01-01T00:00:00Z,2025-01-02T00:00:00Z,S001,temperature_2m,5.0 ICON,2025-01-01T00:00:00Z,2025-01-02T00:00:00Z,S001,temperature_2m,3.9第二個(gè)是觀測文件 observations.csvvalid_time,station_id,variable,value 2025-01-02T00:00:00Z,S001,temperature_2m,3.5這兩個(gè)文件都很小但已經(jīng)包含了排名系統(tǒng)需要的所有核心字段。變量名建議統(tǒng)一使用機(jī)器可讀的小寫字符例如 temperature_2m、wind_speed_10m、precipitation_24h。如果值有單位最好在表里增加 unit 字段或者先在預(yù)處理階段統(tǒng)一單位避免不同模型輸出攝氏度和華氏度造成誤差計(jì)算錯(cuò)誤。注意觀測文件的粒度是“真實(shí)時(shí)刻的客觀值”而預(yù)報(bào)文件還需要保留 model 和 issue_time。合并時(shí)只保留 valid_time 相同的記錄否則樣本對不齊排名會(huì)失真。2.4 數(shù)據(jù)質(zhì)量檢查清單在開始寫評分代碼前建議先對數(shù)據(jù)做一輪質(zhì)量檢查。下面這個(gè)清單可以直接復(fù)用檢查項(xiàng)檢查方式通過標(biāo)準(zhǔn)時(shí)間字段格式打印每個(gè)來源的 unique 時(shí)間樣例全部使用 ISO8601 且?guī)?Z 時(shí)區(qū)valid_time 覆蓋范圍檢查 min/max預(yù)報(bào)和觀測覆蓋同一時(shí)間段station_id 是否一致計(jì)算集合差預(yù)報(bào)和觀測站點(diǎn)集合一致variable 是否一致計(jì)算集合差變量名完全匹配是否存在重復(fù)記錄按 key 分組統(tǒng)計(jì)每個(gè) key 只有一條記錄值是否缺測查看 value 為 NaN 或特殊值缺測數(shù)量低于閾值且被標(biāo)記這個(gè)檢查可以用 pandas 在預(yù)處理階段自動(dòng)完成。不要省略這些步驟因?yàn)闀r(shí)間或空間對齊一旦出錯(cuò)排名結(jié)果會(huì)出現(xiàn)系統(tǒng)性偏差而且肉眼很難發(fā)現(xiàn)。3. 用 Python 實(shí)現(xiàn)最小可運(yùn)行的模型評分與排名流程3.1 項(xiàng)目結(jié)構(gòu)和依賴在本地目錄創(chuàng)建一個(gè)最小項(xiàng)目weather-leaderboard/ ├── data/ │ ├── forecasts.csv │ └── observations.csv ├── leaderboard.py ├── requirements.txt └── run.sh依賴只需要 pandas 和 numpy畫圖可以后面再加。安裝命令pip install pandas numpy matplotlibrequirements.txt 可以寫成pandas2.0 numpy1.24 matplotlib3.7這一步的目的是讓示例可以快速運(yùn)行。真實(shí)項(xiàng)目中如果使用 GRIB2 或 NetCDF需要額外引入 xarray、cfgrib 或 h5netcdf但核心評分邏輯不受影響。3.2 加載和對齊數(shù)據(jù)在 leaderboard.py 中先讀取兩個(gè) CSV 文件import pandas as pd def load_data(forecast_path, observation_path): forecasts pd.read_csv(forecast_path, parse_dates[issue_time, valid_time]) observations pd.read_csv(observation_path, parse_dates[valid_time]) return forecasts, observations然后寫對齊函數(shù)def align_forecast_observation(forecasts, observations): merged pd.merge( forecasts, observations, on[valid_time, station_id, variable], suffixes(_fc, _obs), howinner ) merged[error] merged[value_fc] - merged[value_obs] return merged這里使用內(nèi)連接強(qiáng)制保留預(yù)報(bào)和觀測都存在的樣本。如果某個(gè)模型出現(xiàn)了數(shù)據(jù)缺口它的樣本量會(huì)比其他模型小這時(shí)排名結(jié)果需要謹(jǐn)慎解讀。使用 suffixes 區(qū)分預(yù)報(bào)值和觀測值error 表示預(yù)報(bào)減觀測的差值正偏差表示預(yù)報(bào)偏高。對齊之后需要檢查一下樣本數(shù)量def sample_summary(df): return df.groupby(model).size().reset_index(namesample_size)3.3 計(jì)算誤差指標(biāo)接下來定義三個(gè)連續(xù)變量誤差指標(biāo)平均絕對誤差 MAE、均方根誤差 RMSE、平均偏差 Bias。import numpy as np def mae(obs, fct): obs np.asarray(obs, dtypefloat) fct np.asarray(fct, dtypefloat) return float(np.mean(np.abs(fct - obs))) def rmse(obs, fct): obs np.asarray(obs, dtypefloat) fct np.asarray(fct, dtypefloat) return float(np.sqrt(np.mean((fct - obs) ** 2))) def bias(obs, fct): obs np.asarray(obs, dtypefloat) fct np.asarray(fct, dtypefloat) return float(np.mean(fct - obs))這三個(gè)函數(shù)都比較簡單。MAE 直觀反映平均誤差大小RMSE 對大誤差更敏感因?yàn)槠椒巾?xiàng)放大了離群值的影響B(tài)ias 表示系統(tǒng)性誤差正值偏暖或偏強(qiáng)負(fù)值偏冷或偏弱。排名時(shí)通常以 MAE 或 RMSE 為主指標(biāo)Bias 作為輔助參考。3.4 按模型生成排行榜核心邏輯是按 model 分組對每組樣本計(jì)算指標(biāo)再排序def rank_models(merged, metricmae): records [] for model, group in merged.groupby(model): obs group[value_obs] fct group[value_fc] records.append({ model: model, sample_size: len(group), mae: mae(obs, fct), rmse: rmse(obs, fct), bias: bias(obs, fct), }) result pd.DataFrame(records) # 這里按誤差指標(biāo)從小到大排名bias 不作為主排名指標(biāo) result[rank] result[metric].rank(methodmin, ascendingTrue).astype(int) result result.sort_values([rank, model]).reset_index(dropTrue) return result這段代碼先遍歷每個(gè)模型再生成 DataFrame。metric 參數(shù)可以切換排名用的指標(biāo)默認(rèn) mae。rank 列使用 pandas 的 rank 方法相同誤差時(shí)并列名次。這里要注意樣本量暫時(shí)只作為展示不參與排序。如果兩個(gè)模型樣本量差異很大需要在業(yè)務(wù)層決定是否過濾。主流程如下if __name__ __main__: forecasts, observations load_data(data/forecasts.csv, data/observations.csv) merged align_forecast_observation(forecasts, observations) if merged.empty: raise SystemExit(沒有匹配到有效樣本請檢查時(shí)間、站點(diǎn)、變量字段) table rank_models(merged, metricmae) print(table.to_string(indexFalse))3.5 運(yùn)行結(jié)果示例如果 forecasts.csv 中包含多個(gè)站點(diǎn)和多天數(shù)據(jù)運(yùn)行后可能得到類似這樣的輸出model sample_size mae rmse bias rank ECMWF 2800 1.82 2.41 0.35 1 GFS 2800 2.05 2.78 0.10 2 ICON 2800 2.31 3.02 -0.42 3這個(gè)結(jié)果的業(yè)務(wù)