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