雙偏振雷達差分傳播相移質(zhì)量控制的一種新方法——線性擬合-遞推法

孫 躍，肖 輝，馮 亮，楊慧玲，陳 程，周筠王君

(1.成都信息工程大學大氣科學學院，四川成都610225;2.中國科學院大氣物理研究所云降水物理與強風暴重點實驗室，北京100029)

差分傳播相移ΦDP是雙偏振雷達的一個重要的偏振參量，其不僅可以被用于X-band雙偏振雷達的衰減訂正[1－3]當中，由 ΦDP得到訂正的差分傳播相移率KDP，還可以被應用到雙偏振雷達降水估算[4]、相態(tài)識別[5]當中。但是，ΦDP原始信號中相對于上述用途而言包含了由雷達本身造成的系統(tǒng)噪聲和地物可能引起的觀測噪聲，因此在使用之前需要先對其進行質(zhì)量控制。對于ΦDP的質(zhì)量控制，國內(nèi)外學者分別采用過滑動平均濾波[6－7]、中值濾波[6－8]、Kalman 濾波[2，6，8]、小波濾波[6，8－9]、迭代濾波[10]、分類處理法[7]、FIR 低通濾波[6，8]等多種算法。綜合上述相關(guān)文獻所列舉的濾波個例看，當質(zhì)量控制目標對應的是連續(xù)的云體，且差分傳播相移ΦDP中只存在零星地物噪聲和幅度變化不大的系統(tǒng)噪聲時，用上述算法便能取得較為理想的質(zhì)量控制結(jié)果。但是，當掃描徑向上存在不只一個云體，尤其是對雷暴群掃描時，ΦDP在云體的邊緣位置會產(chǎn)生異常的震蕩且在云體間隔的位置不連續(xù)，或者當觀測噪聲較大以及系統(tǒng)噪聲較大時，使用上述方法均不能得到理想的ΦDP質(zhì)量控制效果，其表現(xiàn)是大幅度雜波只被控制成了小幅度雜波，而沒有完全被濾除。這主要是因為，上述方法大都無法在算法層面將有利用價值的ΦDP變化信息與脈沖式的雜波完全區(qū)分開。杜牧云等[7]提出分類處理法雖然對此問題有了一定的改進，但是其剔除大幅度雜波的算法首先是基于連續(xù)性檢查和大量歷史資料統(tǒng)計得到的，當徑向數(shù)據(jù)不連續(xù)或者大幅度雜波連續(xù)出現(xiàn)時，基于統(tǒng)計意義上的連續(xù)性檢查可能無法正確地從中分離出有價值的ΦDP信號。因此，進一步研究更實用的ΦDP質(zhì)量控制算法是有必要的。

值得注意的是，以往的ΦDP質(zhì)量控制方法均采用一般性的信號處理算法，但ΦDP有別于一般模擬信號/數(shù)字信號的特點。簡而言之，在ΦDP信號中對氣象研究與應用有重要價值的信息不是波動信號或脈沖信號，多數(shù)情況下也不能用單一的解析方式描述。從上述研究中成功的濾波個例[2，6－8，10]來看，有使用價值的ΦDP信號在其掃描經(jīng)向上更像是一個導數(shù)不為負的連續(xù)分段函數(shù)，而這種描述也非常符合ΦDP信號的定義。作為一個距離積分量，ΦDP在粒子直徑較大的探測距離段上增長較快，而在粒子直徑較小的距離段上增長緩慢，甚至有小的波動。從雙偏振雷達的衰減訂正、雨量估算、相態(tài)識別等用途看，也需要使ΦDP曲線的斜率盡量不為負。同時，在以往成功的ΦDP質(zhì)量控制個例中，ΦDP的連續(xù)上升段呈現(xiàn)出明顯的線性增長特征。因此，ΦDP的質(zhì)量控制可以由一個單純的降噪濾波問題轉(zhuǎn)化為一組導數(shù)不為負的連續(xù)線段的確定問題。基于以上論述，設(shè)計了一種針對ΦDP特點的“線性擬合-遞推法”的質(zhì)量控制方法。

1 線性擬合-遞推法介紹

