高 霞, 周 曄, 周淑玥, 杜 捷, 王海江

(1.中國(guó)航空工業(yè)集團(tuán)公司雷華電子技術(shù)研究所,江蘇 無(wú)錫 214000;2.成都信息工程大學(xué)電子工程學(xué)院,四川 成都 610225)

0 引言

雷達(dá)電磁波在傳輸過(guò)程中,能量會(huì)被沿途接觸到的各類(lèi)粒子吸收轉(zhuǎn)變成機(jī)械能、內(nèi)能、熱能等[1],電磁波的衰減與波長(zhǎng)成反比。X波段雙偏振天氣雷達(dá)由于其波長(zhǎng)較短,在探測(cè)降雨時(shí)能量會(huì)被降水粒子吸收等,導(dǎo)致能量衰減,在遠(yuǎn)距離處的弱雷達(dá)回波返回時(shí)很難被雷達(dá)接收機(jī)接收,導(dǎo)致遠(yuǎn)距離處信息丟失,嚴(yán)重影響雷達(dá)的探測(cè)精度以及雷達(dá)探測(cè)資料的應(yīng)用。

X波段雙偏振天氣雷達(dá)存儲(chǔ)有差分傳播相移ΦDP,該參數(shù)為距離累積量不受衰減影響[2],因此被廣泛用于反射率因子的衰減訂正。目前較為廣泛的衰減訂正方法是通過(guò)計(jì)算衰減率來(lái)完成反射率因子的訂正,但衰減率大都是給定一個(gè)固定的區(qū)間范圍[3],存在諸多弊端,如不能靈活適應(yīng)不同降水粒子的衰減。因此可以根據(jù)差分傳播相移ΦDP不受衰減的影響來(lái)設(shè)計(jì)自適應(yīng)算法完成衰減訂正工作[4-5]。

首先對(duì)單部雷達(dá)進(jìn)行反射率因子衰減訂正處理,其次將衰減訂正處理后的3部X波段雙偏振天氣雷達(dá)反射率因子插值到統(tǒng)一的CAPPI網(wǎng)格上,再對(duì)組網(wǎng)反射率因子CAPPI進(jìn)行缺空補(bǔ)值和斑點(diǎn)濾波處理,最后設(shè)計(jì)區(qū)域檢測(cè)算法,對(duì)組網(wǎng)反射率因子進(jìn)行訂正。

1 資料選取

3部X波雙偏振天氣雷達(dá)隸屬于北京市氣象局各區(qū)分局的雷達(dá)站點(diǎn),分別是北京房山站雷達(dá)、北京昌平站雷達(dá)、北京順義站雷達(dá),另外用于驗(yàn)證的1部SA波段雷達(dá)為北京大興站雷達(dá)。各站點(diǎn)雷達(dá)的詳細(xì)參數(shù)如表1所示。

表1 雷達(dá)站點(diǎn)及詳細(xì)參數(shù)

研究數(shù)據(jù)來(lái)源于2020年夏季(8-9月)雷達(dá)探測(cè)數(shù)據(jù)(表1),由于夏季降水豐沛,雷達(dá)聯(lián)合組網(wǎng)后效果更佳。

2 單部雷達(dá)反射率因子衰減訂正

因?yàn)樗p是雙程的累積量,故衰減訂正為

2.1 AH-KDP法

對(duì)于衰減率的計(jì)算,Bringi等[5]通過(guò)散射的模擬實(shí)驗(yàn)發(fā)現(xiàn)反射率因子的衰減率AH與差分傳播相移率KDP基本上呈線性關(guān)系,并在實(shí)驗(yàn)的基礎(chǔ)上提出了衰減訂正公式:

其中,α為衰減系數(shù),Testud指出對(duì)于X波段雷達(dá),衰減系數(shù)α在0.139～0.335[3],本文取其均值0.237。

2.2 自適應(yīng)算法

為更好地適應(yīng)不同降雨類(lèi)型,Bringi等[5]對(duì)上節(jié)算法進(jìn)行改進(jìn),命名為自適應(yīng)約束算法。自適應(yīng)約束算法在計(jì)算中使用了差分傳播相移ΦDP和水平反射率因子ZH,并根據(jù)誤差大小不斷調(diào)整衰減系數(shù)α,將其調(diào)整為最佳值。自適應(yīng)算法流程多,迭代次數(shù)大,算法涉及大量區(qū)間內(nèi)積分運(yùn)算,對(duì)差分傳播相移有較高的質(zhì)量要求,因此需要對(duì)其進(jìn)行質(zhì)控,該算法流程如圖1所示。

