1. 為什么要做一個(gè)基于Python的大氣污染預(yù)測(cè)軟件我先交代一下背景。去年接了個(gè)空氣質(zhì)量監(jiān)測(cè)的課題需要根據(jù)歷史污染物濃度數(shù)據(jù)預(yù)測(cè)未來(lái)幾天的PM2.5、PM10、NO?、O?等指標(biāo)。最開(kāi)始用Excel硬算滑動(dòng)平均、線性回歸都試過(guò)結(jié)果一到污染過(guò)程轉(zhuǎn)折點(diǎn)就嚴(yán)重滯后。后來(lái)決定自己動(dòng)手用Python寫一套完整的時(shí)間序列分析預(yù)測(cè)軟件把數(shù)據(jù)清洗、特征構(gòu)造、模型訓(xùn)練、預(yù)測(cè)、可視化、報(bào)告輸出全部串起來(lái)。這套東西做完以后不僅課題交付了還整理出了一份配套文檔和源碼工程方便后面任何人接手維護(hù)。文章標(biāo)題里的“基于Python的時(shí)間序列分析的大氣污染預(yù)測(cè)軟件”說(shuō)白了就是干這件事的。它能做的事情很明確給定一個(gè)監(jiān)測(cè)站點(diǎn)的歷史濃度數(shù)據(jù)自動(dòng)完成缺失值處理、異常值識(shí)別、趨勢(shì)和周期性分解然后訓(xùn)練時(shí)間序列模型輸出未來(lái)24小時(shí)到72小時(shí)的濃度預(yù)測(cè)并生成圖表和評(píng)估指標(biāo)。適合給環(huán)境監(jiān)測(cè)站、高校實(shí)驗(yàn)室、環(huán)保方向的學(xué)生和科研人員做參考也可以作為課程設(shè)計(jì)或畢業(yè)設(shè)計(jì)的完整項(xiàng)目模板。我見(jiàn)過(guò)太多類似項(xiàng)目最終卡死在兩個(gè)地方一是數(shù)據(jù)質(zhì)量太差二是模型選型過(guò)于隨意。很多教程上來(lái)就教你跑ARIMA但壓根不告訴你數(shù)據(jù)需要滿足什么前提也不告訴你如果數(shù)據(jù)有多個(gè)季節(jié)周期該怎么辦。所以這篇博文我不想只貼一堆代碼而是把完整的設(shè)計(jì)思路和踩坑過(guò)程都講清楚讓你拿到源碼之后不僅能跑通還能根據(jù)實(shí)際數(shù)據(jù)調(diào)整出更好的效果。2. 時(shí)間序列分析在空氣質(zhì)量預(yù)測(cè)里的核心思路2.1 為什么用時(shí)間序列而不是普通回歸大氣污染數(shù)據(jù)天然是帶時(shí)間戳的序列數(shù)據(jù)每小時(shí)的監(jiān)測(cè)值之間存在強(qiáng)烈的自相關(guān)性。舉個(gè)很直觀的例子今天下午3點(diǎn)的PM2.5濃度和昨天下午3點(diǎn)、前幾小時(shí)的數(shù)據(jù)密切相關(guān)而不是和某個(gè)不相關(guān)的特征比如“今天星期幾”線性相關(guān)。普通回歸模型把每個(gè)時(shí)刻當(dāng)作獨(dú)立樣本處理會(huì)徹底丟掉這種時(shí)間依賴關(guān)系。時(shí)間序列分析的核心思路是把“時(shí)間順序”作為模型的一等公民。我們會(huì)先處理三件事趨勢(shì)trend、季節(jié)性seasonality、殘差residual。空氣質(zhì)量數(shù)據(jù)里有明顯的日周期交通早晚高峰導(dǎo)致污染物濃度升高、周周期周末工業(yè)活動(dòng)減少和年周期冬季采暖導(dǎo)致濃度上升。如果不用時(shí)間序列方法這些周期性規(guī)律很難被普通回歸模型捕捉。另一個(gè)原因是預(yù)測(cè)任務(wù)本身的性質(zhì)。我們要做的是“未來(lái)幾小時(shí)/幾天的濃度”不是“給定一堆特征預(yù)測(cè)一個(gè)值”。時(shí)間序列模型天生支持滾動(dòng)預(yù)測(cè)比如用過(guò)去72小時(shí)預(yù)測(cè)未來(lái)1小時(shí)再逐步滾動(dòng)預(yù)測(cè)未來(lái)24小時(shí)。這對(duì)環(huán)境監(jiān)測(cè)、預(yù)警場(chǎng)景極其重要。2.2 數(shù)據(jù)分解看懂序列里藏著的規(guī)律拿到數(shù)據(jù)后第一件事不是建模而是分解。我常用statsmodels的seasonal_decompose把序列拆成趨勢(shì)項(xiàng)、季節(jié)項(xiàng)和殘差項(xiàng)。import pandas as pd import statsmodels.api as sm df pd.read_csv(pm25_hourly.csv, parse_dates[time], index_coltime) decomp sm.tsa.seasonal_decompose(df[pm25], modeladditive, period24) decomp.plot()這里要特別注意period參數(shù)。如果不指定seasonal_decompose默認(rèn)認(rèn)為沒(méi)有季節(jié)周期很多新手會(huì)漏掉這一步?？諝赓|(zhì)量數(shù)據(jù)的日周期是24小時(shí)周周期是168小時(shí)如果做月級(jí)別的預(yù)測(cè)可能還需要考慮年周期8760小時(shí)。分解結(jié)果能直觀告訴你序列是否有明顯趨勢(shì)是否有周期性殘差是否平穩(wěn)這是后續(xù)選擇模型的基礎(chǔ)。2.3 平穩(wěn)性檢驗(yàn)為什么ARIMA要求數(shù)據(jù)平穩(wěn)經(jīng)典ARIMA模型的前提是序列平穩(wěn)也就是均值、方差不隨時(shí)間變化。空氣污染物濃度明顯不滿足——冬季高、夏季低早晚高峰高、午后低。所以直接對(duì)原始數(shù)據(jù)跑ARIMA結(jié)果肯定一塌糊涂。正確做法是先做ADF檢驗(yàn)Augmented Dickey-Fuller test判斷是否平穩(wěn)如果不平穩(wěn)就差分。我常用的代碼from statsmodels.tsa.stattools import adfuller result adfuller(df[pm25].dropna()) print(fADF統(tǒng)計(jì)量: {result[0]:.4f}) print(fp值: {result[1]:.4f}) # p值小于0.05通常認(rèn)為平穩(wěn)但這里有個(gè)陷阱空氣質(zhì)量數(shù)據(jù)單純差分一次往往還不夠因?yàn)檫€有季節(jié)性。如果只做一階差分季節(jié)性依舊存在ADF檢驗(yàn)可能仍然不平穩(wěn)。這時(shí)候需要考慮兩種路一是使用SARIMA對(duì)季節(jié)性部分做季節(jié)差分二是對(duì)原始序列做STL分解后對(duì)殘差項(xiàng)建模。我在項(xiàng)目里通常先做STL分解再對(duì)殘差做ARIMA效果比直接差分好很多。3. 數(shù)據(jù)準(zhǔn)備從爬取到清洗這一步?jīng)Q定了預(yù)測(cè)上限3.1 數(shù)據(jù)來(lái)源選擇與接口對(duì)接大氣污染數(shù)據(jù)來(lái)源一般是兩類一類是環(huán)境監(jiān)測(cè)站提供的CSV或Excel導(dǎo)出另一類是通過(guò)API實(shí)時(shí)抓取。常見(jiàn)的API比如和風(fēng)天氣、AirVisual以及國(guó)內(nèi)的青悅開(kāi)放數(shù)據(jù)平臺(tái)。當(dāng)時(shí)我做的軟件需要支持兩種方式讀本地文件以及定時(shí)從API拉取。以AirVisual API為例請(qǐng)求很簡(jiǎn)單import requests api_key 你的key url fhttps://api.airvisual.com/v2/city?cityBeijingstateBeijingcountryChinakey{api_key} resp requests.get(url).json()注意不同API返回字段差異很大。有的返回的是“當(dāng)前實(shí)時(shí)濃度”有的返回的是“過(guò)去24小時(shí)歷史濃度”。做時(shí)間序列預(yù)測(cè)必須要有足夠長(zhǎng)的歷史數(shù)據(jù)至少需要過(guò)去一個(gè)月以上的小時(shí)級(jí)數(shù)據(jù)否則模型根本學(xué)不到周期性。3.2 缺失值處理不能隨便刪也不能只填均值污染物監(jiān)測(cè)設(shè)備經(jīng)常斷線凌晨數(shù)據(jù)也容易異常。缺失值處理是數(shù)據(jù)處理環(huán)節(jié)最耗時(shí)的部分。簡(jiǎn)單粗暴地刪除缺失行會(huì)打斷時(shí)間的連續(xù)性后續(xù)做滯后特征、傅里葉變換時(shí)都會(huì)出問(wèn)題。簡(jiǎn)單填充均值會(huì)抹掉日周期特征。我實(shí)際使用的策略分三層場(chǎng)景處理方式理由單點(diǎn)缺失1小時(shí)線性插值相鄰時(shí)刻相關(guān)性最接近線性插值足夠連續(xù)缺失3~6小時(shí)時(shí)間加權(quán)插值或前后同周期均值考慮日周期用前后兩天的同一時(shí)刻插值連續(xù)缺失超過(guò)24小時(shí)刪除該段或分段訓(xùn)練長(zhǎng)缺失會(huì)引入大量估計(jì)誤差不如放棄具體代碼df[pm25] df[pm25].interpolate(methodlinear)后來(lái)我改進(jìn)了一下用scipy.interpolate的PchipInterpolator做保形狀插值避免線性插值導(dǎo)致的“尖角”突變。插值之后還要對(duì)所有填充點(diǎn)做標(biāo)記方便后面評(píng)估模型時(shí)區(qū)分真實(shí)值和填充值否則驗(yàn)證集的誤差會(huì)被虛假值拉低。3.3 異常值識(shí)別用滾動(dòng)標(biāo)準(zhǔn)差比閾值更靠譜空氣監(jiān)測(cè)數(shù)據(jù)經(jīng)常出現(xiàn)瞬時(shí)高值比如0點(diǎn)出現(xiàn)一個(gè)10000的值這明顯是設(shè)備故障。如果不去除模型會(huì)被這個(gè)點(diǎn)帶跑偏。我之前試過(guò)固定閾值法比如500就算異常但這個(gè)在污染嚴(yán)重的城市冬季根本不適用——正常濃度就能到400。后來(lái)改用滾動(dòng)窗口的中位數(shù)和MADMedian Absolute Deviation來(lái)檢測(cè)異常。MAD法對(duì)離群值更穩(wěn)健適合非正態(tài)分布的污染物數(shù)據(jù)。import numpy as np def detect_outliers(series, window24, n_sigma5): rolling_median series.rolling(windowwindow, centerTrue).median() mad (series - rolling_median).abs().rolling(windowwindow, centerTrue).median() z_score 0.6745 * (series - rolling_median) / mad return z_score.abs() n_sigma這個(gè)方法的物理意義是如果一個(gè)點(diǎn)的濃度水平顯著偏離周圍24小時(shí)的“正常波動(dòng)程度”就判定為設(shè)備異常用插值替換。注意n_sigma別設(shè)太小空氣質(zhì)量波動(dòng)本身就大設(shè)成3會(huì)誤殺很多真實(shí)的高濃度過(guò)程。3.4 特征工程讓模型看到時(shí)間節(jié)律雖然時(shí)間序列模型本身能捕捉時(shí)間依賴但適當(dāng)構(gòu)造外生特征能顯著提升預(yù)測(cè)精度。我加了幾組特征時(shí)間sin/cos編碼用正弦余弦編碼小時(shí)、星期、月份避免模型誤以為“23”和“0”兩個(gè)小時(shí)內(nèi)差別非常大。滯后值過(guò)去1、3、6、12、24小時(shí)的濃度。氣象因子溫度、濕度、風(fēng)速、氣壓如果API能取到。周邊站點(diǎn)數(shù)據(jù)相鄰站點(diǎn)的濃度對(duì)流場(chǎng)有指示作用。氣象因子的重要性經(jīng)常被忽視。有一次模型預(yù)測(cè)霧霾消散時(shí)間總是滯后2小時(shí)后來(lái)發(fā)現(xiàn)是沒(méi)加入風(fēng)速變化特征。風(fēng)速突然增大后污染物濃度會(huì)在下一小時(shí)迅速下降這個(gè)信號(hào)純靠歷史濃度序列反映不出來(lái)。4. 模型實(shí)現(xiàn)從ARIMA到LSTM怎么選怎么調(diào)4.1 基準(zhǔn)模型SARIMA是效率首選一開(kāi)始我直接上ARIMA結(jié)果被數(shù)據(jù)的光照季節(jié)性和周期性教做人。后來(lái)改用SARIMA即季節(jié)性ARIMA才算是真正入了門。SARIMA的全稱是Seasonal Autoregressive Integrated Moving Average它額外增加了季節(jié)項(xiàng)。模型參數(shù)寫成SARIMA(p,d,q)(P,D,Q,S)其中S是季節(jié)周期PM2.5小時(shí)數(shù)據(jù)的S24。確定p、d、q最常用的方法是我用過(guò)的pmdarima庫(kù)的auto_arima它會(huì)自動(dòng)搜索最佳參數(shù)組合。直接用示例from pmdarima import auto_arima model auto_arima( df[pm25].dropna(), seasonalTrue, m24, start_p0, max_p5, start_q0, max_q5, d1, D1, traceTrue, stepwiseTrue )這里要說(shuō)明兩點(diǎn)。第一d1, D1是我先做ADF檢驗(yàn)、再對(duì)季節(jié)差分后的序列做ADF檢驗(yàn)得到的結(jié)果不是隨便設(shè)的。第二m24意味著模型會(huì)把“24小時(shí)前”的誤差和滯后項(xiàng)納入計(jì)算訓(xùn)練時(shí)間會(huì)成倍增加。數(shù)據(jù)量超過(guò)1萬(wàn)條時(shí)auto_arima可能要跑十幾分鐘正?，F(xiàn)象。4.2 升級(jí)方案Prophet處理多季節(jié)周期SARIMA最多處理一個(gè)季節(jié)周期。但空氣質(zhì)量同時(shí)有日、周、年度周期。Facebook開(kāi)源的Prophet可以同時(shí)擬合多個(gè)周期而且對(duì)缺失值、異常值相當(dāng)魯棒。我第二版軟件把Prophet用作主力模型之一。Prophet的使用特別簡(jiǎn)單from prophet import Prophet model Prophet( yearly_seasonalityTrue, weekly_seasonalityTrue, daily_seasonalityTrue, changepoint_prior_scale0.05 ) model.add_regressor(wind_speed) model.fit(train) future model.make_future_dataframe(periods72, freqH) future[wind_speed] forecasted_weather forecast model.predict(future)雖然Prophet不算嚴(yán)格的“時(shí)間序列分析”工具它本質(zhì)上是可加回歸模型但拆分趨勢(shì)、周期、節(jié)假日的能力非常適合空氣質(zhì)量這種強(qiáng)季節(jié)性數(shù)據(jù)。唯一的問(wèn)題是它對(duì)突變點(diǎn)比如沙塵暴事件的擬合能力偏弱因?yàn)樽凕c(diǎn)默認(rèn)是稀疏先驗(yàn)。我調(diào)高了changepoint_prior_scale到0.1之后能捕捉到一些突發(fā)性污染過(guò)程但代價(jià)是過(guò)擬合風(fēng)險(xiǎn)也增加。4.3 進(jìn)階方向LSTM和它們的現(xiàn)實(shí)局限網(wǎng)上很多教程把LSTM吹得神乎其神實(shí)際用下來(lái)并沒(méi)有想象中那么美。LSTM的優(yōu)勢(shì)是能自動(dòng)學(xué)習(xí)非線性關(guān)系和長(zhǎng)時(shí)間依賴?yán)碚撋虾苓m合空氣質(zhì)量預(yù)測(cè)但需要大量數(shù)據(jù)。我實(shí)測(cè)發(fā)現(xiàn)當(dāng)訓(xùn)練數(shù)據(jù)少于3萬(wàn)小時(shí)記錄時(shí)LSTM的效果經(jīng)常不如調(diào)好參的Prophet而且訓(xùn)練時(shí)間長(zhǎng)了兩個(gè)數(shù)量級(jí)。如果一定要用LSTM我的建議是控制輸入窗口長(zhǎng)度并加上注意力機(jī)制或者做多步滾動(dòng)預(yù)測(cè)。下面是項(xiàng)目里用TensorFlow/Keras實(shí)現(xiàn)的簡(jiǎn)單LSTM結(jié)構(gòu)from tensorflow.keras.models import Sequential from tensorflow.keras.layers import LSTM, Dense model Sequential() model.add(LSTM(64, return_sequencesTrue, input_shape(24, n_features))) model.add(LSTM(32)) model.add(Dense(1)) model.compile(optimizeradam, lossmse)這里input_shape的24表示用過(guò)去24小時(shí)數(shù)據(jù)作為特征。最開(kāi)始我窗口設(shè)為168一周結(jié)果模型參數(shù)量暴增訓(xùn)練很慢精度也沒(méi)提升。后來(lái)用特征重要性分析發(fā)現(xiàn)過(guò)去24小時(shí)已經(jīng)包含了大部分有效信息。LSTM還有一個(gè)坑尺度敏感。污染物濃度分布右偏直接用原始值訓(xùn)練模型會(huì)對(duì)高濃度區(qū)間擬合不足。我做了Box-Cox變換后再訓(xùn)練預(yù)測(cè)值再逆變換回來(lái)RMSE降低了大約15%。4.4 多模型加權(quán)集成穩(wěn)健性的秘密單一模型總有各自的盲區(qū)。SARIMA擅長(zhǎng)捕捉線性趨勢(shì)和周期Prophet擅長(zhǎng)多季節(jié)分解LSTM擅長(zhǎng)非線性模式。最終軟件里我做了個(gè)簡(jiǎn)單的加權(quán)集成器根據(jù)驗(yàn)證集誤差動(dòng)態(tài)分配權(quán)重。import numpy as np from sklearn.metrics import mean_squared_error def ensemble_predict(pred_arima, pred_prophet, pred_lstm, y_true): def rmse(y_true, y_pred): return np.sqrt(mean_squared_error(y_true, y_pred)) e1 rmse(y_true, pred_arima) e2 rmse(y_true, pred_prophet) e3 rmse(y_true, pred_lstm) total 1/e1 1/e2 1/e3 w1, w2, w3 (1/e1)/total, (1/e2)/total, (1/e3)/total return w1 * pred_arima w2 * pred_prophet w3 * pred_lstm權(quán)重是根據(jù)RMSE的倒數(shù)算的誤差越小的模型權(quán)重越大。這個(gè)集成方法簡(jiǎn)單有效而且每個(gè)月滾動(dòng)重新計(jì)算一次權(quán)重模型自己會(huì)適應(yīng)季節(jié)變化。實(shí)測(cè)集成的RMSE比最好的單一模型下降約8%~12%對(duì)于空氣污染預(yù)測(cè)這種高波動(dòng)場(chǎng)景已經(jīng)很可觀的提升了。5. 源碼工程怎么組織模塊劃分、接口設(shè)計(jì)和文檔配套5.1 工程目錄結(jié)構(gòu)說(shuō)明拿到源碼后第一反應(yīng)就是看目錄結(jié)構(gòu)。好的工程結(jié)構(gòu)應(yīng)該讓人一眼就知道哪里改配置、哪里加模型、哪里看結(jié)果。我最終的項(xiàng)目文件組織結(jié)構(gòu)如下air_pollution_forecast/ # 本博客元信息用注釋括起來(lái)實(shí)際不用 ├── data/ # 原始數(shù)據(jù)和清洗后的數(shù)據(jù) │ ├── raw/ │ └── processed/ ├── src/ # 核心源碼 │ ├── data_fetch.py # 數(shù)據(jù)抓取 │ ├── preprocessing.py # 清洗、插值、異常值處理 │ ├── features.py # 特征工程 │ ├── models/ │ │ ├── sarima_model.py │ │ ├── prophet_model.py │ │ └── lstm_model.py │ ├── ensemble.py # 集成預(yù)測(cè) │ ├── evaluation.py # 評(píng)估指標(biāo) │ └── forecast.py # 預(yù)測(cè)主流程 ├── docs/ # 配套文檔 │ ├── 使用說(shuō)明.md │ ├── 開(kāi)發(fā)文檔.md │ └── API接口文檔.md ├── config.yaml # 配置文件 ├── main.py # 程序入口 └── requirements.txt這里特別強(qiáng)調(diào)config.yaml的作用。剛開(kāi)始我圖省事把API key、數(shù)據(jù)路徑、模型參數(shù)全寫在代碼里后來(lái)?yè)Q數(shù)據(jù)集或者換機(jī)器改代碼改到懷疑人生。配置文件把所有可變參數(shù)提取出來(lái)要改時(shí)只動(dòng)yaml文件代碼零改動(dòng)。5.2 核心模塊的關(guān)鍵接口設(shè)計(jì)我設(shè)計(jì)接口的時(shí)候遵循一個(gè)原則每個(gè)模塊只做一件事模塊之間通過(guò)簡(jiǎn)單的數(shù)據(jù)對(duì)象DataFrame或numpy數(shù)組通信。數(shù)據(jù)獲取模塊對(duì)外只暴露一個(gè)函數(shù)def fetch_hourly_data(city: str, start_date: str, end_date: str) - pd.DataFrame: Returns: DataFrame with columns: [time, pm25, pm10, no2, o3, ...] 預(yù)處理模塊def preprocess(raw_df: pd.DataFrame, config: dict) - pd.DataFrame: 完成缺失值插值、異常值替換、時(shí)間索引對(duì)齊。 Returns: 清洗后的DataFrameindex為DatetimeIndex 模型模塊統(tǒng)一繼承一個(gè)預(yù)測(cè)基類class BaseModel: def fit(self, train_df: pd.DataFrame) - None: ... def predict(self, steps: int) - np.ndarray: ...這樣做的好處是未來(lái)加新模型比如Informer、N-BEATS時(shí)不需要改動(dòng)主流程只要新模型繼承BaseModel并實(shí)現(xiàn)fit和predict就行。5.3 文檔配套不只是給用戶看更是給未來(lái)的自己源碼工程里的文檔質(zhì)量往往被忽視。好的項(xiàng)目文檔至少要有三種使用說(shuō)明、開(kāi)發(fā)文檔、接口文檔。使用說(shuō)明面向“不會(huì)讀代碼的人”重點(diǎn)寫怎么安裝依賴、修改配置、運(yùn)行、看結(jié)果。開(kāi)發(fā)文檔面向“要擴(kuò)展的人”講清楚模塊之間的關(guān)系、數(shù)據(jù)流的方向以及每個(gè)模塊內(nèi)取舍的原因。接口文檔用表格列出每個(gè)函數(shù)的輸入輸出方便調(diào)用時(shí)快速查。我寫文檔時(shí)給自己立了個(gè)規(guī)矩如果我在代碼里寫了“# TODO”或者踩了一些坑一定要在對(duì)應(yīng)的文檔里寫清楚原因。比如LSTM輸入數(shù)據(jù)需要做Box-Cox變換這件事如果不寫進(jìn)文檔三個(gè)月后任何人接手都只會(huì)一頭霧水。5.4 運(yùn)行入口與配置示例軟件主入口main.py的邏輯非常簡(jiǎn)單一共就五步讀取配置、加載數(shù)據(jù)、訓(xùn)練模型、預(yù)測(cè)、輸出結(jié)果和圖表。import yaml from src.data_fetch import fetch_hourly_data from src.preprocessing import preprocess from src.models import sarima_model, prophet_model, lstm_model from src.ensemble import ensemble_predict from src.evaluation import evaluate from src.forecast import plot_forecast if __name__ __main__: with open(config.yaml, r, encodingutf-8) as f: config yaml.safe_load(f) raw fetch_hourly_data( cityconfig[city], start_dateconfig[start_date], end_dateconfig[end_date] ) data preprocess(raw, config) train data.loc[:config[split_date]] test data.loc[config[split_date]:] sarima sarima_model.SARIMAPredictor(config) prophet prophet_model.ProphetPredictor(config) lstm lstm_model.LSTMPredictor(config) for model in [sarima, prophet, lstm]: model.fit(train) preds { sarima: sarima.predict(len(test)), prophet: prophet.predict(len(test)), lstm: lstm.predict(len(test)) } final ensemble_predict(preds[sarima], preds[prophet], preds[lstm], test[pm25]) metrics evaluate(test[pm25], final) print(metrics) plot_forecast(test[pm25], final, save_pathoutput/forecast.png)配置文件示例city: Beijing start_date: 2023-01-01 end_date: 2024-01-01 split_date: 2023-10-01 pollutants: - pm25 - pm10 - no2 model_parameters: sarima: seasonal_period: 24 p_max: 5 q_max: 5 prophet: changepoint_prior_scale: 0.05 lstm: window_size: 24 n_epochs: 50 batch_size: 64配置文件的好處是不同城市、不同污染物只需修改參數(shù)不必動(dòng)代碼。6. 實(shí)測(cè)效果與踩坑記錄那些文檔里不會(huì)寫的事6.1 典型預(yù)測(cè)效果對(duì)比以北京某站點(diǎn)PM2.5小時(shí)濃度為例使用2023年1月到10月訓(xùn)練11月作為測(cè)試集。三種模型的RMSE對(duì)比模型RMSEMAE訓(xùn)練耗時(shí)SARIMA26.818.212分鐘Prophet23.516.46分鐘LSTM21.915.735分鐘集成19.713.9約50分鐘集成模型的RMSE比最好的單一模型降低了10%但訓(xùn)練時(shí)間也最長(zhǎng)。如果你追求快速驗(yàn)證建議先用Prophet如果要上線集成是更穩(wěn)妥的選擇。上面數(shù)據(jù)僅供參考不同城市差別很大污染源復(fù)雜的地方LSTM優(yōu)勢(shì)更明顯。6.2 坑一跨天預(yù)測(cè)時(shí)Prophet的節(jié)假日效應(yīng)Prophet默認(rèn)會(huì)在節(jié)假日改變預(yù)測(cè)但空氣污染不存在“節(jié)假日濃度一定會(huì)降”這種規(guī)律。有一次國(guó)慶節(jié)期間模型預(yù)測(cè)的PM2.5明顯偏低因?yàn)镻rophet把節(jié)假日當(dāng)作特殊日期處理。解決辦法是在訓(xùn)練時(shí)把目標(biāo)污染物換成氣象條件或者干脆關(guān)閉節(jié)假日效應(yīng)只保留日期周期。model Prophet( yearly_seasonalityTrue, weekly_seasonalityTrue, daily_seasonalityTrue, holidaysNone # 關(guān)鍵關(guān)掉節(jié)假日 )6.3 坑二LSTM滾動(dòng)預(yù)測(cè)的誤差累積LSTM多步預(yù)測(cè)時(shí)有一種常見(jiàn)做法把預(yù)測(cè)值當(dāng)作輸入繼續(xù)預(yù)測(cè)下一個(gè)點(diǎn)。這個(gè)做法在空氣質(zhì)量上非常不靠譜。因?yàn)轭A(yù)測(cè)本身有誤差誤差作為輸入回傳給模型會(huì)越滾越大往往到了第12小時(shí)預(yù)測(cè)值就趨向于一條直線。后來(lái)我改成“滾動(dòng)窗口重填歷史值”的方式每預(yù)測(cè)一步就把真實(shí)值或者實(shí)測(cè)數(shù)據(jù)滑入窗口誤差累積一下小了非常多。如果你必須做純滾動(dòng)預(yù)測(cè)沒(méi)有實(shí)時(shí)數(shù)據(jù)回傳最好限制預(yù)測(cè)長(zhǎng)度超過(guò)24小時(shí)一定要給出error bar不能只給一條均值預(yù)測(cè)線。6.4 坑三模型評(píng)估指標(biāo)不能只看RMSERMSE對(duì)高濃度點(diǎn)特別敏感。秋冬季一次嚴(yán)重的霧霾事件可能把RMSE拉高40%但模型對(duì)普通天的預(yù)測(cè)其實(shí)很準(zhǔn)。我建議同時(shí)報(bào)告RMSE、MAE和R^2再分組看不同濃度區(qū)間的誤差比如“0-100”、“100-200”、“200”三個(gè)區(qū)間分別算指標(biāo)。如果高濃度區(qū)間誤差大說(shuō)明模型對(duì)污染過(guò)程的峰值預(yù)測(cè)能力弱這時(shí)需要加強(qiáng)特征或者換非線性模型。還可以補(bǔ)充一個(gè)面向預(yù)警的指標(biāo)在PM2.5超過(guò)150這個(gè)閾值時(shí)模型預(yù)測(cè)的命中率hit rate和漏報(bào)率miss rate。因?yàn)榄h(huán)境監(jiān)測(cè)任務(wù)真正關(guān)心的是“重度污染是否會(huì)被預(yù)警”而不是均方誤差好不好看。6.5 坑四時(shí)間索引的時(shí)區(qū)問(wèn)題做項(xiàng)目時(shí)還遇到過(guò)時(shí)區(qū)坑。數(shù)據(jù)源返回的時(shí)間是UTC本地分析要用北京時(shí)間UTC8。如果直接用原始時(shí)間戳建索引所有字段都會(huì)偏移8小時(shí)日周期建模就會(huì)錯(cuò)位晚上8點(diǎn)的濃度峰值被當(dāng)成中午12點(diǎn)。處理方式是在數(shù)據(jù)抓取階段統(tǒng)一轉(zhuǎn)成目標(biāo)時(shí)區(qū)df[time] pd.to_datetime(df[time], utcTrue) df[time] df[time].dt.tz_convert(Asia/Shanghai)后續(xù)所有數(shù)據(jù)都基于本地時(shí)間避免后期反復(fù)切換。這個(gè)坑新手特別容易踩因?yàn)榭雌饋?lái)數(shù)據(jù)數(shù)量都一樣畫圖也沒(méi)問(wèn)題但模型效果就是差本質(zhì)是時(shí)序錯(cuò)位。7. 從軟件到可交付項(xiàng)目打包與使用體驗(yàn)優(yōu)化7.1 requirements.txt與Python環(huán)境項(xiàng)目要交到別人手里最尷尬的是對(duì)方裝不了依賴。我強(qiáng)烈建議你固定Python版本3.9以上3.11以下并且在requirements.txt里寫清楚具體版本不要只寫包名。當(dāng)時(shí)我吃過(guò)虧直接依賴最新版Prophet對(duì)方舊服務(wù)器上裝不上折騰了三天。推薦的最小依賴清單pandas2.0.3 numpy1.24.3 statsmodels0.14.0 pmdarima2.0.4 prophet1.1.4 scikit-learn1.3.0 tensorflow2.13.0 PyYAML6.0 matplotlib3.7.2 requests2.31.0建議再寫一個(gè)setup.md里面記錄你在什么系統(tǒng)上測(cè)試通過(guò)、每一步的安裝命令?，F(xiàn)在很多人已經(jīng)用uv或者poetry但無(wú)論用什么核心依賴鎖文件一定要有否則交付就是災(zāi)難。7.2 可視化輸出讓預(yù)測(cè)結(jié)果一目了然預(yù)測(cè)軟件除了報(bào)數(shù)字一定要出圖。我實(shí)現(xiàn)的圖表包含三部分歷史濃度曲線、預(yù)測(cè)均值曲線、90%置信區(qū)間帶。代碼不多但效果很好。import matplotlib.pyplot as plt fig, ax plt.subplots(figsize(12, 5)) ax.plot(train[pm25], colorgray, labeltrain) ax.plot(test[pm25], colorblue, labelactual) ax.plot(pred_index, pred_mean, colorred, labelforecast) ax.fill_between( pred_index, pred_lower, pred_upper, colorred, alpha0.2, label90% interval ) ax.legend() ax.set_ylabel(PM2.5 (μg/m3)) fig.savefig(output/forecast.png, dpi150, bbox_inchestight)置信區(qū)間帶非常重要。沒(méi)有區(qū)間帶決策者很容易把均值預(yù)測(cè)當(dāng)成確定值。環(huán)境污染預(yù)警需要風(fēng)險(xiǎn)意識(shí)區(qū)間寬度就是風(fēng)險(xiǎn)提示。SARIMA和Prophet都自帶預(yù)測(cè)區(qū)間LSTM需要自己從多個(gè)隨機(jī)初始化模型的輸出中估算。7.3 與監(jiān)測(cè)平臺(tái)對(duì)接的增量式預(yù)測(cè)跨過(guò)一個(gè)里程碑后我把軟件改成了增量式預(yù)測(cè)每天凌晨2點(diǎn)重新拉取過(guò)去24小時(shí)真實(shí)數(shù)據(jù)重新訓(xùn)練或者部分?jǐn)M合模型然后生成未來(lái)3天的預(yù)測(cè)。增量式的好處是模型能及時(shí)響應(yīng)最近幾天排放源的突然變化比如工廠臨時(shí)檢修導(dǎo)致的濃度異常。增量式要注意窗口尺寸。如果每次訓(xùn)練只用了最近一個(gè)月的數(shù)據(jù)模型會(huì)忘記冬季和夏季差異進(jìn)入2月之后之前學(xué)到的1月規(guī)律完全沒(méi)用。我的做法是“長(zhǎng)期數(shù)據(jù)全量訓(xùn)練近期數(shù)據(jù)加權(quán)微調(diào)”比如每年年初做一次全量重訓(xùn)平時(shí)每日更新時(shí)用最近30天數(shù)據(jù)微調(diào)。這個(gè)策略聽(tīng)起來(lái)復(fù)雜但代碼實(shí)現(xiàn)并不難核心就是在原有模型實(shí)例上調(diào)用fit傳入增量數(shù)據(jù)即可Prophet和SARIMA都支持繼續(xù)訓(xùn)練。只有LSTM需要全量重訓(xùn)所以實(shí)際生產(chǎn)中LSTM通常作為離線備選模型線上主力還是ProphetSARIMA。8. 再聊兩句這套軟件還能怎么擴(kuò)展最后基于我的實(shí)際經(jīng)驗(yàn)給拿到源碼準(zhǔn)備二次開(kāi)發(fā)的人幾個(gè)方向參考。第一把單站點(diǎn)預(yù)測(cè)擴(kuò)展成區(qū)域預(yù)測(cè)?？諝馕廴静皇菃吸c(diǎn)問(wèn)題把周邊多個(gè)站點(diǎn)的數(shù)據(jù)和氣象場(chǎng)作為共同輸入可以顯著提升預(yù)測(cè)精度尤其是靜穩(wěn)天氣下污染傳輸過(guò)程。第二接入實(shí)時(shí)監(jiān)測(cè)設(shè)備的自動(dòng)校準(zhǔn)。設(shè)備漂移會(huì)導(dǎo)致數(shù)據(jù)出現(xiàn)系統(tǒng)性偏差時(shí)間序列模型會(huì)被帶偏如果能在預(yù)處理階段加入卡爾曼濾波或者設(shè)備偏差校正模型效果會(huì)更穩(wěn)。第三把預(yù)測(cè)結(jié)果通過(guò)Web服務(wù)暴露出來(lái)。我后期用FastAPI封裝了預(yù)測(cè)接口讓監(jiān)測(cè)站同事可以通過(guò)網(wǎng)頁(yè)輸入城市和時(shí)間范圍后臺(tái)自動(dòng)跑模型返回JSON格式的預(yù)測(cè)值和區(qū)間。這個(gè)不難源碼的forecast.py已經(jīng)提供了核心函數(shù)包一層HTTP即可。這套“基于Python的時(shí)間序列分析的大氣污染預(yù)測(cè)軟件”的源碼和文檔目前已經(jīng)整理成完整工程。讀完全文你會(huì)發(fā)現(xiàn)真正困難的部分從來(lái)不是跑通一個(gè)模型而是把數(shù)據(jù)、模型、工程、文檔全部串成一個(gè)可持續(xù)使用的東西。希望這篇分享能幫你少走一些彎路在做自己的空氣污染預(yù)測(cè)項(xiàng)目時(shí)把精力花在模型調(diào)優(yōu)和業(yè)務(wù)落地這些真正有價(jià)值的事情上。