高 霞, 周 曄, 周淑玥, 杜 捷, 王海江
(1.中國航空工業集團公司雷華電子技術研究所,江蘇 無錫 214000;2.成都信息工程大學電子工程學院,四川 成都 610225)
雷達電磁波在傳輸過程中,能量會被沿途接觸到的各類粒子吸收轉變成機械能、內能、熱能等[1],電磁波的衰減與波長成反比。X波段雙偏振天氣雷達由于其波長較短,在探測降雨時能量會被降水粒子吸收等,導致能量衰減,在遠距離處的弱雷達回波返回時很難被雷達接收機接收,導致遠距離處信息丟失,嚴重影響雷達的探測精度以及雷達探測資料的應用。
X波段雙偏振天氣雷達存儲有差分傳播相移ΦDP,該參數為距離累積量不受衰減影響[2],因此被廣泛用于反射率因子的衰減訂正。目前較為廣泛的衰減訂正方法是通過計算衰減率來完成反射率因子的訂正,但衰減率大都是給定一個固定的區間范圍[3],存在諸多弊端,如不能靈活適應不同降水粒子的衰減。因此可以根據差分傳播相移ΦDP不受衰減的影響來設計自適應算法完成衰減訂正工作[4-5]。
首先對單部雷達進行反射率因子衰減訂正處理,其次將衰減訂正處理后的3部X波段雙偏振天氣雷達反射率因子插值到統一的CAPPI網格上,再對組網反射率因子CAPPI進行缺空補值和斑點濾波處理,最后設計區域檢測算法,對組網反射率因子進行訂正。
3部X波雙偏振天氣雷達隸屬于北京市氣象局各區分局的雷達站點,分別是北京房山站雷達、北京昌平站雷達、北京順義站雷達,另外用于驗證的1部SA波段雷達為北京大興站雷達。各站點雷達的詳細參數如表1所示。

表1 雷達站點及詳細參數
研究數據來源于2020年夏季(8-9月)雷達探測數據(表1),由于夏季降水豐沛,雷達聯合組網后效果更佳。
因為衰減是雙程的累積量,故衰減訂正為

對于衰減率的計算,Bringi等[5]通過散射的模擬實驗發現反射率因子的衰減率AH與差分傳播相移率KDP基本上呈線性關系,并在實驗的基礎上提出了衰減訂正公式:

其中,α為衰減系數,Testud指出對于X波段雷達,衰減系數α在0.139~0.335[3],本文取其均值0.237。
為更好地適應不同降雨類型,Bringi等[5]對上節算法進行改進,命名為自適應約束算法。自適應約束算法在計算中使用了差分傳播相移ΦDP和水平反射率因子ZH,并根據誤差大小不斷調整衰減系數α,將其調整為最佳值。自適應算法流程多,迭代次數大,算法涉及大量區間內積分運算,對差分傳播相移有較高的質量要求,因此需要對其進行質控,該算法流程如圖1所示。

圖1 自適應算法流程
自適應算法流程可以概括為:對每個體掃數據的所有仰角層逐層訂正,每層進行逐徑向訂正,將每個待訂正的徑向分為若干區間,分別計算這些區間內對應的衰減系數α,再根據徑向內所有區間的衰減系數對該徑向進行訂正,然后完成整層數據的訂正,最后再完成整個體掃數據的訂正。算法具體步驟如下:
(1)給定衰減系數α的初始值為0.01。
(2)利用ZH、降水量I和差分傳播相移ΦDP計算衰減率AH。

其中,I為降水量,b為衰減指數常數,取值0.76~0.84,本文取均值0.8,α為區間待求衰減系數,ZH為訂正前水平反射率因子(雷達存儲值)。
(3)根據KDP和AH的關系、KDP和ΦDP的關系估算

(5)衰減系數α以0.05為步長,遞增,重復步驟(4),當兩者的差值最小時,即得到最佳衰減系數。
(6)應用最佳衰減系數α進行訂正。
因S波段天氣雷達波長較長,電磁波受雨區衰減影響較小,因此驗證數據集為北京大興站同時刻(兩雷達探測誤差在1分鐘以內)同區域內的S波段單偏振天氣雷達數據。通過兩種方法對不同地區雷達反射率因子數據進行訂正發現,固定衰減系數的AH-KDP法在近雷達處能夠對反射率因子有較弱訂正,在遠離雷達處效果不佳。AH-KDP法的反射率因子訂正前后距離廓線對比如圖2所示。反射率因子訂正前后對比PPI圖如圖3所示。而自適應算法在近、遠離雷達處效果都較明顯,能較好地恢復真實降雨區數據。自適應算法的反射率因子訂正前后距離廓線對比如圖4所示。反射率因子訂正前后對比PPI圖如圖5所示。