圖1 自適應(yīng)算法流程

自適應(yīng)算法流程可以概括為:對(duì)每個(gè)體掃數(shù)據(jù)的所有仰角層逐層訂正,每層進(jìn)行逐徑向訂正,將每個(gè)待訂正的徑向分為若干區(qū)間,分別計(jì)算這些區(qū)間內(nèi)對(duì)應(yīng)的衰減系數(shù)α,再根據(jù)徑向內(nèi)所有區(qū)間的衰減系數(shù)對(duì)該徑向進(jìn)行訂正,然后完成整層數(shù)據(jù)的訂正,最后再完成整個(gè)體掃數(shù)據(jù)的訂正。算法具體步驟如下:

(1)給定衰減系數(shù)α的初始值為0.01。

(2)利用ZH、降水量I和差分傳播相移ΦDP計(jì)算衰減率AH。

其中,I為降水量,b為衰減指數(shù)常數(shù),取值0.76～0.84,本文取均值0.8,α為區(qū)間待求衰減系數(shù),ZH為訂正前水平反射率因子(雷達(dá)存儲(chǔ)值)。

(3)根據(jù)KDP和AH的關(guān)系、KDP和ΦDP的關(guān)系估算

(5)衰減系數(shù)α以0.05為步長(zhǎng),遞增,重復(fù)步驟(4),當(dāng)兩者的差值最小時(shí),即得到最佳衰減系數(shù)。

(6)應(yīng)用最佳衰減系數(shù)α進(jìn)行訂正。

2.3 兩種衰減訂正方法效果對(duì)比分析

因S波段天氣雷達(dá)波長(zhǎng)較長(zhǎng),電磁波受雨區(qū)衰減影響較小,因此驗(yàn)證數(shù)據(jù)集為北京大興站同時(shí)刻(兩雷達(dá)探測(cè)誤差在1分鐘以內(nèi))同區(qū)域內(nèi)的S波段單偏振天氣雷達(dá)數(shù)據(jù)。通過(guò)兩種方法對(duì)不同地區(qū)雷達(dá)反射率因子數(shù)據(jù)進(jìn)行訂正發(fā)現(xiàn),固定衰減系數(shù)的AH-KDP法在近雷達(dá)處能夠?qū)Ψ瓷渎室蜃佑休^弱訂正,在遠(yuǎn)離雷達(dá)處效果不佳。AH-KDP法的反射率因子訂正前后距離廓線對(duì)比如圖2所示。反射率因子訂正前后對(duì)比PPI圖如圖3所示。而自適應(yīng)算法在近、遠(yuǎn)離雷達(dá)處效果都較明顯,能較好地恢復(fù)真實(shí)降雨區(qū)數(shù)據(jù)。自適應(yīng)算法的反射率因子訂正前后距離廓線對(duì)比如圖4所示。反射率因子訂正前后對(duì)比PPI圖如圖5所示。

圖2 AH-KDP法訂正值與探測(cè)值及驗(yàn)證值對(duì)比圖

圖3 AH-KDP法衰減訂正效果對(duì)比圖

圖4 自適應(yīng)算法訂正值與探測(cè)值及驗(yàn)證值對(duì)比圖

圖5 自適應(yīng)算法衰減訂正效果對(duì)比圖

為更直觀反映兩種算法的優(yōu)劣,對(duì)同種X波段雷達(dá)探測(cè)數(shù)據(jù)分別用兩種算法進(jìn)行訂正,并以同時(shí)刻(誤差在1分鐘左右)S波段雷達(dá)探測(cè)數(shù)據(jù)作驗(yàn)證,結(jié)果表明自適應(yīng)算法明顯優(yōu)于AH-KDP法,對(duì)比圖見(jiàn)圖6。

圖6 AH-KDP法與自適應(yīng)算法效果對(duì)比圖

根據(jù)自適應(yīng)算法和AH-KDP法的效果對(duì)比,本文最終選擇自適應(yīng)算法對(duì)組網(wǎng)前的3部X波段雙偏振天氣雷達(dá)的反射率因子數(shù)據(jù)進(jìn)行第一步衰減訂正處理。

3 網(wǎng)格化處理

3.1 插值處理

