高 霞, 周 曄, 周淑玥, 杜 捷, 王海江
(1.中國航空工業(yè)集團(tuán)公司雷華電子技術(shù)研究所,江蘇 無錫 214000;2.成都信息工程大學(xué)電子工程學(xué)院,四川 成都 610225)
雷達(dá)電磁波在傳輸過程中,能量會被沿途接觸到的各類粒子吸收轉(zhuǎn)變成機(jī)械能、內(nèi)能、熱能等[1],電磁波的衰減與波長成反比。X波段雙偏振天氣雷達(dá)由于其波長較短,在探測降雨時能量會被降水粒子吸收等,導(dǎo)致能量衰減,在遠(yuǎn)距離處的弱雷達(dá)回波返回時很難被雷達(dá)接收機(jī)接收,導(dǎo)致遠(yuǎn)距離處信息丟失,嚴(yán)重影響雷達(dá)的探測精度以及雷達(dá)探測資料的應(yīng)用。
X波段雙偏振天氣雷達(dá)存儲有差分傳播相移ΦDP,該參數(shù)為距離累積量不受衰減影響[2],因此被廣泛用于反射率因子的衰減訂正。目前較為廣泛的衰減訂正方法是通過計算衰減率來完成反射率因子的訂正,但衰減率大都是給定一個固定的區(qū)間范圍[3],存在諸多弊端,如不能靈活適應(yīng)不同降水粒子的衰減。因此可以根據(jù)差分傳播相移ΦDP不受衰減的影響來設(shè)計自適應(yīng)算法完成衰減訂正工作[4-5]。
首先對單部雷達(dá)進(jìn)行反射率因子衰減訂正處理,其次將衰減訂正處理后的3部X波段雙偏振天氣雷達(dá)反射率因子插值到統(tǒng)一的CAPPI網(wǎng)格上,再對組網(wǎng)反射率因子CAPPI進(jìn)行缺空補(bǔ)值和斑點濾波處理,最后設(shè)計區(qū)域檢測算法,對組網(wǎng)反射率因子進(jìn)行訂正。
3部X波雙偏振天氣雷達(dá)隸屬于北京市氣象局各區(qū)分局的雷達(dá)站點,分別是北京房山站雷達(dá)、北京昌平站雷達(dá)、北京順義站雷達(dá),另外用于驗證的1部SA波段雷達(dá)為北京大興站雷達(dá)。各站點雷達(dá)的詳細(xì)參數(shù)如表1所示。
表1 雷達(dá)站點及詳細(xì)參數(shù)
研究數(shù)據(jù)來源于2020年夏季(8-9月)雷達(dá)探測數(shù)據(jù)(表1),由于夏季降水豐沛,雷達(dá)聯(lián)合組網(wǎng)后效果更佳。
因為衰減是雙程的累積量,故衰減訂正為
對于衰減率的計算,Bringi等[5]通過散射的模擬實驗發(fā)現(xiàn)反射率因子的衰減率AH與差分傳播相移率KDP基本上呈線性關(guān)系,并在實驗的基礎(chǔ)上提出了衰減訂正公式:
其中,α為衰減系數(shù),Testud指出對于X波段雷達(dá),衰減系數(shù)α在0.139~0.335[3],本文取其均值0.237。
為更好地適應(yīng)不同降雨類型,Bringi等[5]對上節(jié)算法進(jìn)行改進(jìn),命名為自適應(yīng)約束算法。自適應(yīng)約束算法在計算中使用了差分傳播相移ΦDP和水平反射率因子ZH,并根據(jù)誤差大小不斷調(diào)整衰減系數(shù)α,將其調(diào)整為最佳值。自適應(yīng)算法流程多,迭代次數(shù)大,算法涉及大量區(qū)間內(nèi)積分運算,對差分傳播相移有較高的質(zhì)量要求,因此需要對其進(jìn)行質(zhì)控,該算法流程如圖1所示。
圖1 自適應(yīng)算法流程
自適應(yīng)算法流程可以概括為:對每個體掃數(shù)據(jù)的所有仰角層逐層訂正,每層進(jìn)行逐徑向訂正,將每個待訂正的徑向分為若干區(qū)間,分別計算這些區(qū)間內(nèi)對應(yīng)的衰減系數(shù)α,再根據(jù)徑向內(nèi)所有區(qū)間的衰減系數(shù)對該徑向進(jìn)行訂正,然后完成整層數(shù)據(jù)的訂正,最后再完成整個體掃數(shù)據(jù)的訂正。算法具體步驟如下:
(1)給定衰減系數(shù)α的初始值為0.01。
(2)利用ZH、降水量I和差分傳播相移ΦDP計算衰減率AH。
其中,I為降水量,b為衰減指數(shù)常數(shù),取值0.76~0.84,本文取均值0.8,α為區(qū)間待求衰減系數(shù),ZH為訂正前水平反射率因子(雷達(dá)存儲值)。
(3)根據(jù)KDP和AH的關(guān)系、KDP和ΦDP的關(guān)系估算
(5)衰減系數(shù)α以0.05為步長,遞增,重復(fù)步驟(4),當(dāng)兩者的差值最小時,即得到最佳衰減系數(shù)。
(6)應(yīng)用最佳衰減系數(shù)α進(jìn)行訂正。
因S波段天氣雷達(dá)波長較長,電磁波受雨區(qū)衰減影響較小,因此驗證數(shù)據(jù)集為北京大興站同時刻(兩雷達(dá)探測誤差在1分鐘以內(nèi))同區(qū)域內(nèi)的S波段單偏振天氣雷達(dá)數(shù)據(jù)。通過兩種方法對不同地區(qū)雷達(dá)反射率因子數(shù)據(jù)進(jìn)行訂正發(fā)現(xiàn),固定衰減系數(shù)的AH-KDP法在近雷達(dá)處能夠?qū)Ψ瓷渎室蜃佑休^弱訂正,在遠(yuǎn)離雷達(dá)處效果不佳。AH-KDP法的反射率因子訂正前后距離廓線對比如圖2所示。反射率因子訂正前后對比PPI圖如圖3所示。而自適應(yīng)算法在近、遠(yuǎn)離雷達(dá)處效果都較明顯,能較好地恢復(fù)真實降雨區(qū)數(shù)據(jù)。自適應(yīng)算法的反射率因子訂正前后距離廓線對比如圖4所示。反射率因子訂正前后對比PPI圖如圖5所示。
圖2 AH-KDP法訂正值與探測值及驗證值對比圖
圖3 AH-KDP法衰減訂正效果對比圖
圖4 自適應(yīng)算法訂正值與探測值及驗證值對比圖
圖5 自適應(yīng)算法衰減訂正效果對比圖
為更直觀反映兩種算法的優(yōu)劣,對同種X波段雷達(dá)探測數(shù)據(jù)分別用兩種算法進(jìn)行訂正,并以同時刻(誤差在1分鐘左右)S波段雷達(dá)探測數(shù)據(jù)作驗證,結(jié)果表明自適應(yīng)算法明顯優(yōu)于AH-KDP法,對比圖見圖6。
圖6 AH-KDP法與自適應(yīng)算法效果對比圖
根據(jù)自適應(yīng)算法和AH-KDP法的效果對比,本文最終選擇自適應(yīng)算法對組網(wǎng)前的3部X波段雙偏振天氣雷達(dá)的反射率因子數(shù)據(jù)進(jìn)行第一步衰減訂正處理。
常用的線性插值算法包括:垂直線性插值、垂直水平線性插值和八點線性插值算法。八點線性內(nèi)插(EPI)選取網(wǎng)格點周圍8個相關(guān)聯(lián)的距離庫作為參考點,將參考維度擴(kuò)展至仰角、方位角、徑向上,能夠提升內(nèi)插的可靠性。該算法解釋為選取該網(wǎng)格點對應(yīng)在極坐標(biāo)系下最近的兩個仰角層和方位角(或徑向)上六面體的8個頂點作為參考點[6]。
如圖7所示,p1、p2、p3、p4為位于該網(wǎng)格點p上方仰角層內(nèi)最近的4個距離庫參考點,p5、p6、p7、p8為位于該網(wǎng)格點p下方仰角層內(nèi)最近的4個距離庫參考點;p1、p3、p5、p7 位于同一方位角,p2、p4、p6、p8 位于同一方位角;p1、p2、p5、p6擁有相同的徑向距離庫,p3、p4、p7、p8擁有相同的徑向距離庫。p1~p8對應(yīng)坐標(biāo)為 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對應(yīng)值為V1~V8,待插值點p在空間笛卡兒坐標(biāo)系下的坐標(biāo)為p(x,y,z),需要將其轉(zhuǎn)換至極坐標(biāo)系下的坐標(biāo)p(e,a,r),轉(zhuǎn)換公式如式(7),最后用三個維度的線性插值對待插值點p的值Vp進(jìn)行計算,計算公式如式(8)。
圖7 八點線性插值法示意圖
其中的we、wr、wa分別為各參考點在仰角、徑向距離和方位角上的權(quán)重系數(shù),其計算公式如式(9)~(11)所示。
在得到多部雷達(dá)網(wǎng)格化處理后的雷達(dá)資料后,常見的拼圖算有以下幾種:最臨近法、最大值法、反距離權(quán)重函數(shù)法[7]。反距離權(quán)重函數(shù)法相較于最大值法和最臨近法,綜合考慮了組網(wǎng)所用的多部雷達(dá)反射率因子在同一網(wǎng)格點的貢獻(xiàn),能夠明顯地減少其他兩種方法帶來的色標(biāo)跳變現(xiàn)象,組網(wǎng)圖像更連續(xù)。因此,組網(wǎng)處理部分最終選擇反距離權(quán)重函數(shù)法,其示意圖如圖8所示。
圖8 反距離權(quán)重法示意圖
由圖8可知,3部雷達(dá)形成一個三角形區(qū)域,網(wǎng)格點p的值由該網(wǎng)格點到3部雷達(dá)R1、R2、R3距離決定,最后根據(jù)網(wǎng)格點p到3部雷達(dá)在三維x-y-z空間上的直線距離r1、r2、r3進(jìn)行反距離加權(quán),設(shè) r1、r2、r3大小遞減。3部雷達(dá)的坐標(biāo)為R1(x1,y1,z1)、R2(x2,y2,z2)、R3(x3,y3,z3),p點的坐標(biāo)為p(x,y,z)。3部雷達(dá)在p點的值分別為V1、V2、V3,p點拼圖后的值應(yīng)為Vp。拼圖公式如下:
對2020年8月12日10:00:00北京房山站、北京昌平站、北京順義站X波段雙偏振天氣雷達(dá)的水平反射率因子PPI以八點線性插值方法進(jìn)行3 km高的CAPPI插值與組網(wǎng)處理,其插值結(jié)果、組網(wǎng)拼圖結(jié)果分別如圖9所示。
由圖9可知,EPI法在3個維度(仰角、方位角、距離)上選擇8個參考點,能夠?qū)?部X波段雙偏振天氣雷達(dá)極坐標(biāo)系下的反射率因子資料合適地插值到統(tǒng)一的空間笛卡爾坐標(biāo)系下(CAPPI)。由圖10可知,反距離權(quán)重函數(shù)法基本能將3部雷達(dá)的插值結(jié)果較完整的組網(wǎng)到同一個網(wǎng)格上。至此,完成前期的3部雷達(dá)的反射率因子組網(wǎng)處理。
圖9 三地區(qū)插值處理CAPPI圖
圖10 組網(wǎng)拼圖處理CAPPI圖
預(yù)處理部分主要分為缺空補(bǔ)值和斑點濾波[8-9],關(guān)于反射率因子CAPPI的缺空補(bǔ)值處理,設(shè)計NA×NR窗口濾波函數(shù)[10]對所有網(wǎng)格點上的值進(jìn)行濾波檢測。該窗口濾波函數(shù)可描述為:在每格網(wǎng)格點的第i個經(jīng)度和第j個緯度處的有效值,其經(jīng)度方向上的窗口尺度為NA,其緯度方向上的窗口尺度為NR,當(dāng)窗口內(nèi)的有效數(shù)據(jù)總數(shù)Ti,j在窗口內(nèi)總庫數(shù)占比Pi,j小于閾值Thr1時,則將該庫視為孤立回波[11],并將該庫處的值設(shè)置為無效值NaN,如公式(15)。當(dāng)窗口內(nèi)的有效數(shù)據(jù)總數(shù)Ti,j在窗口內(nèi)總庫數(shù)占比Pi,j大于閾值Thr2時,則將該庫視為缺測值庫,并將該庫處的值從NaN設(shè)置為一個有效值,該有效值F為該窗口內(nèi)所有有效值的均值,如式(16)。
其中閾值Thr1設(shè)置為0.45,閾值Thr2設(shè)置為0.3,對北京房山、昌平、順義地區(qū)經(jīng)過EPI法得到的反射率因子CAPPI進(jìn)行缺空補(bǔ)值、斑點濾波處理前后的對比如圖11所示。
圖11 預(yù)處理前后CAPPI圖
由圖11預(yù)處理前后對比,可知預(yù)處理能夠?qū)吘壙杖敝颠M(jìn)行有效回填,并對組網(wǎng)區(qū)域內(nèi)的斑點雜波進(jìn)行相應(yīng)的去除。
關(guān)于組網(wǎng)反射率因子CAPPI的訂正處理,由于在應(yīng)用反射率因子進(jìn)行CAPPI組網(wǎng)拼圖前,已經(jīng)對3部雷達(dá)的反射率因子以自適應(yīng)算法進(jìn)行訂正處理,但這種處理不夠徹底。因為訂正是基于雷達(dá)電磁波在發(fā)射方向上不斷衰減能量的前提下,運用最佳衰減系數(shù)和差分傳播相移使探測值盡量靠近真實值。因此在反射率因子CAPPI組網(wǎng)之后,也要對網(wǎng)格點上的反射率因子值進(jìn)行適當(dāng)訂正。本文設(shè)計使用區(qū)域檢測算法對CAPPI每個網(wǎng)格點進(jìn)行訂正。該算法整體思路為以S波段雷達(dá)為參考,對拼圖CAPPI和S波段雷達(dá)CAPPI的每個網(wǎng)格點建立維度為N的搜索框,并分別計算其均值Mean_S和Mean_X,當(dāng)前者與后者的差值大于閾值Thr_C時,則判定為偏差嚴(yán)重,該偏差包括系統(tǒng)誤差和衰減,需要進(jìn)行訂正。為防止訂正后的值大于真實值,以盡量靠近真實值的為原則,訂正后的值Vp為當(dāng)前拼圖CAPPI值VX加上Mean_S與Mean_X的差值,其中N取5,閾值Thr_C取2 dBZ。建立的搜索框示意圖如圖12。
圖12 搜索框示意圖
訂正公式如式(17),組網(wǎng)反射率因子CAPPI訂正前后對比效果如圖13。
圖13 訂正處理前后對比圖
由圖13可知,訂正部分主要在39°N~40.3°N、115.5°E~118°E附近,將訂正后的組網(wǎng)反射率因子CAPPI與同時刻同地點的S波段雷達(dá)反射率因子CAPPI進(jìn)行對比,如圖14。
圖14 訂正CAPPI與驗證CAPPI對比圖
由圖 13 可 知,115.5°E~118°E,39°N~40.3°N為主要訂正區(qū)域,區(qū)域檢測算法以待訂區(qū)域內(nèi)同時刻下S波段天氣雷達(dá)探測資料為參考,對每個待訂正的網(wǎng)格點反射率因子建立搜索框,并計算兩個搜索框內(nèi)的區(qū)域均值,當(dāng)兩均值差異較大時,則判定為偏差過大,利用S波段雷達(dá)資料對其進(jìn)行訂正處理。由圖14可知,主要訂正區(qū)域內(nèi)的訂正結(jié)果基本與驗證用S波段雷達(dá)反射率因子CAPPI相吻合,該算法訂正效果表現(xiàn)良好。
在前人研究的基礎(chǔ)上,首先通過分析2種常規(guī)衰減訂正算法的優(yōu)劣,選擇較優(yōu)算法對單部雷達(dá)反射率因子訂正處理;其次運用經(jīng)典八點線性插值算法和反距離權(quán)重函數(shù)法對訂正后的3部X波段雙偏振天氣雷達(dá)反射率因子進(jìn)行網(wǎng)格化處理,最后設(shè)計了區(qū)域檢測算法對組網(wǎng)反射率因子CAPPI進(jìn)行適當(dāng)?shù)挠喺幚?實驗結(jié)果表明,本文設(shè)計的區(qū)域檢測算法在CAPPI網(wǎng)格上表現(xiàn)良好的效果,具有一定的參考價值,但該方法局限于組網(wǎng)區(qū)域內(nèi)需要存在衰減較小的長波長天氣雷達(dá)用以驗證算法效果,極大限制了該方法在業(yè)務(wù)上的應(yīng)用,未來需設(shè)計其他數(shù)值訂正方案來優(yōu)化該算法。