圖2 AH-KDP法訂正值與探測值及驗證值對比圖

圖3 AH-KDP法衰減訂正效果對比圖

圖4 自適應算法訂正值與探測值及驗證值對比圖

圖5 自適應算法衰減訂正效果對比圖
為更直觀反映兩種算法的優劣,對同種X波段雷達探測數據分別用兩種算法進行訂正,并以同時刻(誤差在1分鐘左右)S波段雷達探測數據作驗證,結果表明自適應算法明顯優于AH-KDP法,對比圖見圖6。

圖6 AH-KDP法與自適應算法效果對比圖
根據自適應算法和AH-KDP法的效果對比,本文最終選擇自適應算法對組網前的3部X波段雙偏振天氣雷達的反射率因子數據進行第一步衰減訂正處理。
常用的線性插值算法包括:垂直線性插值、垂直水平線性插值和八點線性插值算法。八點線性內插(EPI)選取網格點周圍8個相關聯的距離庫作為參考點,將參考維度擴展至仰角、方位角、徑向上,能夠提升內插的可靠性。該算法解釋為選取該網格點對應在極坐標系下最近的兩個仰角層和方位角(或徑向)上六面體的8個頂點作為參考點[6]。
如圖7所示,p1、p2、p3、p4為位于該網格點p上方仰角層內最近的4個距離庫參考點,p5、p6、p7、p8為位于該網格點p下方仰角層內最近的4個距離庫參考點;p1、p3、p5、p7 位于同一方位角,p2、p4、p6、p8 位于同一方位角;p1、p2、p5、p6擁有相同的徑向距離庫,p3、p4、p7、p8擁有相同的徑向距離庫。p1~p8對應坐標為 p1(e1,a1,r1)、p2(e1,a2,r1)、p3(e1,a1,r2)、p4(e1,a2,r2)、p5(e2,a1,r1)、p6(e2,a2,r1)、p7(e2,a1,r2)、p8(e2,a2,r2)。p1~p8對應值為V1~V8,待插值點p在空間笛卡兒坐標系下的坐標為p(x,y,z),需要將其轉換至極坐標系下的坐標p(e,a,r),轉換公式如式(7),最后用三個維度的線性插值對待插值點p的值Vp進行計算,計算公式如式(8)。

圖7 八點線性插值法示意圖

其中的we、wr、wa分別為各參考點在仰角、徑向距離和方位角上的權重系數,其計算公式如式(9)~(11)所示。

在得到多部雷達網格化處理后的雷達資料后,常見的拼圖算有以下幾種:最臨近法、最大值法、反距離權重函數法[7]。反距離權重函數法相較于最大值法和最臨近法,綜合考慮了組網所用的多部雷達反射率因子在同一網格點的貢獻,能夠明顯地減少其他兩種方法帶來的色標跳變現象,組網圖像更連續。因此,組網處理部分最終選擇反距離權重函數法,其示意圖如圖8所示。

圖8 反距離權重法示意圖
由圖8可知,3部雷達形成一個三角形區域,網格點p的值由該網格點到3部雷達R1、R2、R3距離決定,最后根據網格點p到3部雷達在三維x-y-z空間上的直線距離r1、r2、r3進行反距離加權,設 r1、r2、r3大小遞減。3部雷達的坐標為R1(x1,y1,z1)、R2(x2,y2,z2)、R3(x3,y3,z3),p點的坐標為p(x,y,z)。3部雷達在p點的值分別為V1、V2、V3,p點拼圖后的值應為Vp。拼圖公式如下:

對2020年8月12日10:00:00北京房山站、北京昌平站、北京順義站X波段雙偏振天氣雷達的水平反射率因子PPI以八點線性插值方法進行3 km高的CAPPI插值與組網處理,其插值結果、組網拼圖結果分別如圖9所示。
由圖9可知,EPI法在3個維度(仰角、方位角、距離)上選擇8個參考點,能夠對3部X波段雙偏振天氣雷達極坐標系下的反射率因子資料合適地插值到統一的空間笛卡爾坐標系下(CAPPI)。由圖10可知,反距離權重函數法基本能將3部雷達的插值結果較完整的組網到同一個網格上。至此,完成前期的3部雷達的反射率因子組網處理。

