夏 凡 張樂堅 張 林
1 山東省氣象科學研究所,濟南 250031 2 中國氣象局氣象探測中心,北京 100081
提 要:為了提高中國X波段雙線偏振雷達(簡稱雙偏振雷達)的數據質量,提升應用水平,基于美國多雷達多傳感器定量降水估測業務系統質量控制(簡稱dpQC)方案,在此基礎上增加了差分反射率閾值判別并對方案中的計算參數進行本地化調整,設計了針對X波段雙偏振雷達的質量控制(簡稱dpXQC)方案,原理主要包括:將雙偏振參量相關系數(CC)區分降水回波與非降水回波的閾值調整為0.9,小于0.9的回波判別為非降水回波;通過計算回波頂高與強回波區位置對冰雹與非均一性波束填充區域進行標識;通過統計探空資料0℃高度上下相關系數的特征,調整了融化層的識別參數;優化算法將雙偏振參量差分反射率絕對值大于5 dB判別為非降水回波;加入了條狀雜波識別、連續性檢驗與斑塊雜波識別算法。利用dpXQC方案對北京五部X波段多普勒雙偏振業務雷達進行試驗,檢驗評估其對非降水回波識別效果,獲得以下結果與結論:大部分非降水回波可以通過CC閾值判別被濾除;對于CC閾值判別無法濾除的條幅狀、孤立點、塊狀非降水回波,dpXQC方案利用條幅狀雜波識別、連續性檢驗與差分反射率閾值判別去濾除;dpXQC方案有效保護了冰雹、非均一性波束填充與融化層區域中CC<0.9的降水回波不被濾除;檢驗結果顯示,dpXQC方案對非降水回波正確識別率為91.8%,錯誤識別率為20.6%,均要優于dpQC方案,但dpXQC方案的錯誤識別率依然很高,主要有兩方面原因,一是由于X波段雷達反射率衰減嚴重,方案會對一些離雷達較遠的非均一性填充區域出現遺漏識別,二是融化層識別算法在所有方位角的位置與厚度是固定的,與實際觀測不符,這均會導致被遺漏標識區域中CC<0.9的降水回波被誤判為非降水回波??傮w來看,dpXQC方案在濾除非降水回波的同時,最大程度保護了降水回波不被濾除,為雷達資料定量化應用等相關業務提供保障。
我國目前正逐步將CINRAD單線偏振雷達(簡稱單偏振雷達)升級改造為雙線偏振雷達(簡稱雙偏振雷達),與單偏振雷達相比,雙偏振雷達可以接收水平和垂直方向的回波信號,除了反射率(ZH)、徑向速度(Vr)、譜寬(SPW),雙偏振雷達還可以探測接收零滯后相關系數(CC)、差分反射率(ZDR)、差分傳播相移(ФDP)等雙偏振參量。雙偏振雷達資料在強對流天氣預警(馮晉勤等,2018)、定量降水估測(Zhang et al, 2016;陳超等,2019)、水凝物粒子分類(馮亮等,2018)領域發揮了重要的作用,我國自主研發的GRAPES(Global and Regional Assimilation and Prediction System)數值預報模式(Chen et al,2015),也正嘗試將雙偏振雷達數據同化進模式初始場中,提高強對流天氣的預報效果。然而這些不可避免地會受到非降水回波的影響,高質量的數據是保證雷達自動化應用的關鍵,因此雷達基數據在應用前需要進行質量控制,識別并濾除非降水回波。
模糊邏輯法(Kessinger et al,2003;劉黎平等,2007;江源等,2009)與神經網絡法(Lakshmanan et al,2007)是較為成熟的雷達質量控制算法,在區分降水回波與非降水回波中取得了較好的效果。隨著雙偏振參量的引入,研究人員可以在原有算法的基礎上增加更多判別條件。杜牧云等(2013a)在模糊邏輯算法中增加了CC、ZDR與ФDP紋理隸屬成員函數,改善了算法對零速度區降水回波錯誤識別問題,同時提升了混合型回波的識別效果,但是仍遺漏了不少地物雜波;Lakshmanan et al(2013)在原有神經網絡算法的基礎上,訓練數據中增加了CC、ZDR與ФDP參量來區分降水與非降水回波,結果顯示識別技巧要高于原有的神經網絡法,對混合型回波的識別效果較好,然而該方法質量控制效果依賴于搭建網絡模型的訓練數據集。Tang et al(2014)基于ZH、CC與探空溫度數據,針對雙偏振雷達設計了一套質量控制(dual polarization quality control,dpQC)方案,其核心是利用CC識別大部分非降水回波與降水回波,然后利用回波頂高等參數處理CC無法區分的回波。統計結果顯示,dpQC方案的質量控制技巧高于神經網絡法,而且計算時間遠少于神經網絡方法。這一方法已經在美國多雷達多傳感器(MRMS)定量降水估測(QPE)系統中業務應用。在國內,張林和楊洪平(2018)借鑒了dpQC算法流程,分類識別上海南匯業務雙偏振雷達的降水回波與非降水回波。還有一些研究人員單純對一個雙偏振參量進行質量控制(杜牧云等,2013b;吳林林等,2015;肖艷姣等,2012;趙川鴻等,2019;馬建立等,2019),為后續反射率質量控制做好鋪墊。
目前已經業務化的雙偏振雷達以S與C波段為主,很多質量控制工作也是圍繞這兩種波段的雷達數據展開,但是兩者在局地強對流的精細觀測方面有所不足。與S、C波段雷達相比,X波段雷達觀測數據具有更高的時空分辨率,可以提供更細微水凝物信息,在災害天氣預警中發揮了更重要的補充作用。如北京自2015年搭建了5部X波段雙偏振雷達并投入到科研與業務中,在隨后的幾年應用中發現X波段雷達容易受到環境噪聲與信號衰減等因素影響,觀測數據摻雜了晴空回波、電磁干擾雜波與生物回波等非降水回波。目前國內X波段雙偏振雷達非降水回波識別的相關研究較少(郭春輝等,2019;王超等,2019),研究人員基于CC去識別非降水回波,這可能將非均一填充與融化層區域中的降水回波錯誤識別為非降水回波。為了彌補以往研究的不足,本文對Tang et al(2014)設計的dpQC方案進行優化改進,設計了針對我國X波段雙偏振雷達的質量控制方案來識別并濾除非降水回波,為業務化應用X波段雙偏振雷達數據提供更好的質量控制方案。
本文使用了北京昌平站、房山站、密云站、順義站與通州站5部X波段雙偏振天氣雷達在2019年5—7月觀測的基數據資料,雷達采用VCP21體掃模式,每個體掃有9個仰角,大約3 min完成一次體掃,距離庫長為75 m。用于計算融化層高度的0℃高度來自北京探空數據。
根據探空站距離昌平、房山、密云、順義與通州雷達站距離大致為42、33,70、39與31 km,因此,探空站0℃高度可以代表5部雷達。
由于北京X波段雷達在生成基數據時已經做了地物回波的濾除處理,生物回波個例較難收集,本文設計的質量控制方案識別的非降水回波主要是電磁干擾回波與晴空回波。圖1給出了方案的整體流程,包含了7個步驟,分別是:(1)CC閾值判別;(2)冰雹、非均一性波束填充(NBF)區域識別;(3)融化層識別;(4)ZDR閾值判別;(5)條狀雜波濾除;(6)水平連續性檢驗;(7)斑塊雜波濾除并在降水回波區域彌補小面積空洞。雷達基數據依次經過這7步自動化處理,最后保留的即為降水回波。