為了說明線性擬合-遞推法的具體應用，下面以中國科學院大氣物理研究所車載X波段雙偏振雷達2008年6月23日在首都機場觀測的多單體強對流天氣為例[11]，選取14:59 分方位角 296°RHI掃描 0.34°仰角的雷達水平偏振反射率ZH與差分傳播相移ΦDP作為研究對象(圖1a)(徑向庫長150 m)。從圖1(a)的ZH曲線可以看出，在掃描線上前400個距離庫上有兩個主要云體對應的兩段回波，而兩段回波并非是連續(xù)的。因此需要先將對應兩段的ΦDP進行拼接(合并)，即填補中間不連續(xù)的數(shù)據(jù)段，再按照介紹的方法對ΦDP進行質(zhì)量控制。具體步驟如下:

(1)滑動線性平滑

采用的滑動線性平滑的原理類似于經(jīng)典的多項式平滑濾波，如五點二次濾波、五點三次濾波[12]，即在滑動窗口D范圍內(nèi)通過最小二乘法求解該段的線性擬合系數(shù)aj與bj，進而求得點j處的線性平滑結(jié)果 ^yj。

其中xj代表第j個距離庫距雷達原點的距離，系數(shù)aj與bj的計算公式即在滑動窗口內(nèi)利用典型的最小二乘法[13]求解一次多項式系數(shù)的公式

xl與yl分別代表第l個點距雷達原點的距離和觀測值。

與一般多項式濾波的不同之處是采用的多項式為一次多項式且滑動窗口D更長，以求得原ΦDP信號各點的線性特征估計值(圖1b)，并以每次線性擬合的殘差平方和RSS的算式定義該j點ΦDP的線性擬合偏差σj(圖1c)，該值越小表明該點符合線性變化特征的程度越好，同時記錄下該點線性擬合的斜率k(圖1d)。單從線性滑動平滑的結(jié)果看，其平滑出的圖形與一般性的滑動平均類似，即波動只在幅度上存在不同程度的減小而并沒有消失，但滑動線性平滑過程中得到的線性擬合偏差和斜率可以在下面的步驟中為遞推求得目標點提供依據(jù)。

其中，σj為線性擬合偏差，D為滑動窗口，D=21。RSS為每次線性擬合的殘差平方和。

圖1 D=21時的線性滑動平滑

(2)對線性擬合偏差進行排序

從經(jīng)驗的角度來說，針對圖1(b)中的原始信號，一般先確定圖中波動幅度較小且線性特征明顯的數(shù)據(jù)段(如虛線框確定的兩段位置)，再依據(jù)波動幅度和斜率不為負的原則依次確定其他數(shù)據(jù)段。對比圖1(b)和圖1(c)中虛線框位置的對應關(guān)系，可以將“先確定波動幅度較小且線性特征明顯的數(shù)據(jù)段”這一目標轉(zhuǎn)化為“先確定線性擬合偏差小的點”。這就需要將上一步驟中得到的線性擬合偏差序列進行升序排序，并求得排序索引，因為希望先確定噪聲小且線性趨勢明顯的數(shù)據(jù)段，再依據(jù)“ΦDP質(zhì)量控制后斜率不明顯為負”的實用目標去確定噪聲大且線性趨勢不明顯的數(shù)據(jù)段，而不是像一般性的滑動平均、中值濾波、低通濾波等方法那樣按照距離庫的探測順序確定質(zhì)量控制后的最終計算值，因為上述這些按照探測順序計算的方法顯然沒有一種機制是專門保留正斜率而抑制負斜率的，即無法考慮后未計算的遠端距離庫是否還有更可信的但是比近端已經(jīng)確定的數(shù)據(jù)還要顯著低的數(shù)據(jù)。因此對于反應噪聲和線性趨勢的線性擬合偏差σj進行排序并按照重新排序的順序進行逐庫遞推計算是非常必要的。

(3)依據(jù)排序結(jié)果遞推目標點

為了確定導數(shù)不為負的連續(xù)線段，給出逐點遞推判斷的規(guī)則為

依據(jù)上述步驟求得排序索引，按照式(1)的規(guī)則逐點遞推判斷該點是否應該為最終求得的點，其中Yi代表將^yj升序排序所得的n個排序索引中第i個索引對應的距離庫上線性滑動平滑值，Yleft代表Yi左側(cè)第一個已被賦值的數(shù)據(jù)值，Yright代表Yi右側(cè)第一個已被賦值的數(shù)據(jù)值，－999和999代表初始的左右邊界。ξ代表松弛變量或稱損失函數(shù)，表征遞推過程中容忍波動的程度。Y0和Yn+1分別代表遞推起始的左右邊界，k代表線性擬合斜率。在每次遞推中，Yleft和Yright中至少有一個會改變，以此判斷新的Yi是否應該屬于最終求得的“導數(shù)不為負的連續(xù)線段”上的點。表1給出D=21、ξ=1時，前15次遞推的過程列表結(jié)果。

