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