李巖松,孔龍時,魏艷強,齊濤,特韋塔
(1.北京無線電測量研究所,北京 100864;2.中國氣象局氣象探測中心,北京 100081;3.西仰光科技大學,仰光11082)
2019 年年初至2020 年8 月,應北京人影辦要求,北京無線電測量研究所生產了X 波段相控陣天氣雷達,并在張家口進行氣象保障試驗。設備架設地點位于張家口市東花園鎮大隊院中,海拔50 m,與北京市區雷達高度接近。
這次試驗持續時間為20 個月,覆蓋了北京整個冬天、春天、夏天,設備運行穩定,獲取到了北京當年所有天氣過程。
X 波段相控陣雙偏振天氣雷達是北京無線電測量研究所自主研制的相控陣天氣雷達[1-4],雷達采用全固態發射接收模塊[5],工作穩定可靠。天線采用雙偏振體制,可以有效獲取多種雙偏振參數[6-7],擴充雷達二次產品類型。雷達在站點長期運行,經受了多次突然斷電、低溫、降雪、高溫、降水等過程的考驗,全程運行穩定,獲得了人影辦用戶的肯定。
這次試驗中,在分析降水過程時(這次采用2020年7 月2 日的數據),畫出了差分傳播相移[8]ΦDP分布散點圖,如圖1 所示,可以看出,ΦDP分布區間較廣,數值跨度較大,數據連續性較低,整體數據質量不高。抽取畫出其中一個單徑向ΦDP曲線(如圖2 所示)后可以看出,在近距離庫,ΦDP數值明顯屬于異常值,推測為信號處理將近距離雜波統一放入計算中,導致ΦDP數據質量降低。

圖1 ΦDP 散點圖

圖2 ΦDP 單徑向廓線圖
通過分析數據,可以找到當前信號處理輸出ΦDP所面臨的幾個問題:
1)在低仰角的雷達近距離庫,ΦDP數據沒有進行質量控制,有部分明顯為雜波的數據沒有得到很好的剔除。
2)在徑向數據中,存在少量明顯有較大偏差的數據沒有剔除,導致局部ΦDP曲線抖動較大。
針對以上問題,該文提出了數據質量控制算法,重新計算雷達ΦDP數據,主要方法為:
1)根據雷達數據間的關聯關系,加入針對ΦDP的數據質量控制,剔除明顯為地雜波的ΦDP數據,保留可反映降水的ΦDP數據。
2)根據ΦDP自身的物理特性,提出ΦDP的平滑算法,剔除局部變化較大的ΦDP數值,使整體ΦDP變化更加平滑。
雙極化參量ΦDP定義為垂直通道相位與水平通道相位的差值,主要是因為水平極化電磁波與垂直極化電磁波在空間行進,在經過小的非球形的液滴時,垂直極化的相位變化與水平極化相位變化速度不同,而ΦDP主要刻畫兩個通道在同一位置相位的差值。根據實際降水過程,當天氣雷達發射的雙極化波經過標準圓形液滴時,ΦDP值不會有任何變化,而當液滴不是標準圓形液滴時(較多天氣過程中液滴為橢球形),ΦDP會出現逐漸增加的情況。ΦDP是一個累積量,當經過降水過程區域時,ΦDP逐漸增加,而增加的速度與降水區域中液滴的橢球形態相關,一般情況下,強度越強的降水,ΦDP的增加速度也會提升。
根據ΦDP的特點,可以看出ΦDP本身屬于一個較平滑的逐漸遞增量,而當遇到地雜波等非氣象目標時,ΦDP則會出現不規則的抖動且偏離正常值較大的情況。
根據定義,地雜波處[9-11]的ΦDP沒有參考價值,可以利用相關系數[12-13]以及譜寬[14]的關系對地雜波數據進行剔除。
首先設置門限:
RhoHV>0.9;W=0。
根據相關系數的特性,在降水區域,雷達獲取到的相關系數數值一般均在0.95 以上,而雜波由于沒有相關性,相關系數一般較低,因此設置相關系數門限為0.9。
根據譜寬特性,一般降水區域,由于內部降水顆粒受局部上升氣流等因素的影響,速度整體呈現一定程度的不一致性,導致譜寬一般不為0,而地雜波由于無任何內部運動,譜寬一般為0。
根據以上兩個雷達參量,可以給出聯合的地雜波判決門限,即當譜寬W=0 且相關系數RhoHV<0.9時,判決同一位置所測到的數據為地雜波,對對應位置的ΦDP進行剔除操作。
根據2.1 節所述,在穿過降水或其他氣象過程時,ΦDP是一個逐漸增加的量,而實際中,由于收到其他非氣象目標、雷達系統噪聲或偶然誤差的影響,ΦDP的數值可能會出現小幅度的抖動。
針對ΦDP徑向數據出現的抖動情況,提出了徑向方向上的平滑濾波算法[15],利用濾波算法平滑ΦDP徑向方向上的變化,同時剔除可能出現的幅度變化較大的ΦDP數值。
在平滑時,為了不影響ΦDP本身的物理性質,參考2.1 節所述,ΦDP在經過強度數據較高的區域時,上升情況明顯,而在經過強度數據一般的區域時,上升情況放緩,因此設定以強度數據為依據的動態平滑窗。
針對強度數據從弱到中等的ΦDP數據回波(≤40 dBz),選取較長的滑動窗(設定滑動窗距離為6 km),而對于強度數據較強區域的ΦDP數據回波,設定較短的滑動窗(設定滑動窗距離為3 km)。
具體平滑流程如下:
1)選取對應長度的ΦDP數據,并計算數據的平均值。
2)當窗內的某一數據偏離平均值大于10°時,認為該數據有明顯偏差,屬于抖動異常值,需剔除。
3)向后移動數據窗,再次進行步驟2)。
該文選用的濾波算法為中值濾波算法,算法的主要原理:利用一個奇數的滑動窗,將滑動窗的中心點的數值用滑動窗中數值的中值代替,達到平滑濾波的目的,該濾波算法是一種非線性的濾波技術,可以很好地匹配ΦDP參數的物理特點。
利用2.3 節所述依據強度數據(Z)的滑動窗動態選取方法,具體濾波算法步驟如下:
1)在徑向數據方向上,選取對應長度的ΦDP數據,并對窗內的數據進行排序,形成新數組I。
2)找到滑動窗數據中中間位置的數值,將該數值替換為新數組I中間位置的數值。
3)向后移動數據窗,再次進行步驟2),直至到達徑向數據末尾。
根據2.2 節提出的地雜波剔除算法,對2020 年7 月2 日11:25 分的體掃第一層數據進行仿真。
從圖3 中第一層體掃強度可以看出,雷達300°~360°方向上有地雜波,而240°~300°方向上有降水回波,而ΦDP數值均顯示所有區域,容易造成數據質量較差。