表1 部分遞推過程舉例

為了結(jié)合遞推規(guī)則說明參數(shù)ξ與D在該方法中的選擇依據(jù)，下面首先給出D=21時ξ分別為0、0.5、1、10時的遞推結(jié)果(圖2)，可以看出，隨著ξ的增大，遞推判斷后得到的點的數(shù)量增多，但是，波動的幅度也隨之增大，尤其是當ξ=10時，150距離庫附近開始出現(xiàn)明顯的導數(shù)值正負波動，并出現(xiàn)導數(shù)為負值的點，這時遞推過程逐漸泛化為滑動線性平滑;而ξ=0遞推出的符合條件的點過少，最終以這些點做線性插值可能不具有代表性。因此，對“盡可能抑制波動及負斜率”與“盡可能多的反應數(shù)據(jù)細節(jié)特征”兩個目標折中，經(jīng)過一系列實踐，選取ξ=1。但ξ究竟取什么值才具有絕對的普適性，是一個很難論證的問題。不過此處可以設(shè)計成按照經(jīng)驗給定一個ξ初值，當極少數(shù)計算達不到某些特定要求時，迭代增加1到2此ξ即可。

圖2 取不同松弛函數(shù)ξ值的遞推結(jié)果

而對于D的選取依據(jù)，總的來說其不宜過小或過大。因為當D過小時，某些云體邊緣帶有趨勢性的異常噪聲可能會被錯誤的識別為線性趨勢，從而導致線性遞推錯誤。圖3以另一條掃描徑向為例，給出了ξ=1，D分別為11和21的情況，可以看到，700至800距離庫間存在一處有連續(xù)上升趨勢的異常噪聲，過短的滑動窗口會使該段計算出的線性擬合偏差很小，從而錯誤的將該段識別為正常的線性趨勢，這將影響它前后位置的遞歸過程，譬如在該例中，700距離庫之前的數(shù)據(jù)點由于線性平滑后的值以及線性擬合斜率均比該處大，按照遞推關(guān)系將不被保留，從而導致數(shù)據(jù)斷裂(圖3a)。當D=21時，滑動窗口足夠大到貫穿這個云體邊緣的異常噪聲區(qū)，便不存在上述問題(圖3b)。而為了盡可能多的保留數(shù)據(jù)趨勢變化的細節(jié)，D也不宜過大。保守起見，選取D=21(3.15 km)以保證這個滑動窗口能貫穿大部分云體邊緣的噪聲異常區(qū)域。

圖3 不同滑動窗口D對遞推結(jié)果的影響

(4)對遞推得到的目標點以外的點進行線性插值

將給定ξ時得到的遞推結(jié)果中的空缺值進行線性插值，即可得到完整連續(xù)的ΦDP質(zhì)量控制結(jié)果。圖4(a)為D=21，ξ=1時經(jīng)過最終質(zhì)量控制得到的ΦDP結(jié)果與原信號的比較。由圖4(a)可見，質(zhì)量控制后ΦDP圖形基本濾去了所有噪聲，并盡可能地反映了原始信號中的線性變化趨勢，且沒有明顯的負斜率波動，實現(xiàn)了預設(shè)的“一組導數(shù)不為負的連續(xù)線段的確定”的目標。圖4(b)為將質(zhì)量控制后的數(shù)據(jù)按照最初拼合數(shù)據(jù)段填補空白的方式再展開，以和ZH進行對比。從圖中可以看到，在整個距離徑向上，當雷達波束經(jīng)過連續(xù)的ZH大值的較強回波區(qū)時，ΦDP的增長是非常顯著的，這符合ΦDP變化的基本原理。

圖4 D=21，ξ=1時的質(zhì)量控制結(jié)果、原始ΦDP和ZH