圖1 北京X波段雙偏振雷達質量控制方案流程圖
與原始dpQC方案(Tang et al,2014)相比,本文的質量控制方案移除了CC紋理閾值判別,增加了ZDR閾值判別,這是由于通過大量質量控制試驗發現:CC紋理閾值判別會濾除大量降水回波,ZDR閾值判別會去除面積較大的電磁干擾雜波;方案還基于統計結果對CC閾值判別與融化層區域識別的計算參數進行調整。本文將新的質量控制方案稱為dpXQC(dual polarization X-band quality control),下面詳細介紹各步驟原理。
1.2.1 CC閾值判別
CC表示回波的水平偏振分量與垂直偏振分量的相關程度。研究表明,單一類型水凝物降水回波,其對應的CC較大,通常接近于1.0,而非降水回波的CC較小(Kumjian,2013)?;贑IMISS數據庫中國地面逐小時資料的降水量與總云量、地面分鐘降水資料的降水量數據選取了非降水回波與降水回波的個例(表1),統計兩者CC的概率分布(圖2)。由圖可見,降水回波CC主要分布在0.8~1.0區間內,非降水回波分布范圍較廣,在0.2~1.0。

表1 2019年5—7月降水過程時間及體掃數量
CC閾值的選取直接影響方案對非降水回波與降水回波的識別效果,由圖2可知,閾值選取較低可減少對降水回波的錯誤識別,但會使少量非降水回波無法被濾除;反之則可以有效地識別非降水回波,但降水回波可能被錯誤識別。綜合考慮對非降水回波濾除效果與降水回波錯誤識別因素,本文以0.05為步長,依次計算CC區間[0.5,1]對應的降水回波概率分布閾值右側的累計概率值與非降水回波概率分布閾值左側累計概率值之和,從選出最大值,最終將0.9定為閾值。如果雷達回波對應CC<0.9,將其判定為非降水回波。
1.2.2 冰雹或者非均一性波束填充區域判別
從圖2可以看出,小部分降水回波對應的CC<0.9。Ryzhkov et al(2005)研究發現,冰雹區經常出現部分回波對應的CC較小,Ryzhkov(2007)還發現NBF區域也會出現較小的CC,這主要是由于隨著雷達探測波束的延長,單位距離庫內所探測的體積不斷增大,強降水回波中心占據單位探測體積的一部分或者不同相態的水凝物會出現在同一單位探測體積中,使得ФDP梯度變大導致CC降低。由于冰雹或者NBF的低CC區域通常出現在強降水天氣中,其對應的回波較強,回波頂較高,因此利用公式對其進行識別:

