簡介這份資料包面向氣象與水文領(lǐng)域的研究人員及數(shù)據(jù)分析學(xué)習(xí)者聚焦最大風(fēng)速序列的均一化訂正問題幫助消除因儀器更換、測量方法調(diào)整或站點遷移帶來的系統(tǒng)性偏差使不同站點、不同時期的風(fēng)速記錄具備可比性。包內(nèi)共3個文件包含1個Python腳本、1個CSV數(shù)據(jù)文件和1份PDF說明文檔壓縮包約377KB分別對應(yīng)訂正算法實現(xiàn)、示例風(fēng)速數(shù)據(jù)與代碼使用說明便于直接運行與對照理解。目前已有591人學(xué)習(xí)下載。通過腳本與數(shù)據(jù)讀者可完整走通數(shù)據(jù)讀取、缺失值與異常值預(yù)處理、訂正因子計算、均一化序列生成及結(jié)果驗證等環(huán)節(jié)并借助可視化對比訂正前后的風(fēng)速變化趨勢掌握Thorne-Wyatt、Renfrew、HOM等方法的實現(xiàn)思路為氣候研究中的非均一數(shù)據(jù)處理提供可復(fù)用的實踐參考。1. 風(fēng)速均一化訂正為什么最大風(fēng)速序列不能直接拿來做趨勢分析如果你手上有某個氣象站 1960 年至今的日最大風(fēng)速序列直接畫一條趨勢線大概率會得到一個“風(fēng)速在顯著下降”的結(jié)論。但這個結(jié)論很可能是假的。原因不在天氣而在觀測系統(tǒng)本身測風(fēng)儀器的型號換過、塔高變過、周邊建筑物長起來了、觀測時次從 4 次變成 24 次、甚至站址搬過。這些非氣候因素造成的跳變會淹沒真實的氣候信號。風(fēng)速均一化訂正要解決的就是把這些“人為斷點”識別出來并訂正掉讓訂正后的序列只保留氣候意義上的變化。這篇筆記以日最大風(fēng)速為例把均一化訂正的完整鏈路走一遍從數(shù)據(jù)準備、斷點檢測、到訂正量計算和效果驗證每一步都給可復(fù)現(xiàn)的代碼和參數(shù)說明。適合手里有長序列風(fēng)速數(shù)據(jù)、準備做趨勢分析或極值統(tǒng)計的氣象水文從業(yè)者也適合剛接觸均一化、想找一個完整案例上手的人。需要先說清楚一個邊界均一化訂正不是“把數(shù)據(jù)修好看”它是有假設(shè)的。核心假設(shè)是——參考序列鄰近站或再分析與被檢站經(jīng)歷相同的氣候變化但參考序列本身是均一的。如果參考站也換過儀器那訂正結(jié)果就是錯的。所以選參考站這件事比后面跑什么檢測算法都重要。我一般會花一半時間在選站和元數(shù)據(jù)核對上剩下的一半才交給算法。2. 均一化訂正的技術(shù)路線從元數(shù)據(jù)到統(tǒng)計檢驗怎么選2.1 先搞清楚斷點從哪來元數(shù)據(jù)優(yōu)先統(tǒng)計檢驗兜底風(fēng)速序列的非均一性來源按影響從大到小排通常是這幾類儀器更換尤其從風(fēng)杯換成超聲啟動風(fēng)速閾值變了、觀測高度變化10m 換到 2m 或反之、站址遷移、周邊環(huán)境變化新建高樓、樹木生長、觀測時次和統(tǒng)計方法變化日最大風(fēng)速的統(tǒng)計窗口變了。這些信息最可靠的來源是臺站元數(shù)據(jù)——歷史沿革表、儀器更換記錄、站址變更文件。但現(xiàn)實是很多老站的元數(shù)據(jù)殘缺甚至根本沒有。這時候才輪到統(tǒng)計檢驗上場。統(tǒng)計檢驗的邏輯是如果序列在某點前后相對于參考序列的差值發(fā)生了顯著變化那這個點就可能是斷點。注意是“相對于參考序列”不是看序列本身的均值跳變——因為氣候變化本身也會造成均值變化只有“差值”的變化才指向非均一性。常見做法是元數(shù)據(jù)和統(tǒng)計檢驗結(jié)合先用元數(shù)據(jù)圈定可疑年份再用統(tǒng)計檢驗確認或者反過來統(tǒng)計檢驗找出候選斷點再回元數(shù)據(jù)里找解釋。兩者對上了斷點可信度就高。2.2 四種主流方法的選擇SNHT、PMT、MASH、RHtest 各適合什么場景風(fēng)速均一化訂正常用的方法有這么幾類選哪個取決于你的數(shù)據(jù)條件和目標方法全稱/來源適合場景局限SNHT標準正態(tài)均一性檢驗單斷點、序列較長、參考序列質(zhì)量好多斷點場景容易漏檢PMTPairwise Multiple Test多斷點、有多個參考站計算量大參考站需均一MASHMultiple Analysis of Series for Homogenization多站聯(lián)合檢測適合臺站網(wǎng)實現(xiàn)復(fù)雜對參考站數(shù)量有要求RHtestR 語言 rhTest 包單站、可結(jié)合元數(shù)據(jù)、支持多種檢驗需要 R 環(huán)境參數(shù)較多我一般會這樣選如果只有一個站、參考序列是再分析資料用 RHtest 的 PMT 或 SNHT 模式如果是區(qū)域臺站網(wǎng)、有多個鄰近站用 MASH 或 PMT 做多站聯(lián)合檢測。風(fēng)速這個變量比氣溫噪聲大單站檢測的漏檢率不低所以有條件盡量用多站方法。2.3 參考序列怎么選鄰近站、再分析、還是區(qū)域平均參考序列的質(zhì)量直接決定訂正成敗。三種常見選擇鄰近站選距離近、海拔相近、氣候背景一致的站。要求參考站本身經(jīng)過均一化檢驗或者至少有完整元數(shù)據(jù)證明沒換過儀器。距離一般控制在 50km 以內(nèi)山區(qū)要更近。再分析資料比如 ERA5 的 10m 風(fēng)速。優(yōu)點是時空連續(xù)、沒有斷點缺點是分辨率粗對局地風(fēng)速的刻畫能力有限尤其是復(fù)雜地形。用再分析做參考時通常要先做尺度匹配——把再分析風(fēng)速插值到站點或者用回歸建立兩者關(guān)系。區(qū)域平均把多個鄰近站平均成一個區(qū)域序列。平均能削弱單站噪聲但如果區(qū)域內(nèi)有站也非均一會污染參考序列。我的習(xí)慣是優(yōu)先用經(jīng)過檢驗的鄰近站沒有就用再分析兩者都有就做交叉驗證——用鄰近站訂正一遍用再分析訂正一遍看結(jié)果差多少。差得大說明參考序列本身有問題得回頭查。3. 用 Python 跑通最大風(fēng)速均一化訂正的最小流程3.1 數(shù)據(jù)準備日最大風(fēng)速序列的讀取與質(zhì)控假設(shè)你拿到的是 CSV 格式的日最大風(fēng)速數(shù)據(jù)列包括日期、風(fēng)速值、站點號。第一步不是直接跑檢測而是質(zhì)控。風(fēng)速數(shù)據(jù)的常見問題缺測、異常大值比如臺風(fēng)天記錄到 50m/s 但實際是儀器故障、連續(xù)相同值儀器卡死、單位不統(tǒng)一m/s 和 km/h 混用。import pandas as pd import numpy as np # 讀取日最大風(fēng)速數(shù)據(jù) df pd.read_csv(daily_max_wind.csv, parse_dates[date]) df df.sort_values(date).reset_index(dropTrue) # 基本質(zhì)控 # 1. 風(fēng)速為負或超過合理上限的標記為缺測 df.loc[(df[wind_speed] 0) | (df[wind_speed] 60), wind_speed] np.nan # 2. 連續(xù)相同值超過5天的標記為可疑儀器卡死 same_count (df[wind_speed] df[wind_speed].shift()).astype(int) same_group (same_count 0).cumsum() group_sizes df.groupby(same_group)[wind_speed].transform(size) df.loc[group_sizes 5, wind_speed] np.nan # 3. 統(tǒng)計缺測率 missing_rate df[wind_speed].isna().mean() print(f缺測率: {missing_rate:.2%}) # 4. 按月統(tǒng)計檢查是否有整月缺測 monthly_count df.set_index(date)[wind_speed].resample(M).count() print(monthly_count[monthly_count 20])這段代碼做了四件事把物理上不可能的值置為缺測、識別儀器卡死造成的連續(xù)相同值、統(tǒng)計整體缺測率、檢查是否有整月缺測。參數(shù)說明風(fēng)速上限 60m/s 是通用閾值沿海臺風(fēng)影響區(qū)可以放寬到 70連續(xù)相同值閾值 5 天是經(jīng)驗值干燥地區(qū)可以放寬到 7 天。缺測率超過 20% 的年份后續(xù)訂正要謹慎因為斷點檢測對缺測敏感。3.2 斷點檢測用 RHtest 的 SNHT 模式找候選斷點RHtest 是加拿大環(huán)境部開發(fā)的均一化檢驗工具有 R 包和命令行版本。這里用 Python 調(diào)用 R 的方式演示也可以直接用 R 跑。核心是rhTest函數(shù)支持 SNHT、PMT 等多種檢驗。import subprocess import os # 準備 RHtest 輸入文件兩列第一列參考序列第二列待檢序列 # 這里假設(shè)參考序列是鄰近站已經(jīng)過質(zhì)控 ref pd.read_csv(ref_station.csv, parse_dates[date]).set_index(date)[wind_speed] target df.set_index(date)[wind_speed] # 對齊日期 combined pd.concat([ref, target], axis1, joininner).dropna() combined.columns [ref, target] combined.to_csv(rhtest_input.txt, sep , indexFalse, headerFalse) # 調(diào)用 RHtest假設(shè)已安裝 R 和 rhTest 包 r_script library(rhTest) data - read.table(rhtest_input.txt) # SNHT 檢驗返回候選斷點 result - rhTest(data$V2, data$V1, testSNHT) print(result$breaks) write.csv(result$breaks, breaks.csv, row.namesFALSE) with open(run_rhtest.R, w) as f: f.write(r_script) subprocess.run([Rscript, run_rhtest.R], checkTrue) breaks pd.read_csv(breaks.csv) print(breaks)邏輯說明RHtest 需要兩列輸入——參考序列和待檢序列按日期對齊。rhTest函數(shù)返回候選斷點及其顯著性。參數(shù)說明testSNHT指定用標準正態(tài)均一性檢驗適合單斷點如果懷疑多斷點改成testPMT。顯著性水平默認 0.05可以通過p.value參數(shù)調(diào)整。注意RHtest 對缺測敏感輸入前要確保兩列都沒有缺測或者用插值補齊——但插值本身會引入誤差缺測多的年份建議單獨處理。3.3 訂正量計算用分位數(shù)匹配做逐日訂正找到斷點后下一步是計算訂正量。風(fēng)速的訂正不能簡單加減均值差因為風(fēng)速分布是偏態(tài)的不同分位數(shù)的偏差不一樣。常用做法是分位數(shù)匹配Quantile Matching在斷點前后各取一段窗口分別計算參考序列和待檢序列的分位數(shù)關(guān)系然后把這個關(guān)系外推到整個時段。def quantile_matching_correction(target, ref, break_date, window5): 用分位數(shù)匹配計算訂正量 target: 待檢序列Series索引為日期 ref: 參考序列Series索引為日期 break_date: 斷點日期 window: 斷點前后各取多少年做擬合 break_year pd.Timestamp(break_date).year # 斷點前窗口 pre_mask (target.index.year break_year - window) (target.index.year break_year) # 斷點后窗口 post_mask (target.index.year break_year) (target.index.year break_year window) # 計算斷點前后的分位數(shù)關(guān)系 quantiles np.arange(0.05, 1.0, 0.05) pre_target_q target[pre_mask].quantile(quantiles) pre_ref_q ref[pre_mask].quantile(quantiles) post_target_q target[post_mask].quantile(quantiles) post_ref_q ref[post_mask].quantile(quantiles) # 斷點前的比值關(guān)系target/ref pre_ratio pre_target_q / pre_ref_q post_ratio post_target_q / post_ref_q # 訂正系數(shù)把斷點后的關(guān)系調(diào)整到斷點前 correction_factor pre_ratio / post_ratio # 對斷點后的數(shù)據(jù)做逐日訂正 corrected target.copy() post_all target.index pd.Timestamp(break_date) # 用插值把分位數(shù)訂正系數(shù)映射到每個值 for i, val in target[post_all].items(): # 找到 val 在斷點后分位數(shù)中的位置 q_rank (post_target_q val).mean() # 插值得到對應(yīng)的訂正系數(shù) factor np.interp(q_rank, quantiles, correction_factor) corrected.loc[i] val * factor return corrected # 應(yīng)用訂正 corrected_series quantile_matching_correction(target, ref, 1985-01-01, window5)邏輯說明這段代碼的核心思想是——如果斷點前后 target 和 ref 的比值關(guān)系發(fā)生了變化那這個變化就是非均一性造成的需要把斷點后的比值調(diào)整回斷點前的水平。參數(shù)說明window5表示用斷點前后各 5 年做擬合窗口太短擬合不穩(wěn)太長會混入氣候變化信號一般 5 到 10 年。quantiles從 0.05 到 0.95步長 0.05覆蓋了風(fēng)速的主要分布范圍。注意分位數(shù)匹配假設(shè)斷點前后的偏差是乘性的如果偏差是加性的應(yīng)該用差值而不是比值。風(fēng)速一般用乘性更合理因為高風(fēng)速的絕對偏差通常更大。3.4 效果驗證訂正前后序列對比與趨勢檢驗訂正完不能直接信要做驗證。驗證分兩步一是看訂正后斷點是否消除二是看訂正對趨勢的影響。import matplotlib.pyplot as plt from scipy import stats # 1. 訂正前后對比圖 fig, axes plt.subplots(2, 1, figsize(12, 8)) axes[0].plot(target.index, target.values, label原始序列, alpha0.7) axes[0].plot(corrected_series.index, corrected_series.values, label訂正后, alpha0.7) axes[0].axvline(pd.Timestamp(1985-01-01), colorr, linestyle--, label斷點) axes[0].legend() axes[0].set_ylabel(日最大風(fēng)速 (m/s)) # 2. 訂正前后趨勢對比 def trend_test(series): 用 Mann-Kendall 檢驗計算趨勢 n len(series) s 0 for i in range(n-1): for j in range(i1, n): s np.sign(series.iloc[j] - series.iloc[i]) var_s n*(n-1)*(2*n5)/18 z (s - np.sign(s)) / np.sqrt(var_s) if s ! 0 else 0 p 2 * (1 - stats.norm.cdf(abs(z))) # Sens slope slopes [] for i in range(n-1): for j in range(i1, n): slopes.append((series.iloc[j] - series.iloc[i]) / (j - i)) slope np.median(slopes) return slope, p slope_orig, p_orig trend_test(target.dropna()) slope_corr, p_corr trend_test(corrected_series.dropna()) print(f原始序列趨勢: {slope_orig:.4f} m/s/年, p{p_orig:.4f}) print(f訂正后趨勢: {slope_corr:.4f} m/s/年, p{p_corr:.4f})邏輯說明Mann-Kendall 是非參數(shù)趨勢檢驗不要求數(shù)據(jù)正態(tài)分布適合風(fēng)速這種偏態(tài)變量。Sens slope 給出趨勢的穩(wěn)健估計。參數(shù)說明p 值小于 0.05 認為趨勢顯著。如果訂正前后趨勢方向或顯著性發(fā)生根本變化說明訂正起了作用——但也要警惕訂正過度。我一般會對比訂正前后的趨勢如果訂正后趨勢從顯著變不顯著或者斜率變化超過 50%就要回頭檢查斷點是否找多了。4. 風(fēng)速均一化訂正的避坑清單五個血淚教訓(xùn)4.1 坑一參考站自己就不均一訂正結(jié)果全歪現(xiàn)象訂正后序列在斷點處確實平滑了但整體趨勢變得很奇怪和鄰近區(qū)域其他站對不上。原因參考站本身在同期也換過儀器或遷過站但沒做均一化檢驗。用非均一的參考去訂正待檢站相當于用一把不準的尺子去量另一把尺子。解決參考站必須先做均一化檢驗。如果參考站也有斷點要么換站要么先訂正參考站。我一般會要求參考站的元數(shù)據(jù)完整或者至少用 RHtest 跑一遍確認沒有顯著斷點。4.2 坑二斷點檢測對缺測太敏感缺測多的年份誤報現(xiàn)象檢測出一堆斷點但元數(shù)據(jù)里這些年份沒有任何儀器或站址變化。原因缺測導(dǎo)致序列的統(tǒng)計特性在缺測前后發(fā)生變化被算法誤判為斷點。風(fēng)速數(shù)據(jù)在早期1960-1980缺測率普遍偏高。解決檢測前先統(tǒng)計逐月缺測率缺測率超過 30% 的年份標記出來檢測時排除或單獨處理。RHtest 有處理缺測的選項但效果有限。我的做法是缺測率高的時段不參與斷點檢測但訂正時仍然處理——用鄰近時段的關(guān)系外推。4.3 坑三訂正窗口選太長把氣候信號也訂正掉了現(xiàn)象訂正后序列的趨勢比原始序列還大或者趨勢方向反轉(zhuǎn)。原因分位數(shù)匹配的窗口如果取 10 年以上窗口內(nèi)本身包含氣候變化趨勢導(dǎo)致訂正系數(shù)里混入了氣候信號訂正時把真實趨勢也放大了。解決窗口一般取 5 年最多不超過 8 年。如果斷點密集窗口還要縮短。另外訂正前先對序列做去趨勢處理訂正后再把趨勢加回去——這個做法更穩(wěn)妥但實現(xiàn)復(fù)雜一些。4.4 坑四多斷點場景用單斷點方法漏檢嚴重現(xiàn)象序列明顯有多個跳變但 SNHT 只檢測出一個斷點。原因SNHT 是單斷點檢驗序列有多個斷點時它只能找到最顯著的那個其他斷點被掩蓋。解決用 PMT 或 MASH 做多斷點檢測。如果只能用 SNHT就分段檢測——先找最顯著的斷點把序列分成兩段再在每段里繼續(xù)找直到?jīng)]有顯著斷點。但分段檢測會累積誤差斷點多了要謹慎。4.5 坑五訂正后不做獨立驗證自說自話現(xiàn)象訂正報告里只有訂正前后的對比圖沒有獨立驗證。原因?qū)Ρ葓D只能說明訂正起了作用不能說明訂正對了。訂正可能過度也可能不足。解決至少做兩種獨立驗證。一是用另一套參考序列比如再分析重復(fù)訂正看結(jié)果是否一致二是留出部分時段不參與訂正用訂正后的關(guān)系去預(yù)測看預(yù)測誤差。我一般會用 ERA5 做交叉驗證——如果鄰近站訂正和 ERA5 訂正的結(jié)果差在 10% 以內(nèi)就認為訂正可信。5. 進階技巧用元數(shù)據(jù)約束斷點檢測把誤報率壓下來前面講的流程是純統(tǒng)計的但實際業(yè)務(wù)里元數(shù)據(jù)才是最強的約束。我現(xiàn)在的做法是先用元數(shù)據(jù)圈定“可疑年份”——儀器更換、站址遷移、觀測時次變化的年份然后在這些年份附近做統(tǒng)計檢驗而不是全序列盲掃。這樣誤報率能降一大半。具體操作把元數(shù)據(jù)整理成一張表列包括年份、變更類型、變更描述。然后用 Python 把這張表和斷點檢測結(jié)果做匹配。# 元數(shù)據(jù)表 metadata pd.DataFrame({ year: [1975, 1985, 1998, 2005], change_type: [儀器更換, 站址遷移, 觀測時次變化, 儀器更換], description: [風(fēng)杯換型, 搬遷至新址, 4次改24次, 換超聲風(fēng)速儀] }) # 斷點檢測結(jié)果 detected_breaks pd.DataFrame({ break_date: [1975-06-01, 1985-03-01, 2005-08-01], p_value: [0.01, 0.03, 0.02] }) detected_breaks[year] pd.to_datetime(detected_breaks[break_date]).dt.year # 匹配斷點年份與元數(shù)據(jù)變更年份相差不超過1年 matched detected_breaks.merge(metadata, onyear, howleft) matched[confirmed] matched[change_type].notna() print(matched)邏輯說明這段代碼把統(tǒng)計檢測出的斷點和元數(shù)據(jù)變更記錄做匹配。如果斷點年份附近有元數(shù)據(jù)變更記錄這個斷點就是“有解釋的”可信度高如果沒有就是“無解釋的”需要進一步排查——可能是參考序列的問題也可能是氣候變化造成的假斷點。參數(shù)說明匹配窗口設(shè)為 ±1 年因為元數(shù)據(jù)記錄的年份和實際變更時間可能有偏差。如果元數(shù)據(jù)精確到月窗口可以縮小到 ±3 個月。還有一個技巧對“無解釋的斷點”不要急著訂正先檢查參考序列在同期有沒有異常。如果參考序列在同期也有跳變那這個斷點很可能是參考序列的問題不是待檢站的。這時候應(yīng)該換參考序列而不是硬訂正。最后說一個我自己的習(xí)慣每次訂正完我都會把訂正前后的序列、斷點位置、元數(shù)據(jù)變更記錄畫在一張圖上人工過一遍。算法再先進也不如人眼掃一遍來得可靠。有些斷點算法覺得顯著但你看一眼就知道是臺風(fēng)年造成的假信號——這種就得手動排除。均一化訂正這件事算法是工具判斷在人。希望幫到你。本文還有配套的精品資源點擊獲取