為了體現(xiàn)上述質(zhì)量控制結(jié)果的優(yōu)越性，下面分別列舉滑動平均、中值濾波、FIR低通濾波3種具有簡單參數(shù)的質(zhì)量控制方法進行對比(圖5)。從圖中可以看到，這些傳統(tǒng)的濾波方法對于云區(qū)內(nèi)幅度不大且較為規(guī)則的噪聲或者零星噪聲都起到了一定的濾除效果，但在云體邊緣仍然有異常噪聲沒有被濾除。而參數(shù)更加復雜的濾波方法如卡爾曼濾波和小波濾波雖然在此無法通過個例來說明其質(zhì)量控制的必然結(jié)果，但這些傳統(tǒng)的信號處理方法對ΦDP質(zhì)量控制的局限性是可以預見的。

圖5 其他質(zhì)量控制方法得到的結(jié)果

2 在衰減訂正中的應用與分析

衰減訂正是差分傳播相移(ΦDP)的主要用途之一。X-band雙偏振雷達由于波長較短，雷達波經(jīng)過強降水區(qū)后反射率會出現(xiàn)明顯的衰減，因此，需要對雷達探測到的反射率進行衰減訂正。從雙偏振雷達的原理上看，ΦDP是不因探測粒子尺度的大小而減弱的，因此，許多雙偏振雷達反射率衰減訂正方法[1，2，9]都是基于ΦDP進行的。這里列舉一種雙偏振雷達衰減訂正的自適應算法[1，14－15]，該方法通過計算每個距離庫上的衰減量AH并積分，再與觀測值Z'h相加得到衰減訂正后的反射率 Zh，10lg[Zh(r)] =10lg[Z'h(r)] +，衰減量 AH為

其中計算參數(shù)I

式中，r0、r1是掃描徑向上雨區(qū)的起止距離庫，ΔΦDP是r1、r0位置上的ΦDP之差。所謂“自適應”算法指的是通過枚舉參數(shù)α和b可能的取值，得到許多組AH(r，α，b)并利用其重構(gòu)出 ΦcalDP(r，α，b):

再通過雨區(qū)徑向范圍內(nèi) ΦcalDP(r，α，b)與 ΦDP(r)的差值最小作為約束條件確定最優(yōu)的參數(shù)α和b。

以第一節(jié)中ΦDP質(zhì)量控制個例的掃描徑向數(shù)據(jù)為例，參考畢永恒等[15]同樣是針對X波段雙偏振雷達的衰減訂正研究，選取α取值范圍0.13～0.35，b值固定為0.8，對于第242至357距離庫上的雨區(qū)，采用上述衰減訂正方法對反射率ZH進行衰減訂正，得到圖6結(jié)果。從圖中可以看到，在300距離庫以后ΦDP增長明顯的部分，ZH在被訂正以后也明顯增大了。

圖6 采用經(jīng)過質(zhì)量控制后ΦDP進行雷達回波強度的衰減訂正前后對比

為了進一步說明新的ΦDP質(zhì)量控制方法在衰減訂正中的優(yōu)勢，下面列舉一個典型的在ΦDP質(zhì)量控制效果并不理想的情況下衰減訂正失敗的例子。圖7為采用傳統(tǒng)的滑動平均濾波方法(這里取13點滑動平均)濾波后的ΦDP代入與圖6相同的自適應算法得到的衰減訂正結(jié)果。由圖7可見，訂正后的ZH反而減小了，這是因為以滑動平均做ΦDP質(zhì)量控制并不能有效地濾除ΦDP在云體邊緣異常振幅的噪聲，因而在計算242到357距離庫雨區(qū)內(nèi)的 ΦDP變化時，ΔΦDP的符號為負，這將使式4不能正確地反映衰減量，從而完全影響了該算法的衰減訂正結(jié)果。針對這一情況，雖然可以通過改變計算雨區(qū)范圍的大小減小這種影響，但只要質(zhì)量控制后的ΦDP仍然存在較大量級的噪聲或者較多負斜率的情況，這種問題就無法完全避免。

圖7 采用滑動平均后ΦDP根據(jù)采用自適應算法進行衰減訂正

同時還應注意到，其他與ΦDP有關(guān)的衰減訂正方法，比如以 KDP計算 AH的方法[1，9]、利用 ΦDP與總衰減量成線性關(guān)系的方法[2]等，對ΦDP的質(zhì)量控制結(jié)果是同樣敏感的。因此，能否進行有效的ΦDP質(zhì)量控制是雙偏振氣象雷達衰減訂正結(jié)果好壞的重要因素。