圖3 第一層體掃強度
利用2.2 節提到的相關系數及譜寬門限,將地雜波數據進行剔除,得到的圖像如圖4 所示。

圖4 剔除地雜波及其他雜波后的ΦDP 圖像
從圖4 可以看到,雷達近區的雜波有了明顯改善。
根據2.3節提出的數據平滑算法,針對第2節同一徑向數據,利用Matlab同樣對2020年7月2日11:25分體掃第一層數據進行仿真,得到的結果如圖5 所示。

圖5 270°方向徑向ΦDP 濾波后廓線
對比圖2、圖5 可以看出,經過平滑濾波后,徑向上ΦDP有波動值的情況得到了明顯改善,大范圍抖動的ΦDP被濾除,ΦDP數值基本處于同一個區間,數據質量可信,通過徑向數據濾除統計,86%的地雜波數據被濾除。
根據2.4 節給出的ΦDP平滑濾波滑動窗濾波算法,對平滑后的數據進行中值濾波。這次選取2020年7月2日11:25分體掃第六層數據進行仿真,如圖6-10 所示。

圖6 第六層強度數據

圖7 第六層濾波前ΦDP 數據
根據圖6 可以看出,第六層(高仰角)數據中,地雜波影響有了明顯降低,數據主要反映了降水過程,大部分ΦDP數據均應處于合理范圍內,對比圖6、圖8 可以看出,在降水區域(強度>20 dBz),經過濾波算法后,數據很好地保留了ΦDP數據,經過計算統計,無效數據濾除率高于82%。對比圖9、圖10 可以看出,經過濾波后,數據抖動情況有了很好的改善,數據質量有了很大提高。

圖8 第六層濾波后ΦDP 數據

圖9 第六層全部ΦDP 數據濾波前散點圖

圖10 第六層全部ΦDP 數據濾波后散點圖
文中討論了X 波段相控陣天氣雷達在降水過程中ΦDP數據質量改善的問題,通過研究ΦDP數據的物理意義,利用相關系數以及譜寬等物理量,剔除了部分地雜波引起的ΦDP數值偏差問題。同時根據ΦDP自身的物理特性,提出數據的平滑算法以及中值濾波算法,更加有效地提高了ΦDP的數據質量,同時利用仿真驗證了濾波算法的效果,地雜波濾除率大于或等于85%,而無效數據濾除率大于或等于80%,ΦDP整體數據質量有明顯提升。
該文濾波采用的平滑濾波算法為中值濾波算法,后期應引入更好、更高效的濾波[16]和擬合算法,提升數據平滑效率。同時根據數據分析出現的穿過雨區的ΦDP隨距離逐漸減小的問題,需要檢查雷達設備信號處理算法符號,保證后期ΦDP數據的正確性。