胡義成,張正陽,王秋香*
(1.中國氣象局烏魯木齊沙漠氣象研究所,新疆 烏魯木齊830002;2.新疆氣象信息中心,新疆 烏魯木齊830002)
目前存檔的探空資料中蘊含了氣候變化的重要信息,但是由于探空資料中人為因素造成的不均一性和氣候因素造成的變化混淆,使得研究者很難利用探空資料準(zhǔn)確地判斷氣候變化趨勢。為了使探空資料很好地運用到氣候研究中,很多學(xué)者致力于檢驗和訂正探空資料中時間序列的不均一性。和地面氣象要素均一性檢驗相類似的是[1-9],對于探空資料不均一性的檢驗需要識別由于人為因素造成的間斷點,估計間斷點偏離大小并對間斷點進(jìn)行訂正[10]。我國高空觀測站建站時間較早,最早出現(xiàn)在1950 年,大部分高空站從1960 年陸續(xù)開始建站,站點遍布全國。均一性研究在20 世紀(jì)90 年代開始進(jìn)行,但只有少數(shù)臺站進(jìn)行了試驗性研究[10]。近期對全國多數(shù)探空溫度序列的均一性分析證實了輻射訂正、儀器換型等方法的修改造成了序列的非均一性,這種非均一性在各站探空資料中不同程度的存在,均一化過程也存在著相當(dāng)?shù)牟淮_定性[11,12]。
喀什高空氣象觀測站(以下簡稱“喀什站”)經(jīng)歷了多次儀器變更、臺站位置變動、觀測規(guī)范變更、觀測時間變化以及輻射訂正方法變化(以下簡稱:輻射訂正)。這些活動是否對高空資料的均一性產(chǎn)生了影響,高空資料是否可以連續(xù)使用,都是值得去研究和探討的問題。本文利用加拿大環(huán)境部氣候研究中心研發(fā)的RHtest[13,14]均一化檢驗方法,對喀什站1958—2016 年各個標(biāo)準(zhǔn)等壓面層的年平均溫度資料進(jìn)行均一性檢驗與訂正,對喀什站高空資料的完整性、連續(xù)性和均一性進(jìn)行了分析,有利于幫助臺站提高單站高空資料的質(zhì)量,也可以為新疆高空氣象觀測業(yè)務(wù)的運行提供參考。
喀什站始建于1954 年7 月,自1957 年6 月1日承擔(dān)著北京時間08、20 時常規(guī)探空和測風(fēng)的同步觀測任務(wù)。建站以來,喀什站址多次變動。1955 年5月1 日因不明原因發(fā)生了站址位置變動(圖1);1957年7 月1 日,站址由疏附縣浩罕莊遷至喀什市西南郊;1964 年1 月1 日又因不明原因?qū)е抡局方?jīng)緯度信息發(fā)生變化;因2014 年1 月1 日觀測環(huán)境已遠(yuǎn)遠(yuǎn)不能滿足氣象觀測的要求。
圖1 歷次遷站過程圖
使用的基礎(chǔ)資料來自國家氣象信息中心“中國高空規(guī)定等壓面定時值數(shù)據(jù)集(V2.1)”及高空臺站元數(shù)據(jù)集。資料時段為喀什站1958—2016 年,包括北京時間08、20 時(世界時00 時和12 時)各規(guī)定等壓面位勢高度、氣溫、露點溫度觀測數(shù)據(jù)和風(fēng)向風(fēng)速。數(shù)據(jù)集中1951—2010 年數(shù)據(jù)來源于高空觀測紙質(zhì)月報表數(shù)字化資料,2010—2016 年12 月的數(shù)據(jù)來源于高空全月觀測數(shù)據(jù)文件(高空月報G 文件)。該數(shù)據(jù)集在研制過程中借鑒了《無線電探空資料質(zhì)量控制(QX/T 123-2011)》中相關(guān)質(zhì)量控制方法及參數(shù),以多源數(shù)據(jù)對比檢測為主,對可疑資料實施人工核查與訂正。不僅更正了紙質(zhì)資料本身存在的觀測和抄錄錯誤數(shù)據(jù),還更正了信息化過程導(dǎo)致的資料錯誤和無數(shù)據(jù)現(xiàn)象,最大程度地確保了研究數(shù)據(jù)的可靠性。
參考序列是ERA-Interim,ERA-Interim 是歐洲中期天氣預(yù)報中心(ECMWF)發(fā)布的最新的全球大氣數(shù)值預(yù)報再分析資料,它是繼早期資料(如ERA-15,ERA-40)后的新產(chǎn)品,目的是對ERA-40 等資料進(jìn)行完善,從而逐漸取代ERA-40。ERA-Interim所采用的同化系統(tǒng)是ECMWF 集成預(yù)報系統(tǒng)(IFS,32R2 版本的全球譜模式)。相對于ERA-40 在很多的關(guān)鍵部分都進(jìn)行了改進(jìn),包括模式的改進(jìn)、四維空間變量分析的運用、對衛(wèi)星資料變量偏差的修正以及對觀測系統(tǒng)中其他數(shù)據(jù)的處理等。
RHtestsV5 軟件包是在RHtestsV4 分位數(shù)匹配訂正(Quantile-Matching adjustments,QM)中增加使用參考序列的基礎(chǔ)上,在序列檢驗中提供了一階滯后自相關(guān),還提供了數(shù)據(jù)序列的均值訂正、QM 訂正和序列曲線圖以及結(jié)果的回歸模擬。可以對日或者小時數(shù)據(jù)序列進(jìn)行QM 訂正。檢驗方法包括兩種,一種是懲罰最大t 檢驗(PMT)方法,其檢驗過程中需要建立參考序列,待檢序列與參考序列的差值(差值序列)是被檢驗的對象;另一種是懲罰最大F 檢驗(PMF),適用于無參考序列的檢驗過程。
PMT 方法檢驗過程如下:
假設(shè){Xt}(t=1,…,N)代表一個ⅡD 正態(tài)分布序列,要檢驗時間序列{Xt}是否存在斷點。
備擇假設(shè):
μ1≠μ2,并且“ { Xt}~‖DN(μ,σ2)”代表{ Xt}遵循正態(tài)分布,其均值為μ,方差為σ2。當(dāng)Ha為真時,t=k 時的點被稱作間斷點,Δ=被稱作平均突變的大小,最可能的間斷點服從以下分布:
其中:P(κ)是建立的經(jīng)驗性的懲罰因子,其建立方法可參考文獻(xiàn)[13]。
通過回歸檢驗算法來檢驗出多個間斷點,首先在t∈{Nmin,Nmin+1,…,(N-Nmin)}時段找出最可能的間斷點C0,Nmin指的是時間序列中被分割片段的最小長度,也就是說,兩個相鄰斷點之間子序列的長度或者第一(最后一)個斷點到序列起點(終點)之間的長度,類似的,分別在時段t={Nmin,…,(C0-Nmin)}和t={(C0+1+Nmin),…,(N-Nmin)}找出次可能的間斷點C01和C03,接下來在C01和C03之間的時段找出最可能的間斷點C02,然后在C01,C02和C03之間尋找最可能的間斷點,即依次分段找出序列中各段最可能的間斷點,計算所有間斷點的統(tǒng)計量PTmax,估計lag-1 自相關(guān)和p值的大小,找出最大的PTmax值對應(yīng)的間斷點,如果是顯著的,該突變就被認(rèn)為是找到的第一個間斷點,該間斷點被記為NC=1。尋找該間斷點位置之后每段最可能的間斷點,估計其顯著性,找出下一個可能的間斷點,重復(fù)該過程,將間斷點按照顯著性由大到小排列,形成間斷點列表,判斷最小的間斷點是否顯著,當(dāng)不顯著時剔除該間斷點,需要再次評估剩余間斷點的顯著性,最終保留都統(tǒng)計顯著的間斷點即為序列檢驗得到的間斷點,受分段檢驗的影響,該方法能夠檢驗的最短序列長度是20,詳細(xì)過程見文獻(xiàn)[14]。
PMF 方法檢驗過程如下:
假設(shè)εt代表變量,該變量的均值為0,方差為σ2,對于存在的線性趨勢β,其時間序列為{Xt},要檢查在t=k 時刻是否存在一個平均突變,假設(shè):
備選假設(shè):
μ1≠μ2,當(dāng)Ha為真時,t=k 時的點被稱作間斷點。 Δ=被稱作平均突變的大小,最可能的間斷點服從以下分布:
其中,P(k)是建立的經(jīng)驗性的懲罰因子,其建立方法可參考[14]。
其中μ^0和β0是在μ1=μ2=μ 時的估計值。通過回歸檢驗算法來檢驗出多個間斷點,首先在t={Nmin,Nmin+1,…,(N-Nmin)}時段找出最可能的間斷點C0。類似的,分別在時段t={Nmin,…,(C-Nmin)}和t={(C0+1+Nmin),…,(N-Nmin)}找出次可能的間斷點C01和C03,接下來在C01和C03之間的時段找出最可能的間斷點C02,然后在C01、C02和C03之間尋找最可能的間斷點,即依次分段找出序列中各段最可能的間斷點,計算所有間斷點的統(tǒng)計量PFmax,估計自相關(guān)和P 的大小,找出最大的PFmax值對應(yīng)的間斷點,如果是顯著的,該突變就被認(rèn)為是找到的第一個間斷點,該間斷點被記為NC=1。尋找該間斷點位置之后每段最可能的間斷點,估計其顯著性,找出下一個可能的間斷點,重復(fù)該過程,將間斷點按照顯著性由大到小排列,形成間斷點列表。
PMF 方法可以使用元數(shù)據(jù):當(dāng)利用回歸檢驗檢測出多個間斷點后,輸入元數(shù)據(jù),再次判斷輸入的時間點是否存在間斷點(即依次分段找出序列中各段最可能的間斷點),計算所有間斷點的統(tǒng)計量,確定第一個間斷點;尋找該間斷點位置之后每段最可能的其他間斷點,分別估計其顯著性,找出下一個可能的間斷點,重復(fù)該過程,分步找出所有的間斷點。然后再將間斷點按照顯著性由大到小排列,形成新的有元數(shù)據(jù)支持的間斷點列表。判斷最小的間斷點是否顯著,當(dāng)不顯著時剔除該間斷點,再次評估剩余間斷點的顯著性,最終保留的統(tǒng)計顯著的間斷點即為序列檢驗得到的間斷點。該方法能夠檢驗的最短序列長度是20[14]。
RHtests 方法已被人們廣泛地應(yīng)用于氣溫、降水、風(fēng)速及相對濕度等氣候序列的均一性分析中。因喀什站高空各層次日溫度資料有大量缺測值,且待檢序列(1958—2016 年)與對應(yīng)的ERA-Interim 參考序列(1979—2013 年)時間不一致,在對喀什站高空日、年溫度資料進(jìn)行均一性檢驗和訂正的過程中發(fā)現(xiàn),在利用RHtestsV5 軟件包進(jìn)行訂正時,大量缺省值的代入造成了訂正序列與實際觀測序列的偏差量過大,導(dǎo)致訂正序列不可信。故本文利用RHtestsV5 軟件包提供的PMT 方法,采用對應(yīng)的ERA-Interim 作為參考序列,對喀什站高空日溫度資料進(jìn)行均一性檢驗和訂正;利用RHtestsV5 軟件包提供的PMF 方法,對喀什站高空年溫度資料進(jìn)行均一性檢驗和訂正。
喀什站常規(guī)探空和測風(fēng)的同步觀測任務(wù)起始于1957 年6 月1 日,以有整年觀測數(shù)據(jù)(1958 年1 月)以來的平均溫度序列按照850、700、600、500、400、300、200 hPa 和100 hPa 共8 個規(guī)定等壓面的月值(共708 個月)來統(tǒng)計缺測率及數(shù)據(jù)缺測年月分布狀況,統(tǒng)計結(jié)果見表1 和表2,由表可見,喀什站8 個規(guī)定層1967 年9—10 月均缺測。
根據(jù)表1 和表2 給出的缺測統(tǒng)計結(jié)果,08 時8個層次的月平均溫度除100 hPa 的缺測率>3%,完整性較差外,其他層次完整性較好;20 時100 hPa的月平均氣溫完整性較差,200 hPa 的完整性一般,其他6 個層次的完整性較好。綜合分析認(rèn)為喀什站850、700、600、500、400 hPa 和300 hPa 共6 個層次的高空溫度觀測資料序列完整性較好。
表1 喀什站高空溫度序列月值缺測情況(08 時)
表2 喀什站高空溫度序列月值缺測情況(20 時)
利用RHtestsV5 軟件包提供的PMT 方法,采用對應(yīng)的ERA-Interim 資料作為參考序列,分別對喀什站08 時、20 時對應(yīng)的850、700、500、400、300、200 hPa 和100 hPa 共14 條(因搜集的ERA-Interim資料中無對應(yīng)的600 hPa 序列,故此處沒有對喀什站600 hPa 的日氣溫序列進(jìn)行斷點檢測和分析。)溫度序列作置信度為0.95 的PMT 檢驗。保留其中的兩類斷點,一類是在年或月序列中檢驗出的顯著間斷點,并且有元數(shù)據(jù)的支持,如果該斷點出現(xiàn)的時間與臺站元數(shù)據(jù)記錄信息相差一年以內(nèi),根據(jù)元數(shù)據(jù)記錄時間替換該斷點位置;另一類是沒有元數(shù)據(jù)支持,但是在日和月(或年)序列中被同時檢驗出的相同年份的顯著間斷點,斷點位置依據(jù)日平均溫度序列被檢驗出的斷點時間。
檢驗得到的斷點信息顯示,喀什站14 個溫度序列被檢驗出的斷點共有221 個(圖2),其中有元數(shù)據(jù)支持的斷點占89.14%。從檢測結(jié)果和元數(shù)據(jù)分析情況看,喀什站始建于1954 年7 月,而20 世紀(jì)50 年代以來中國高空觀測儀器和規(guī)范變化較大。1957 年4 月高空觀測時間從北京時間11 時和23 時變更為08 時和20 時。1959 年前后、1966 年以及1999 年7 月1 日后又經(jīng)歷了輻射訂正方法的變化,1962 年5 月31 日、1968 年10 月1 日、1970年7 月1 日和2005 年喀什站又陸續(xù)更換探空儀器。綜合分析認(rèn)為,喀什站溫度資料非均一性主要來源于臺站位置變動、觀測儀器換型、觀測規(guī)范變更和輻射訂正。
圖2 喀什站溫度序列檢測出的斷點數(shù)量
通過PMT 方法和元數(shù)據(jù)信息相結(jié)合的方法來綜合判斷資料的斷點,所確定的斷點均通過了0.05的顯著性檢驗。其中,有一部分比例的斷點雖然有元數(shù)據(jù)支持,但調(diào)查數(shù)據(jù)時發(fā)現(xiàn)由于各種原因,在斷點處、前或后的數(shù)據(jù)為缺測。如,14 條數(shù)據(jù)序列均在1967 年9 月12 日—1967 年11 月10 日缺測,而檢測出得到的1967 年9 月12 日和1967 年11 月10日兩個斷點及其間缺測的數(shù)據(jù)不進(jìn)行訂正。另外,在對日原始數(shù)據(jù)序列進(jìn)行訂正時,僅對有參考序列對應(yīng)的時段進(jìn)行訂正。通過判斷,20 時700 hPa、08時500 hPa 和20 時100 hPa 的非均一性最為明顯,分別檢測出32 個、26 個和25 個斷點(表3),經(jīng)查詢臺站元數(shù)據(jù)信息,這3 條序列斷點的元數(shù)據(jù)支持率分別為84.38%、57.69%和88%,綜合考慮以上各種情況,這3 條序列最終確定需要訂正的斷點分別為2 個、6 個和2 個,占所檢測出的斷點比例僅分別為6.25%、23.08%和8%,做到了絕不無根據(jù)地訂正資料。
表3 2 個觀測時次7 個標(biāo)準(zhǔn)等壓面上斷點數(shù)
通過檢驗出的各規(guī)定層氣溫序列斷點的位置和訂正幅度來看,14 條序列均檢驗出了斷點,除了08時100 hPa 以外,其他13 條序列均檢驗出一個相同位置(1997-12-09)的斷點,經(jīng)查詢臺站元數(shù)據(jù)文件發(fā)現(xiàn),在上世紀(jì)90 年代前后,喀什站經(jīng)歷了輻射訂正方法的變化。由此可以認(rèn)為,輻射訂正方法的變化導(dǎo)致13 條序列產(chǎn)生了間斷點。另外,08 時850~400、100 hPa 和20 時400 hPa 共7 條序列因臺站位置的變動(1980-05-01 和1993-07-01)而產(chǎn)生了間斷點。在利用RHtestsV5 軟件包進(jìn)行斷點檢測和訂正時,序列中的缺測值需用缺省值“-999.99”來代替。在對各序列進(jìn)行訂正時發(fā)現(xiàn),300 hPa 以上各序列的訂正幅度非常大,其中20 時100 hPa 的日氣溫序列訂正量高達(dá)-184 ℃,這么大的訂正量是不可信的。通過查找原因,發(fā)現(xiàn)300~100 hPa 這6 條日值序列的缺測率較高(表4),而08 時和20 時各序列訂正量的大小與對應(yīng)的序列缺測率之間呈顯著的正相關(guān)關(guān)系,相關(guān)系數(shù)分別高達(dá)0.97 和0.99,由此不難理解,大量缺省值的代入,會導(dǎo)致訂正序列與觀測序列的偏差過大。為了保證訂正序列的可信度,本研究僅對缺測率小于1%的溫度日值序列進(jìn)行訂正,即僅對08 時和20 時850~400 hPa 共8 條序列進(jìn)行了訂正,具體訂正情況見表5。
利用RHtestsV5 軟件包中的最大懲罰F 檢驗(PMFT)方法對喀什站1958 年以來觀測的850~100 hPa 共14 條平均氣溫的年值序列進(jìn)行斷點檢驗。該方法經(jīng)驗性的考慮了時間序列的滯后一階自相關(guān),并嵌入多元線性回歸算法,能夠用于檢驗、訂正包含一階自回歸誤差數(shù)據(jù)序列的多個間斷點(平均突變),而且不需要使用參考序列。
表4 各規(guī)定等壓面日氣溫序列缺測率
表5 日溫度序列檢驗出的斷點位置和訂正量
針對PMFT 方法檢測出的年序列中的顯著間斷點,基于喀什站詳實的高空歷史沿革記錄,挑出儀器變更、臺站位置變動、觀測規(guī)范變更、觀測時間變化以及輻射訂正方法變化等有記錄的時間點(表6),來判斷在該時間是否為間斷點。如果該斷點出現(xiàn)的時間與臺站元數(shù)據(jù)記錄信息相差一年以內(nèi),根據(jù)元數(shù)據(jù)記錄時間替換該斷點位置。若檢測出的間斷點不顯著,也沒有元數(shù)據(jù)支持,則舍棄該間斷點,也不進(jìn)行訂正,做到絕不無根據(jù)地訂正資料。
表6 喀什站存在變化的時間點及原因
年平均氣溫資料序列均一性檢驗結(jié)果如圖3~8所示。
圖3 喀什站08 時100 hPa 年平均氣溫檢測結(jié)果
(不連續(xù)點:1997 年,輻射訂正方法變化)
圖4 喀什站08 時300 hPa 年平均氣溫檢測結(jié)果
圖5 喀什站08 時500 hPa 年平均氣溫檢測結(jié)果
圖6 喀什站08 時700 hPa 年平均氣溫檢測結(jié)果
圖7 喀什站08 時850 hPa 年平均氣溫檢測結(jié)果
(不連續(xù)點:1967 年,兩個月無觀測;1977 年,海拔高度變化、觀測手冊變化;1997 年,輻射訂正方法變化;2005 年,雷達(dá)和探空儀換型、L 波段系統(tǒng)運行、數(shù)據(jù)文件換型、觀測規(guī)范變化)
圖8 喀什站20 時100 hPa 年平均氣溫檢測結(jié)果
喀什站14 條溫度序列中有6 條序列即08 時100a、300、500、700、850 hPa 和20 時100 hPa 存在不連續(xù)點,其中,08 時850 hPa 的不連續(xù)點高達(dá)4個,20 時除100 hPa 外,其他序列均無不連續(xù)點。結(jié)合元數(shù)據(jù)分析,造成喀什站觀測資料不連續(xù)的原因主要有:1967 年有近2 個月停止觀測,1974 年觀測規(guī)范變化,1977 年海拔高度變化和觀測手冊變化,1997 年輻射訂正方法的變化,2005 年雷達(dá)和探空儀換型、L 波段系統(tǒng)運行、數(shù)據(jù)文件換型、觀測規(guī)范變化。其中,08 時850 hPa 年平均溫度在2005 年出現(xiàn)的斷點,很可能是由于L 波段與59-701 系統(tǒng)的感應(yīng)元件和測溫原理不同而存在的器差所致[15]。
經(jīng)過均一性分析和趨勢分析發(fā)現(xiàn),喀什站14 條高空溫度序列中,08 時850、700、500、300、100 hPa和20 時100 hPa 共6 條序列不連續(xù),需要訂正,其他8 條序列連續(xù),則不需訂正。通過PMFT 方法和元數(shù)據(jù)信息相結(jié)合的方法確定斷點后,主要采用了RHtestsV5 軟件包中的平均值訂正法對最終被采納的斷點進(jìn)行了訂正,即以序列的最新資料為標(biāo)準(zhǔn),把斷點前后的平均值差作為偏差,對斷點前序列進(jìn)行訂正,訂正結(jié)果見表7。最終采納的斷點數(shù)占PMFT方法檢測出斷點數(shù)的百分比為76.9%,也就是說有23.1%的斷點沒有被采納。
通過訂正,6 條氣溫序列訂正后與訂正前均呈現(xiàn)出相反的變化趨勢,其中,08 時100、300 hPa 和500 hPa 均為副訂正,而且訂正前氣溫均為下降趨勢,訂正后的氣溫均為上升趨勢;08 時700 hPa 和20 時100 hPa 則為正訂正,訂正后的氣溫均為下降趨勢??梢院苊黠@地發(fā)現(xiàn)訂正幅度隨著規(guī)定層高度的變化呈現(xiàn)出兩頭大中間小的現(xiàn)象,其中08 時和20 時100 hPa 的訂正幅度分別為-1.13 ℃和-0.80 ℃,而08 時850 hPa 的訂正幅度則高達(dá)2.56℃。08 時需要訂正的序列要多于20 時,08 時7 條序列中有6 條序列做了訂正,訂正率達(dá)到了85.71%,而20 時僅有100 hPa 做了訂正。
表7 年溫度序列檢驗出的斷點位置和訂正量
(1)喀什站08 時100 hPa 的高空月平均氣溫缺測率大于3%,完整性較差,20 時100 hPa 的完整性最差,缺測率高達(dá)11.44%;其他層次的完整性較好。
(2)喀什站14 條日溫度觀測序列共檢驗出221個斷點,14 條年平均氣溫序列共檢驗出10 個斷點。結(jié)合元數(shù)據(jù)分析認(rèn)為,喀什站高空溫度資料的非均一性主要來源于臺站位置變動、觀測儀器換型、觀測規(guī)范變更和輻射誤差訂正方法的變化。
(3)利用RHtestsV5 軟件包對喀什站日平均氣溫序列進(jìn)行斷點檢測和訂正時發(fā)現(xiàn),因缺測值需用缺省值“-999.99”來代替,大量缺省值的代入,會導(dǎo)致訂正序列與實際觀測序列的偏差過大。為了保證訂正序列的可信度,僅對缺測率小于1%的溫度日值序列進(jìn)行了訂正。
(4)喀什站高空年平均氣溫序列的訂正幅度隨著高度的變化呈現(xiàn)出兩頭大中間小的現(xiàn)象,08 時的訂正率明顯多于20 時,其中,08 時的訂正率達(dá)到了85.71%,而20 時僅有100 hPa 一條序列做了訂正。