另外，文中所介紹的線性擬合-遞推法在其他觀測個例中的ΦDP質(zhì)量控制與ZH衰減訂正工作中也取得了較好的結(jié)果，在這里就不再全面介紹了。

3 小結(jié)與討論

提出了一種雙參數(shù)的線性擬合-遞推法來解決雙偏振雷達差分傳播相移ΦDP的質(zhì)量控制問題，并通過多個掃描徑向上的ΦDP質(zhì)量控制應用說明該新算法的有效性，結(jié)論如下:

(1)提出了ΦDP質(zhì)量控制的線性擬合-遞推法。該方法通過滑動線性平滑得到每個距離庫上ΦDP的線性擬合偏差與線性擬合斜率，并根據(jù)ΦDP斜率非負性的質(zhì)量控制目標構(gòu)造遞推關(guān)系式，依據(jù)線性擬合偏差的排序索引逐點遞推，再經(jīng)過線性插值，最終得到斜率基本非負、可濾除所有噪聲且保留有價值變化趨勢的ΦDP質(zhì)量控制結(jié)果。

(2)多個實例應用表明，該線性擬合-遞推法均取得了較好的ΦDP質(zhì)量控制結(jié)果，不同量級幅度的噪聲基本得到了有效的濾除，基本抑制了ΦDP的負斜率，且體現(xiàn)了ΦDP增長的部分與ZH連續(xù)大值區(qū)的對應關(guān)系，說明該方法物理意義明確，而且操作簡單。但對于該方法以及兩個參數(shù)D和ξ的普適性，應在今后的工作中開展更多的實踐與討論。

[1] Park S G，V N Bringi，V Chandrasekar，et al.Correction of Radar Reflectivity for Rain Attenuation at X Band.Part I:Theoretical and Empirical Basis[J].J.Atmos.Oceanic Techno.，2005，22:1621－1632.

[2] 何宇翔，呂達仁，肖輝，等.X波段雙線極化雷達反射率的衰減訂正[J].大氣科學，2009，33(5):1027－1037.

[3] 何宇翔，呂達仁，肖輝.X波段雙線極化雷達差分反射率的衰減訂正[J].高原氣象，2009，28(3):607－616.

[4] 劉黎平，葛潤生，張沛源.雙線偏振多普勒天氣雷達遙測降水強度和液態(tài)含水量的方法和精度研究[J].大氣科學，2002，26(5):709－719.

[5] 郭鳳霞，馬學謙，王濤，等.基于X波段雙線偏振天氣雷達的雷暴云粒子識別[J].氣象學報，2014，72(6):1－15.

[6] 魏慶，胡志群，劉黎平.雙偏振雷達差分傳播相移的五種濾波方法對比分析[J]，成都信息工程學院學報，2014，29(6):596－602.

[7] 杜牧云，劉黎平，胡志群，等.雙線偏振雷達差分傳播相移的質(zhì)量控制[J]，應用氣象學報，2012，23(6):710－720.

[8] Hu Z Q，Liu L P.Applications of wavelet analysis in differential propagation phase shift data denoising[J].Adv.Atmos.Sci.，2014，31(4):824－834.

[9] 杜牧云，劉黎平，胡志群，等.雙線偏振雷達差分傳播相移的小波濾波初探[J].暴雨災害，2012，31(3):248－254.

[10] Hubbert J，Bringi V N.An iterative filtering technique for the analysis of copolar differential phase and dual-frequency radar measurements[J].J.Atmos.Oceanic Technol.，1995，12(3):643－648.

[11] 陳程，肖輝，馮亮，等.北京地區(qū)強冰雹風暴的雙偏振特征觀測分析[J].安徽農(nóng)業(yè)科學，2014，42(22):7511－7515，7516.

[12] 徐士良，馬爾妮.常用算法程序集(C/C++描述)(第五版)[M].北京:清華大學出版社，2013:591.

[13] 賈俊平.統(tǒng)計學(第四版)[M].北京:中國人民大學出版社，2011:170－173.

[14] Bringi V N，Keenan T D，Chandrasekar V.Correcting C-band radar reflectivity and differential reflectivity data for rain attenuation:a self-consistent method with constraints[J].IEEE Transactions on Geosience and Remote Sensing，2001，39(9):1906－1915.

[15] 畢永恒，劉錦麗，段樹，等.X波段雙線偏振氣象雷達反射率的衰減訂正[J].大氣科學，2012，36(3):495－506.