圖2 北京X波段雷達降水回波與非降水回波對應的CC的概率分布
式中:ETOP18 dBz與ETOP0 dBz分別表示被檢測的距離庫之上18 dBz與0 dBz反射率的高度,rstormcore表示為強回波區到雷達站的徑向距離,定義強回波區為某一徑向大于45 dBz回波累計距離大于1 km的區域。被標識為冰雹與NBF區域中CC<0.9的回波不被濾除。
1.2.3 融化層亮帶識別
融化層包含了融化的雪、霰與雹等不同相態的粒子,因此雷達探測到部分回波對應的CC較小。為了使融化層融區域中CC<0.9的降水回波不被濾除,借鑒了dpQC方案,即基于探空資料的0℃高度與CC來對其進行識別,步驟如下:
(1)首先將每個徑向的雷達數據分為三組,融化層之下(0℃高度其下1 km至其下2 km),融化層(0℃高度至其下1 km)與融化層之上(0℃高度至其上1 km),計算每組CC平均數(分別為CCb、CCin、CCa),公式如下:
(3)
式中:H0℃表示探空資料0℃高度;h(i)表示第i個距離庫的高度;n表示每個分組中距離庫的總數。
(2)當某個徑向上三組數據CC的平均值滿足下式時,這個徑向距離庫高度在[H0℃-1 km,H0℃)區間內被判定為融化層
CCin≥0.9∩{[CCin<(CCb-0.02)∩CCin<
(CCa-0.02)]∪CCin<(CCb-0.05)}
(4)
(3)在判定的融化層區域中,小于0.7的回波被標識為非降水回波并濾除,CC在[0.7,0.9)剩余的回波暫不處理,需要在接下來的步驟中繼續判別回波屬性。
本文統計了177個出現融化層的仰角,計算了融化層之下、融化層之中與融化層之上的CC的平均值(圖3),可以發現,融化層組的CC平均值在0.85~0.97波動,而且部分仰角融化層之上、之下CC平均值與融化層差較小,如果基于dpQC方案對融化層進行識別,會出現遺漏,基于統計結果,dpXQC修改了式(4)中的閾值,新的判別條件如下:

圖3 177個仰角融化層之上、融化層與融化層之下CC的平均值
CCin≥0.85∩{[CCin<(CCb-0.01)∩CCin<
(CCa-0.01)]∪CCin<(CCb-0.03)}
(5)
1.2.4ZDR閾值判別
這一步是在dpQC方案基礎上新增的,因為在實際觀測中,會出現面積較大并且CC較大的電磁干擾雜波,這些回波無法利用CC閾值判別去濾除,也無法利用水平連續性檢驗去濾除。通過大量觀測個例發現,ZDR可以較為有效識別這類非降水回波,ZDR反映了粒子偏離球形的情況,這些電磁干擾雜波的ZDR通常偏大或者偏小。本文統計了這類電磁干擾回波與降水回波二者的ZDR概率分布(圖4)發現,這類回波的ZDR分布在-7~8 dB,大概率分布在絕對值較大的ZDR區間,降水回波則分布在-4~5 dB,為了盡可能保留降水回波,本文將︱ZDR|>5 dB做為區分此類電磁干擾回波與降水回波的判別閾值。

圖4 北京X波段雷達降水回波與非降水回波對應ZDR的概率分布
1.2.5 條狀雜波濾除
條狀雜波通常由太陽射線或者電磁干擾造成,有的可以利用CC閾值判別進行濾除,但是也存在整條徑向回波的CC>0.9的情況。條狀回波的主要特征是徑向上連續性較好,并且條狀雜波通常出現在低層仰角,在垂直方向缺乏連續性。基于這些特征,本文質量控制算法中利用下式來識別條狀雜波:
(6)

1.2.6 連續性檢驗
實際觀測中常出現散點狀回波,這些回波可能是晴空回波,也可能是電磁干擾回波,散點狀回波的特點主要是具有孤立性,在方位與徑向上都不連續,這些回波有時CC>0.9或者對應CC為缺測。算法利用連續檢驗可以濾除這類非降水回波,其檢驗原理是,以被檢測的距離庫為中心,設定一個窗口并統計范圍內的ZH,如果多于一半是缺測值,或者有效回波ZH的平均值小于被檢測回波ZH的25%,那么該距離庫的回波被判定為非降水回波,本文將這個窗口大小設置為0.75 km×2.0°。
1.2.7 斑塊雜波濾除
經過前面幾步質量控制后,計算剩余的回波的面積,如果小于10 km2,那么這個區域被判別為斑塊雜波并濾除。這一步可以濾除利用CC、ZDR閾值判別與連續性檢驗無法濾除的非降水回波。
下面利用上節介紹dpXQC方案中各算法進行試驗,為了驗證效果,每一步算法單獨進行試驗,不與其他算法協同試驗。
圖5給出了2019年5月26日04:42密云雷達站回波分布,雷達站周圍存在大量散點狀電磁干擾回波與晴空回波(圖5a),對應的CC主要分布在0.1~0.5(圖5b),圖5c顯示了用CC閾值判別后的回波分布,可以發現雷達站周圍大部分非降水回波被濾除。但是仍然存在很多散點狀或小面積非降水回波保存下來,這些回波對應的CC為缺測或者大于0.9,后期可以利用連續性檢驗或者斑塊回波判別進行濾除。

圖5 2019年5月26日04:42密云雷達站0.5°仰角觀測的(a)質量控制前ZH,(b)CC,(c)CC閾值判別后ZH
圖6a給出了2019年5月17日22:29北京順義站雷達0.5°仰角質量控制前的ZH分布,在距離雷達100 km西南偏南方向有些回波的反射率達到45~60 dBz,實際觀測中這一區域出現冰雹,這些區域有的回波對應的CC較低(圖6b),只利用CC閾值判別,可以濾除雷達站周圍大部分非降水回波(圖6c上部橢圓),但是部分冰雹區與其附近區域(圖6c下部橢圓)的降水回波也會被濾除;對冰雹區域進行識別并標識后(圖6d),使得區域中的降水回波(圓圈區域)不被濾除,但是冰雹區附近區域降水回波沒有保留下來,這里主要由于不同相態的降水粒子出現在在同一距離庫中形成NBF區域,加入NBF區域識別(圖6e),有效保護了這一區域的降水回波。