常用的線性插值算法包括:垂直線性插值、垂直水平線性插值和八點(diǎn)線性插值算法。八點(diǎn)線性內(nèi)插(EPI)選取網(wǎng)格點(diǎn)周?chē)?個(gè)相關(guān)聯(lián)的距離庫(kù)作為參考點(diǎn),將參考維度擴(kuò)展至仰角、方位角、徑向上,能夠提升內(nèi)插的可靠性。該算法解釋為選取該網(wǎng)格點(diǎn)對(duì)應(yīng)在極坐標(biāo)系下最近的兩個(gè)仰角層和方位角(或徑向)上六面體的8個(gè)頂點(diǎn)作為參考點(diǎn)[6]。

如圖7所示,p1、p2、p3、p4為位于該網(wǎng)格點(diǎn)p上方仰角層內(nèi)最近的4個(gè)距離庫(kù)參考點(diǎn),p5、p6、p7、p8為位于該網(wǎng)格點(diǎn)p下方仰角層內(nèi)最近的4個(gè)距離庫(kù)參考點(diǎn);p1、p3、p5、p7 位于同一方位角,p2、p4、p6、p8 位于同一方位角;p1、p2、p5、p6擁有相同的徑向距離庫(kù),p3、p4、p7、p8擁有相同的徑向距離庫(kù)。p1～p8對(duì)應(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對(duì)應(yīng)值為V1～V8,待插值點(diǎn)p在空間笛卡兒坐標(biāo)系下的坐標(biāo)為p(x,y,z),需要將其轉(zhuǎn)換至極坐標(biāo)系下的坐標(biāo)p(e,a,r),轉(zhuǎn)換公式如式(7),最后用三個(gè)維度的線性插值對(duì)待插值點(diǎn)p的值Vp進(jìn)行計(jì)算,計(jì)算公式如式(8)。

圖7 八點(diǎn)線性插值法示意圖

其中的we、wr、wa分別為各參考點(diǎn)在仰角、徑向距離和方位角上的權(quán)重系數(shù),其計(jì)算公式如式(9)～(11)所示。

3.2 組網(wǎng)拼圖處理

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

圖8 反距離權(quán)重法示意圖

由圖8可知,3部雷達(dá)形成一個(gè)三角形區(qū)域,網(wǎng)格點(diǎn)p的值由該網(wǎng)格點(diǎn)到3部雷達(dá)R1、R2、R3距離決定,最后根據(jù)網(wǎng)格點(diǎn)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點(diǎn)的坐標(biāo)為p(x,y,z)。3部雷達(dá)在p點(diǎn)的值分別為V1、V2、V3,p點(diǎn)拼圖后的值應(yīng)為Vp。拼圖公式如下:

3.3 網(wǎng)格化處理結(jié)果

對(duì)2020年8月12日10:00:00北京房山站、北京昌平站、北京順義站X波段雙偏振天氣雷達(dá)的水平反射率因子PPI以八點(diǎn)線性插值方法進(jìn)行3 km高的CAPPI插值與組網(wǎng)處理,其插值結(jié)果、組網(wǎng)拼圖結(jié)果分別如圖9所示。

由圖9可知,EPI法在3個(gè)維度(仰角、方位角、距離)上選擇8個(gè)參考點(diǎn),能夠?qū)?部X波段雙偏振天氣雷達(dá)極坐標(biāo)系下的反射率因子資料合適地插值到統(tǒng)一的空間笛卡爾坐標(biāo)系下(CAPPI)。由圖10可知,反距離權(quán)重函數(shù)法基本能將3部雷達(dá)的插值結(jié)果較完整的組網(wǎng)到同一個(gè)網(wǎng)格上。至此,完成前期的3部雷達(dá)的反射率因子組網(wǎng)處理。

圖9 三地區(qū)插值處理CAPPI圖

圖10 組網(wǎng)拼圖處理CAPPI圖

4 組網(wǎng)反射率因子訂正

4.1 預(yù)處理

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

其中閾值Thr1設(shè)置為0.45,閾值Thr2設(shè)置為0.3,對(duì)北京房山、昌平、順義地區(qū)經(jīng)過(guò)EPI法得到的反射率因子CAPPI進(jìn)行缺空補(bǔ)值、斑點(diǎn)濾波處理前后的對(duì)比如圖11所示。

圖11 預(yù)處理前后CAPPI圖

由圖11預(yù)處理前后對(duì)比,可知預(yù)處理能夠?qū)吘壙杖敝颠M(jìn)行有效回填,并對(duì)組網(wǎng)區(qū)域內(nèi)的斑點(diǎn)雜波進(jìn)行相應(yīng)的去除。

4.2 區(qū)域檢測(cè)算法