圖9 三地區插值處理CAPPI圖

圖10 組網拼圖處理CAPPI圖
預處理部分主要分為缺空補值和斑點濾波[8-9],關于反射率因子CAPPI的缺空補值處理,設計NA×NR窗口濾波函數[10]對所有網格點上的值進行濾波檢測。該窗口濾波函數可描述為:在每格網格點的第i個經度和第j個緯度處的有效值,其經度方向上的窗口尺度為NA,其緯度方向上的窗口尺度為NR,當窗口內的有效數據總數Ti,j在窗口內總庫數占比Pi,j小于閾值Thr1時,則將該庫視為孤立回波[11],并將該庫處的值設置為無效值NaN,如公式(15)。當窗口內的有效數據總數Ti,j在窗口內總庫數占比Pi,j大于閾值Thr2時,則將該庫視為缺測值庫,并將該庫處的值從NaN設置為一個有效值,該有效值F為該窗口內所有有效值的均值,如式(16)。

其中閾值Thr1設置為0.45,閾值Thr2設置為0.3,對北京房山、昌平、順義地區經過EPI法得到的反射率因子CAPPI進行缺空補值、斑點濾波處理前后的對比如圖11所示。

圖11 預處理前后CAPPI圖
由圖11預處理前后對比,可知預處理能夠對邊緣空缺值進行有效回填,并對組網區域內的斑點雜波進行相應的去除。
關于組網反射率因子CAPPI的訂正處理,由于在應用反射率因子進行CAPPI組網拼圖前,已經對3部雷達的反射率因子以自適應算法進行訂正處理,但這種處理不夠徹底。因為訂正是基于雷達電磁波在發射方向上不斷衰減能量的前提下,運用最佳衰減系數和差分傳播相移使探測值盡量靠近真實值。因此在反射率因子CAPPI組網之后,也要對網格點上的反射率因子值進行適當訂正。本文設計使用區域檢測算法對CAPPI每個網格點進行訂正。該算法整體思路為以S波段雷達為參考,對拼圖CAPPI和S波段雷達CAPPI的每個網格點建立維度為N的搜索框,并分別計算其均值Mean_S和Mean_X,當前者與后者的差值大于閾值Thr_C時,則判定為偏差嚴重,該偏差包括系統誤差和衰減,需要進行訂正。為防止訂正后的值大于真實值,以盡量靠近真實值的為原則,訂正后的值Vp為當前拼圖CAPPI值VX加上Mean_S與Mean_X的差值,其中N取5,閾值Thr_C取2 dBZ。建立的搜索框示意圖如圖12。

圖12 搜索框示意圖
訂正公式如式(17),組網反射率因子CAPPI訂正前后對比效果如圖13。


圖13 訂正處理前后對比圖
由圖13可知,訂正部分主要在39°N~40.3°N、115.5°E~118°E附近,將訂正后的組網反射率因子CAPPI與同時刻同地點的S波段雷達反射率因子CAPPI進行對比,如圖14。

圖14 訂正CAPPI與驗證CAPPI對比圖
由圖 13 可 知,115.5°E~118°E,39°N~40.3°N為主要訂正區域,區域檢測算法以待訂區域內同時刻下S波段天氣雷達探測資料為參考,對每個待訂正的網格點反射率因子建立搜索框,并計算兩個搜索框內的區域均值,當兩均值差異較大時,則判定為偏差過大,利用S波段雷達資料對其進行訂正處理。由圖14可知,主要訂正區域內的訂正結果基本與驗證用S波段雷達反射率因子CAPPI相吻合,該算法訂正效果表現良好。
在前人研究的基礎上,首先通過分析2種常規衰減訂正算法的優劣,選擇較優算法對單部雷達反射率因子訂正處理;其次運用經典八點線性插值算法和反距離權重函數法對訂正后的3部X波段雙偏振天氣雷達反射率因子進行網格化處理,最后設計了區域檢測算法對組網反射率因子CAPPI進行適當的訂正處理,實驗結果表明,本文設計的區域檢測算法在CAPPI網格上表現良好的效果,具有一定的參考價值,但該方法局限于組網區域內需要存在衰減較小的長波長天氣雷達用以驗證算法效果,極大限制了該方法在業務上的應用,未來需設計其他數值訂正方案來優化該算法。