圖6 2019年5月17日22:29順義雷達站0.5°仰角觀測的(a)質量控制前ZH,(b)CC,(c)CC閾值判別后ZH(上部橢圓內為經過CC閾值判別濾除的非降水回波,下部橢圓內為錯誤濾除的降水回波),(d)冰雹區識別后保留的ZH(圓圈內為經過冰雹識別后保留的降水回波),(e)冰雹與NBF區域識別保留的ZH
圖7a和7b分別給出了2019年7月5日20:03房山站雷達2.5°仰角質量控制前的ZH與CC分布,在距離雷達50~90 km、45°~180°方位角,大量回波的ZH分布在25~35 dBz,CC分布在0.6~0.9,結合當日的探空數據,這些區域出現在0℃層高度附近。單純利用CC閾值判別會濾除大量降水回波(圖7c橢圓區域內),對融化層區域進行識別后(圖7d),保留了其中大部分CC<0.9降水回波。

圖7 2019年7月5日20:03昌平雷達站2.5°仰角觀測的(a)質量控制前ZH,(b)CC,(c)CC閾值判定后ZH(橢圓內白色區域表示被錯誤濾除的降水回波),(d)融化層識別后保留的ZH
圖8a給出了2019年7月13日12:15密云站雷達站0.5°仰角觀測,在距離雷達站50~200 km、180°~240°方位角區域出現大量電磁干擾雜波,這些回波對應的CC大部分在0.9以上(圖8b),無法通過CC閾值判別來濾除,而其對應的ZDR>5 dB(圖8c),利用ZDR閾值判別后可以大部分電磁干擾回波被剔除(圖8d)。
圖9a給出了2019年5月25日15:21房山站雷達站0.5°仰角回波分布,在雷達站南側出現條狀雜波,其對應的CC>0.9(圖9b),無法利用CC閾值判別進行濾除,利用了條狀雜波識別算法,部分雜波可以被濾除(圖9c),剩余的回波主要由不連續的散點狀或者小面積回波構成,可以利用連續性檢驗或者斑塊回波識別進行濾除。
圖10a和10b分別給出了圖5c與圖8c經過連續性檢驗后回波分布,通過對比可見,利用CC閾值判別或者條狀雜波判別無法濾除的非降水回波,利用連續性檢驗可以將其濾除。

圖8 2019年7月13日12:15密云站雷達0.5°仰角觀測的(a)質量控制前ZH,(b)CC,(c)ZDR,(d)ZDR閾值判別后ZH

圖9 2019年7月16日15:21房山雷達站0.5°仰角觀測的(a)質量控制前ZH,(b)CC,(c)質量控制后ZH

圖10 經過連續性檢驗后2019年(a)5月26日04:42北京密云雷達站0.5°仰角觀測和(b)7月16日15:21房山雷達站0.5°仰角觀測的ZH分布
圖11給出了2019年5月26日11:54昌平站雷達0.5°仰角觀測經過斑塊狀雜波濾除前后的ZH分布,對比圖11a與11b可以發現,距離雷達站100~200 km的90°~120°方位角區域的斑塊狀雜波被濾除。

圖11 2019年5月26日11:54昌平站雷達0.5°仰角ZH分布(a)質量控制前,(b)質量控制后
為了評估質量控制方案對非降水回波的識別效果,借鑒了文浩等(2016)提出的評估方案,以一個仰角層的觀測數據為單位,先人工觀測進行判別,選出有、無非降水回波的仰角層數據,當質量控制方案識別非降水回波的面積超過人工判別非降水回波面積的90%,定為正確識別,當質量控制方案識別非降水回波區域在人工觀測判別為降水回波,且其面積占降水回波的10%,定為錯誤識別,從北京5部X波段雙偏振天氣雷達站選取了320個體掃、共811個仰角層進行測試。
統計指標基于二分類事件評估方法,主要包括正確識別率(Hr)與錯誤識別率(Far),公式如下:
(7)
(8)
式中:a為人工觀測與質量控制算法均判別出現非降水回波的仰角數;b為人工觀測判別沒有出現非降水回波而質量控制算法判別錯誤的仰角數;c為人工觀測判別出現非降水回波而質量控制方法判別遺漏的仰角數;d為人工觀測與質量控制算法均沒有判別出現非降水回波的仰角數。
在選取所有仰角層數據中,含有非降水回波的仰角層為123個,降水回波的仰角層為688個(有的仰角層也有非降水回波,但是比例非常小,忽略不計)。表2給出了dpQC方案與dpXQC方案質量控制效果統計,由表可見,非降水回波識別正確的仰角層數,dpXQC方案多于dpQC方案,dpXQC方案Hr為91.8%,dpQC方案的為85.4%,這主要是由于dpXQC方案增加了ZDR閾值判別,有效濾除了塊狀電磁干擾回波;降水回波被誤識別為非降水回波的仰角層數,dpXQC方案少于dpQC方案,dpXQC方案的Far為20.6%,dpQC方案的為28.8%,主要是由于dpXQC方案調整了融化層識別參數,使區域中CC<0.9的降水回波更多地被正確識別,并且老方案中包含了CC紋理閾值判別,這一項有時會將降水回波判別為非降水回波濾除,dpXQC方案刪除了這一判別項。