關(guān)于組網(wǎng)反射率因子CAPPI的訂正處理,由于在應(yīng)用反射率因子進(jìn)行CAPPI組網(wǎng)拼圖前,已經(jīng)對(duì)3部雷達(dá)的反射率因子以自適應(yīng)算法進(jìn)行訂正處理,但這種處理不夠徹底。因?yàn)橛喺腔诶走_(dá)電磁波在發(fā)射方向上不斷衰減能量的前提下,運(yùn)用最佳衰減系數(shù)和差分傳播相移使探測(cè)值盡量靠近真實(shí)值。因此在反射率因子CAPPI組網(wǎng)之后,也要對(duì)網(wǎng)格點(diǎn)上的反射率因子值進(jìn)行適當(dāng)訂正。本文設(shè)計(jì)使用區(qū)域檢測(cè)算法對(duì)CAPPI每個(gè)網(wǎng)格點(diǎn)進(jìn)行訂正。該算法整體思路為以S波段雷達(dá)為參考,對(duì)拼圖CAPPI和S波段雷達(dá)CAPPI的每個(gè)網(wǎng)格點(diǎn)建立維度為N的搜索框,并分別計(jì)算其均值Mean_S和Mean_X,當(dāng)前者與后者的差值大于閾值Thr_C時(shí),則判定為偏差嚴(yán)重,該偏差包括系統(tǒng)誤差和衰減,需要進(jìn)行訂正。為防止訂正后的值大于真實(shí)值,以盡量靠近真實(shí)值的為原則,訂正后的值Vp為當(dāng)前拼圖CAPPI值VX加上Mean_S與Mean_X的差值,其中N取5,閾值Thr_C取2 dBZ。建立的搜索框示意圖如圖12。

圖12 搜索框示意圖

訂正公式如式(17),組網(wǎng)反射率因子CAPPI訂正前后對(duì)比效果如圖13。

圖13 訂正處理前后對(duì)比圖

由圖13可知,訂正部分主要在39°N～40.3°N、115.5°E～118°E附近,將訂正后的組網(wǎng)反射率因子CAPPI與同時(shí)刻同地點(diǎn)的S波段雷達(dá)反射率因子CAPPI進(jìn)行對(duì)比,如圖14。

圖14 訂正CAPPI與驗(yàn)證CAPPI對(duì)比圖

4.3 訂正結(jié)果分析

由圖 13 可 知,115.5°E～118°E,39°N～40.3°N為主要訂正區(qū)域,區(qū)域檢測(cè)算法以待訂區(qū)域內(nèi)同時(shí)刻下S波段天氣雷達(dá)探測(cè)資料為參考,對(duì)每個(gè)待訂正的網(wǎng)格點(diǎn)反射率因子建立搜索框,并計(jì)算兩個(gè)搜索框內(nèi)的區(qū)域均值,當(dāng)兩均值差異較大時(shí),則判定為偏差過(guò)大,利用S波段雷達(dá)資料對(duì)其進(jìn)行訂正處理。由圖14可知,主要訂正區(qū)域內(nèi)的訂正結(jié)果基本與驗(yàn)證用S波段雷達(dá)反射率因子CAPPI相吻合,該算法訂正效果表現(xiàn)良好。

5 結(jié)束語(yǔ)

在前人研究的基礎(chǔ)上,首先通過(guò)分析2種常規(guī)衰減訂正算法的優(yōu)劣,選擇較優(yōu)算法對(duì)單部雷達(dá)反射率因子訂正處理;其次運(yùn)用經(jīng)典八點(diǎn)線性插值算法和反距離權(quán)重函數(shù)法對(duì)訂正后的3部X波段雙偏振天氣雷達(dá)反射率因子進(jìn)行網(wǎng)格化處理,最后設(shè)計(jì)了區(qū)域檢測(cè)算法對(duì)組網(wǎng)反射率因子CAPPI進(jìn)行適當(dāng)?shù)挠喺幚?實(shí)驗(yàn)結(jié)果表明,本文設(shè)計(jì)的區(qū)域檢測(cè)算法在CAPPI網(wǎng)格上表現(xiàn)良好的效果,具有一定的參考價(jià)值,但該方法局限于組網(wǎng)區(qū)域內(nèi)需要存在衰減較小的長(zhǎng)波長(zhǎng)天氣雷達(dá)用以驗(yàn)證算法效果,極大限制了該方法在業(yè)務(wù)上的應(yīng)用,未來(lái)需設(shè)計(jì)其他數(shù)值訂正方案來(lái)優(yōu)化該算法。