表2 dpQC方案與dpXQC方案的質量控制效果統計
然而,dpXQC方案對非降水回波錯誤識別率仍然很高,統計發現主要有兩點原因:一是由于X波段雷達反射率隨著距離的延長較S、C波段衰減更為嚴重,因此對于離雷達站較遠的NBF區域,利用算法的原有參數去識別會出現錯誤,因而這些NBF區域CC<0.9的降水回波會被識別為非降水回波,目前由于缺乏個例累積,還無法對原始計算參數做出修訂;二是融化層的識別算法是基于探空資料0℃層高度數據并統計其附近的CC特征來識別,探空數據一天只有三次,不能反映各時刻0℃層真實高度,同時算法中設定每個方位角的融化層位置相同且厚度固定為1 km,這與真實的融化層并不相符,因此識別區域外的真實融化層區域CC<0.9的降水回波被錯誤識別為非降水回波。
本文借鑒了MRMS-QPE系統dpQC方案原理,并根據北京X波段多普勒雙偏振業務雷達大量觀測個例對其進行優化改進,設計了dpXQC方案,用其進行質量控制試驗并且批量檢驗了dpXQC方案對非降水回波的識別效果,獲得以下結果與結論:
(1)根據降水回波與非降水回波CC概率分布特征,dpXQC方案將區分二者的閾值調整為0.9,CC<0.9的回波判別為非降水回波,試驗結果顯示這可以有效去除雷達站周圍的晴空回波與電磁干擾回波等非降水回波。
(2)冰雹區與NBF區域容易出現CC<0.9的降水回波,通過計算回波頂高與強回波區位置來對這些區域進行標識,標識區域內不進行CC閾值判別,可以有效保護這些區域CC<0.9的降水回波。
(3)融化層由于混合了多種相態水凝物粒子,部分降水回波的CC<0.9,dpXQC方案通過統計探空資料0℃層高度上下CC的特征來調整融化層的識別參數,以此保護融化層區域內CC<0.9的降水回波。
(4)dpXQC方案增加了ZDR閾值判別進行優化,通過統計降水回波與非降水回波ZDR概率分布特征,將|ZDR|>5 dB的回波判別為非降水回波,試驗分析表明ZDR閾值判別可以去除CC>0.9的塊狀電磁干擾回波。
(5)對于無法利用CC與ZDR閾值判別濾除的非降水回波,如條狀、散點與面積小于10 km2的非降水回波,算法通過條狀雜波識別、連續性檢驗與斑塊雜波識別分別來濾除。
(6)統計檢驗結果表明,dpXQC方案對非降水回波正確識別率為91.8%,錯誤識別率為20.6%,優于dpQC方案;其錯誤識別率依然較高,主要有兩方面原因,一是NBF區域的識別算法的參數沒有進行本地化調整,導致部分NBF區域沒有被標識,這些區域內CC<0.9降水回波被誤判為非降水回波;二是融化層識別算法在所有方位角的位置與厚度是固定的,與實際觀測不相符,這會造成被遺漏識別的融化層區域中CC<0.9的降水回波被錯誤識別為非降水回波。
總體來看,本文設計的dpXQC方案,可以濾除條狀、散點狀與塊狀的電磁干擾回波和晴空回波等非降水回波,并最大程度地保護了冰雹、NBF與融化層區域中的降水回波不被濾除。同時也要看到,在部分NBF區域與融化層區域,降水回波會被錯誤判別為非降水回波。對于離雷達站較遠的NBF區域識別,需要通過積累大量個例去修訂原有的計算參數。對于融化層的識別,還要研究或借鑒更新的識別算法,這是后續改進的方向。