一種基于譜連續(xù)的極化多普勒氣象雷達(dá)信號重構(gòu)方法

摘" 要：

譜極化濾波技術(shù)廣泛應(yīng)用于氣象雷達(dá)雜波抑制中，但是當(dāng)?shù)仉s波與降雨目標(biāo)混疊時(shí)，這些方法在濾除雜波的同時(shí)會(huì)損失氣象信息。針對這一問題，提出一種基于譜連續(xù)的極化多普勒氣象雷達(dá)信號重構(gòu)方法。該方法利用氣象回波和雜波在距離多普勒圖上的特征差異實(shí)現(xiàn)氣象目標(biāo)的保留和雜波的濾除。具體是利用氣象目標(biāo)的多普勒速度和譜寬的空間連續(xù)性，采用多項(xiàng)式擬合和高斯擬合的方法補(bǔ)償降雨區(qū)域。業(yè)務(wù)氣象雷達(dá)實(shí)測數(shù)據(jù)驗(yàn)證了所提方法的有效性。

關(guān)鍵詞：

極化多普勒氣象雷達(dá); 雜波抑制; 信號恢復(fù); 譜連續(xù)

中圖分類號：

TN 95

文獻(xiàn)標(biāo)志碼： A""" DOI：10.12305/j.issn.1001-506X.2024.07.12

A method for signal reconstruction of polarization-Doppler weather

radar based on spectral continuity

AN Mengyun， YIN Jiapeng*， HUANG Jiankai， LI Yongzhen， WANG Xuesong

（School of Electronic Science， National University of Defense Technology， Changsha 410073， China）

Abstract：

Spectral polarization filtering technology is widely used in weather radar clutter suppression， but when ground clutter is mixed with precipitation， these methods will lose precipitation while filtering out clutter. To address this issue， a polarization-Doppler weather radar signal reconstruction method based on spectral continuity is proposed. This method utilizes the characteristic differences between precipitation and clutter on the range-Doppler map to achieve the preservation of precipitation and the filtering of clutter. Specifically， the spatial continuity of Doppler velocity and spectral width of precipitation is utilized to compensate for precipitation using polynomial fitting and Gaussian fitting method. The effectiveness of the proposed method is verified by the measured data of the operational weather radar.

Keywords：

polarization-Doppler weather radar; clutter suppression; signal recovery; spectral continuity

0" 引" 言

氣象雷達(dá)是進(jìn)行高時(shí)間和空間分辨率大氣觀測的有效工具［12］。利用極化多普勒氣象雷達(dá)獲取的信息可以進(jìn)行水凝物分類［3］、定量降水估計(jì)［26］、雨滴譜反演［3］、短臨預(yù)報(bào)［7］。然而，氣象雷達(dá)獲得的數(shù)據(jù)經(jīng)常受到各種雜波、噪聲［8］以及無意的電磁干擾［910］的影響，這會(huì)影響對氣象目標(biāo)的觀測、分類、量化和預(yù)報(bào)［3］。因此，探索和研究有效的雜波抑制技術(shù)尤為重要。

國內(nèi)外學(xué)者提出了多種雜波抑制方法來解決氣象雷達(dá)雜波干擾問題。根據(jù)多普勒速度的不同，可以將影響氣象雷達(dá)數(shù)據(jù)質(zhì)量的雜波分為靜態(tài)的地物雜波和動(dòng)態(tài)雜波。地物對雷達(dá)有兩種類型的干擾，一種是來自固定目標(biāo)（如建筑、樹木）正常傳播（normal propagation， NP）時(shí)的雜波，另一種是與大氣條件有關(guān)的電磁波（如超折射）異常傳播（anomalous propagation， AP）時(shí)的地物雜波，NP雜波和AP雜波具有類似的回波特征［3］。通常來說，地雜波數(shù)據(jù)可通過晴空條件下得到的靜態(tài)雜波圖來識別［11］。然而，由于AP雜波的存在，地物雜波圖會(huì)隨著天氣條件的不同而改變［3］，因此需要通過濾波算法來濾除地雜波。傳統(tǒng)的方法是用以零速度為中心的凹口濾波器來濾除地雜波［12］，但是該方法會(huì)損失氣象信號，進(jìn)而影響降雨目標(biāo)參量的估計(jì)精度。Siggia等提出高斯模型自適應(yīng)處理（Gaussian model adaptive processing， GMAP）算法以降低地雜波對降雨目標(biāo)參量估計(jì)的影響［13］。然而，為了不影響未被雜波污染的降雨目標(biāo)參量的估計(jì)［14］，一般會(huì)先檢測地雜波所在位置，然后再利用此方法處理該位置的信號。目前，常用譜雜波識別（spectrum clutter identification， SCI）算法［15］和雜波抑制決策（clutter mitigation decision， CMD）算法［1617］檢測地雜波，但是這些算法會(huì)對降雨目標(biāo)造成漏判或誤判。為了提高對地雜波的識別率，Li等提出一種模糊邏輯算法，但是該方法也會(huì)造成降雨目標(biāo)的缺失［18］。為了在抑制地雜波的同時(shí)保證降雨目標(biāo)的測量精度，Warde等將地雜波的檢測和濾波集成到一種算法中，提出一種自適應(yīng)雜波環(huán)境分析法（clutter environment analysis using adaptive processing， CLEAN-AP），該方法可以通過改變凹口濾波器降低對氣象目標(biāo)的估計(jì)偏差［19］。Hubbert等提出一種時(shí)域回歸濾波方法以在濾波的同時(shí)保留較完整的降雨目標(biāo)，但是該算法還未利用實(shí)測數(shù)據(jù)進(jìn)行驗(yàn)證［20］。

近年來，譜極化濾波技術(shù)廣泛應(yīng)用于氣象雷達(dá)雜波抑制中。Unal提出利用雙譜線性去極化比（double spectral linear depolarization ratio， DsLDR）方法和以0 m/s為中心的窄凹口濾波器來抑制IDRA（IRCTR drizzlc radar）的雜波［21］，但是凹口濾波器可能會(huì)使降雨區(qū)域不連續(xù)。為了解決在抑制地雜波時(shí)導(dǎo)致的降雨區(qū)域不連續(xù)的問題，Yin等提出在濾波的同時(shí)可以利用數(shù)學(xué)形態(tài)學(xué)方法補(bǔ)償缺失的降雨［22］，但是所提方法只適用于全極化雷達(dá)。為了能解決雙極化雷達(dá)的雜波抑制問題，文獻(xiàn)［23］提出面向?qū)ο蟮淖V極化（object-orientated spectral polarimetric， OBSpol）濾波方法。該方法利用降雨目標(biāo)和雜波的譜極化特性及兩者在距離多普勒（range Doppler， RD）（平面位置顯示器（plan position indicator， PPI）每個(gè)方位向）圖上的連通特性差異，基于面向?qū)ο蟮乃枷?，篩選出所有氣象對象掩模，在盡可能多地保留降雨目標(biāo)的同時(shí)濾除非降雨目標(biāo)。但是，由于地雜波和降雨目標(biāo)的譜相關(guān)系數(shù)相似，該方法需采用窄凹口濾波器去除地雜波，故在去除地雜波的同時(shí)也會(huì)損失降雨目標(biāo)。當(dāng)雜波和降雨重疊時(shí)，移動(dòng)雙譜線去極化比（moving double spectral linear depolarization ratio， MDsLDR）和OBSpol濾波器會(huì)保留被雜波污染的降雨，從而影響對降雨的估計(jì)?？紤]到不同方位向的降雨和雜波存在差異，文獻(xiàn)［24］提出一種自適應(yīng)譜極化濾波方法，該方法可以根據(jù)不同方位向回波的譜去極化比分布情況自適應(yīng)地濾除雜波，但是該方法需要多組同一站點(diǎn)的晴空數(shù)據(jù)做支撐。

盡管上述提及的譜極化濾波器可以取得較好的雜波濾除效果，但是無法解決氣象目標(biāo)與雜波混疊情況下氣象目標(biāo)的恢復(fù)問題，不能完全滿足業(yè)務(wù)氣象雷達(dá)中雜波抑制的需要。為了解決零速度氣象信號的損失問題，Yin等人提出比GMAP算法更穩(wěn)健的Kriging算法［25］，但是該算法計(jì)算量較大，對于實(shí)時(shí)性較強(qiáng)的氣象雷達(dá)適用性不夠。Nai等人針對降雨目標(biāo)和風(fēng)車雜波混疊的情況，提出一種RD遞歸法來抑制風(fēng)車雜波［26］，但該方法需要預(yù)先知道雜波的位置信息。

針對業(yè)務(wù)極化多普勒氣象雷達(dá)中氣象目標(biāo)與雜波的混疊問題，本文提出一種基于譜連續(xù)的極化多普勒氣象雷達(dá)信號重構(gòu)方法。首先，根據(jù)氣象目標(biāo)的譜極化特征進(jìn)行譜極化濾波。其次，本文利用數(shù)學(xué)形態(tài)學(xué)方法和面向?qū)ο蟮乃枷?，在RD圖上生成一個(gè)二元掩模，并將該掩模應(yīng)用于原始譜功率可得濾波后的功率譜。然后，本文利用降雨目標(biāo)在相鄰距離庫上速度和譜寬的連續(xù)性進(jìn)行氣象窗口的補(bǔ)償。最后，根據(jù)降雨目標(biāo)的功率譜特性實(shí)現(xiàn)降雨區(qū)域的補(bǔ)償。將所提方法與OBSpol濾波器以及業(yè)務(wù)中常用的CMD+GMAP算法進(jìn)行了比較。結(jié)果表明，本文所提方法具有較好的雜波抑制與降雨保留性能。同時(shí)，本文所提方法計(jì)算時(shí)間較短，有望在業(yè)務(wù)氣象雷達(dá)中實(shí)時(shí)應(yīng)用。

本文的結(jié)構(gòu)如下：第1節(jié)主要闡述了氣象雷達(dá)的典型參量;第2節(jié)詳細(xì)介紹所提算法的設(shè)計(jì)步驟;第3節(jié)利用業(yè)務(wù)氣象雷達(dá)實(shí)測數(shù)據(jù)驗(yàn)證所提算法的性能;最后對本文工作進(jìn)行總結(jié)。

1" 氣象雷達(dá)參量

由于結(jié)合了多普勒和極化信息，譜極化技術(shù)有利于反演氣象目標(biāo)抑制雜波［2224］。假設(shè)極化氣象雷達(dá)發(fā)射x極化電磁波、接收y極化電磁波，在后向散射的約束下，譜反射功率［2223］定義為

sZxy（r，v）=C·sPxy（r，v）·r2=C·|sSxy（r，v）|2·r2（1）

式中：x，y∈{h，v}，h表示水平極化，v表示垂直極化;r為距離;v為多普勒速度;C表示雷達(dá)常數(shù);sPxy（r，v）表示譜功率;sSxy（r，v）表示二維的復(fù)RD譜。在此基礎(chǔ)上，譜差分反射率sZdr（r，v）和譜共極化相關(guān)系數(shù)sρco（r，v）［2223］可以表示為

sZdr（r，v）=10lgsZhh（r，v）sZvv（r，v）（2）

sρco（r，v）=|〈sShh（r，v）·sSvv（r，v）〉|〈|sShh（r，v）|2〉〈|sSvv（r，v）|2〉（3）

式中：sρco（r，v）中的〈·〉表示在多普勒速度維的滑動(dòng)平均;表示復(fù)共軛。當(dāng)存在地雜波時(shí)，由于地雜波和降雨的sρco相似，僅使用sρco不足以區(qū)分降雨。本文引入雜波相位排序（clutter phase alignment， CPA）［15］來輔助識別地雜波。CPA可以用來衡量后向散射相位的穩(wěn)定性，由于地物目標(biāo)后向散射相位基本是固定的，故在特定的方位和距離上，相干處理時(shí)間內(nèi)回波時(shí)間序列的相位是基本相同的［1516］。對于分布式目標(biāo)（如降雨），脈沖間的相位會(huì)隨著觀測體積內(nèi)目標(biāo)的速度和譜寬的變化而變化［1516］。故可以根據(jù)地雜波和降雨目標(biāo)的后向散射相位的穩(wěn)定性差異，利用參量CPA輔助識別地雜波。設(shè)某一距離的CPA為

CPA=∑Nn=1xn∑Nn=1|xn|（4）

式中：xn=In+j·Qn是該距離上的回波信號;N是積累脈沖數(shù)。0≤CPA≤1，地雜波的相位相對較為穩(wěn)定，CPA接近于1。噪聲的相位是隨機(jī)的，CPA接近于0。降雨目標(biāo)的CPA一般小于0.9［15］。特別地，如果降雨目標(biāo)的平均徑向速度不為0，那么CPA一般小于0.5［15］。

譜極化濾波技術(shù)可以在濾除雜波和噪聲的同時(shí)，盡可能多地保留RD圖中的降雨目標(biāo)，因此可以更準(zhǔn)確地估計(jì)降雨目標(biāo)的測量參量。雷達(dá)反射率可以表示為

Zhh（r）=C·∑v∈pre（sPhh（r，v）－sN）·r2

（5）

式中：v表示RD圖上屬于降雨目標(biāo)（precipitation， pre）的點(diǎn)，v∈pre;sN表示整個(gè)RD圖的譜噪聲［23］。本文將出現(xiàn)頻率最高的譜功率視為譜噪聲sN，譜噪聲的具體估計(jì)方法見文獻(xiàn)［27］。

利用雷達(dá)反射率可以得到徑向速度v-和譜寬σv為

v-（r）=1Zhh（r）∑v∈prev·sZhh（r，v）（6）

σv（r）=1Zhh（r）∑v∈pre（v－v（r））2·sZhh（r，v）（7）

本文將用反射率和徑向速度這兩個(gè)雷達(dá)參量評估濾波和信號恢復(fù)性能。

2" 算法設(shè)計(jì)

當(dāng)存在雜波時(shí)，降雨目標(biāo)的功率會(huì)被高估，如果該雜波是地雜波，那么降雨目標(biāo)的平均徑向速度會(huì)向零速偏移。為了提高降雨目標(biāo)參量的測量精度，本文提出一種基于譜連續(xù)的極化多普勒氣象雷達(dá)信號重構(gòu)方法，該方法可以有效地去除雜波，同時(shí)盡可能多地保留降雨目標(biāo)，其信號處理流程如圖1所示。本文所提方法主要包含濾波和信號重構(gòu)兩部分，其中濾波的目的是去除盡可能多的雜波，保留純凈的降雨目標(biāo)，但是當(dāng)雜波較強(qiáng)時(shí)，在去除雜波的同時(shí)可能會(huì)去除部分降雨目標(biāo)，這會(huì)使降雨估計(jì)出現(xiàn)偏差。本文在濾波之后加入信號重構(gòu)部分以恢復(fù)降雨目標(biāo)。

2.1" 濾" 波

最近，Yin等人提出的OBSpol濾波器在濾除雜波和噪聲上有很好的性能［23］。該濾波器分為4個(gè)步驟。第1步，利用譜極化參量進(jìn)行預(yù)濾波處理以得到RD二值圖。第2步，利用數(shù)學(xué)形態(tài)學(xué)方法恢復(fù)RD圖中的降雨信號。第3步，將RD圖中連通的區(qū)域標(biāo)記為不同的對象，并將這些對象按照面積大小進(jìn)行排序。第4步，針對每一個(gè)對象利用譜寬信息進(jìn)一步濾除每個(gè)目標(biāo)對象中的雜波和噪聲點(diǎn)，然后把這些參量合成得到濾除雜波和噪聲后的RD圖。當(dāng)作用于雙極化氣象雷達(dá)時(shí)，OBSpol濾波器選用sρco參數(shù)濾除非氣象回波。由于降雨目標(biāo)和地雜波的sρco較為相近［23］，該濾波器使用一個(gè)以零為中心的凹口濾波器輔助去除地雜波，這雖然可以濾除大部分強(qiáng)度較高的地雜波，但是也會(huì)去除降雨目標(biāo)。

文獻(xiàn)［25］指出降雨目標(biāo)的譜差分反射率標(biāo)準(zhǔn)差（standard deviation of spectral differential reflectivety， SDsZdr）集中在0～3 dB，雜波的SDsZdr集中在2～10 dB。其中，SDsZdr是對RD圖中譜反射率應(yīng)用一個(gè)3×3滑動(dòng)窗口來計(jì)算的標(biāo)準(zhǔn)差。本文利用SDsZdr來識別被地雜波污染的距離單元。濾波器的設(shè)計(jì)主要包括以下4個(gè)步驟。

步驟 1" 譜極化濾波。降雨目標(biāo)的sρco值一般遠(yuǎn)大于除地雜波以外的其他雜波，采用0.98的sρco閾值進(jìn)行譜極化濾波，其中sρco的閾值是根據(jù)降雨和雜波的去除比例確定的，具體選取方法見文獻(xiàn)［22］。為了抑制地雜波，在零頻附近對SDsZdr的RD圖進(jìn)行閾值濾波。若此區(qū)域內(nèi)SDsZdr小于等于3 dB，則認(rèn)為該區(qū)域?yàn)榻涤昴繕?biāo)，濾波掩模設(shè)為1，若此區(qū)域內(nèi)SDsZdr大于3 dB，則認(rèn)為該區(qū)域?yàn)榈仉s波，濾波掩模設(shè)為0。其中，SDsZdr的閾值也是根據(jù)降雨和雜波的去除比例確定的，具體選取方法見文獻(xiàn)［25］。該步驟充分利用了雜波和噪聲以及降雨目標(biāo)的譜極化特性差異，因此可以在去除雜波和噪聲的同時(shí)保留較多的降雨目標(biāo)。

步驟 2" 數(shù)學(xué)形態(tài)學(xué)法恢復(fù)降雨。利用sρco和SDsZdr兩參數(shù)進(jìn)行閾值濾波可能會(huì)濾除降雨目標(biāo)內(nèi)部或邊緣的弱氣象目標(biāo)。作為一種分布式目標(biāo)，降雨目標(biāo)在RD圖上是連續(xù)分布的，本文采用形態(tài)學(xué)閉運(yùn)算來恢復(fù)缺失的降雨目標(biāo)，該步驟主要填充二值圖像內(nèi)的小孔，彌合邊緣區(qū)域的小裂縫。

步驟 3" 標(biāo)記對象，生成對象掩模。一般情況下，降雨區(qū)域在RD圖上是連續(xù)分布的且面積和數(shù)量有限，本文按照面積大小進(jìn)行排序，并將其歸納為降雨對象和雜波對象。

步驟 4" 利用譜寬去除雜波。根據(jù)以上步驟可以得到濾波后的“降雨目標(biāo)區(qū)域”，但是此“降雨目標(biāo)區(qū)域”可能會(huì)包含部分雜波，根據(jù)降雨目標(biāo)的空間幾何異向性進(jìn)一步去除雜波。降雨目標(biāo)空間幾何異向性指降雨目標(biāo)沿著距離維度和多普勒維度有不同的空間幾何結(jié)構(gòu)。通過求解主要變化方向和異向比可以判斷該RD圖中是否含有雜波，從而進(jìn)一步去除雜波，具體如下。

對求得的濾波掩模在距離維求和可得掩模寬度為

wide（n）=∑Mm=1mask（m，n）（8）

式中：mask（m，n）是步驟1～步驟3求得的濾波掩模;m表示距離維;n表示多普勒維;M表示距離單元數(shù)。計(jì)算掩模寬度的整體平均寬度為

W1=1N∑Nn=1wide（n）（9）

式中：N表示多普勒單元總數(shù)。

對求得的掩模寬度按照升序排序，計(jì)算出20%到50%的平均寬度W2。若W2/W1gt;10，掩模寬度均值取W2。取掩模寬度最大值兩側(cè)的且在W1之下的極小值，認(rèn)為位于兩個(gè)極小值兩側(cè)的區(qū)域可能為剩余的雜波，濾波掩模設(shè)為0。

利用上述步驟可以得到較純凈的降雨目標(biāo)，但是利用較高閾值的譜極化濾波會(huì)損失弱降雨目標(biāo)，后面將利用選出的降雨目標(biāo)進(jìn)行信號恢復(fù)。

2.2" 信號恢復(fù)

利用較大的濾波閾值濾波會(huì)損失一些降雨目標(biāo)，從而使得在某個(gè)距離單元上沒有降雨目標(biāo)，降雨區(qū)域不連續(xù)。為了彌補(bǔ)這種不足，利用降雨目標(biāo)譜矩之間的關(guān)系以及降雨區(qū)域的分布特性進(jìn)行降雨區(qū)域的補(bǔ)償。

2.2.1" 選擇降雨區(qū)域

降雨信號通常占據(jù)較大的空間，其空間分布在距離上一般是連續(xù)的，與典型的氣象雷達(dá)分辨率相比，氣象信號在鄰近距離單元上的速度和譜寬是緩變的［25，28］，根據(jù)此假設(shè)進(jìn)行降雨目標(biāo)區(qū)域的估計(jì)。降雨目標(biāo)區(qū)域的估計(jì)主要分為3個(gè)步驟，如圖2所示。

步驟 1" 計(jì)算速度與譜寬。將第2.1節(jié)得到的濾波后RD二值圖作用于原始的功率譜，再利用式（6）和式（7）計(jì)算濾波后降雨目標(biāo)的速度和譜寬。

步驟 2" 去除速度和譜寬的異常值。利用第2.1節(jié)濾波步驟選出來的區(qū)域可能存在偏差，為了降低這些偏差對降雨區(qū)域的選擇造成影響，根據(jù)濾波之后速度和譜寬的分布情況去除速度和譜寬中的異常值，該想法來源于箱型圖去除異常值。

異常值的判斷方法如下：將求得的速度（譜寬）按照升序進(jìn)行排列，找到速度（譜寬）的上四分位數(shù)Q1與下四分位數(shù)Q3，記上邊緣Qup與下邊緣Qdown為

Qup=Q1－1.5（Q3－Q1）（10）

Qdown=Q3+1.5（Q3－Q1）（11）

將不在［Qup，Qdown］區(qū)間之內(nèi)的速度（譜寬）視為異常值。

步驟 3" 多項(xiàng)式擬合速度與譜寬。在多項(xiàng)式擬合中，擬合序列越長，擬合誤差越大。為了提高擬合精度，本文采用分段擬合的方法進(jìn)行速度和譜寬的擬合與補(bǔ)償。將序列分成4個(gè)子序列，利用不同的加權(quán)以平滑過渡區(qū)域［29］。子序列的分段方式如圖3所示。

在本次分段過程中，每一段都和上一段有所重疊，重疊的部分利用不同的加權(quán)系數(shù)進(jìn)行組合以提高每一段邊緣點(diǎn)擬合的準(zhǔn)確度。具體加權(quán)系數(shù)按照重疊的點(diǎn)數(shù)確定。例如，第1段和第2段重疊的點(diǎn)數(shù)為L，此時(shí)重疊部分的計(jì)算方式為

val（l）=L+1－lL+1a（l）+lL+1b（l）（12）

式中：l為該段中需要擬合的距離庫;a（l）為第1段擬合得到的重疊部分的結(jié)果;b（l）為第2段擬合得到的重疊部分的結(jié)果。這里考慮到擬合的精度和擬合的效率，將重疊的點(diǎn)數(shù)設(shè)置為子序列點(diǎn)數(shù)的1/5。

步驟 4" 選擇降雨區(qū)域。利用多項(xiàng)式擬合得到每個(gè)距離單元的速度和譜寬后，利用降雨目標(biāo)速度與譜寬之間的關(guān)系進(jìn)行降雨區(qū)域的補(bǔ)償，距離r處的降雨區(qū)域［26］為［v-（r）－K·σv（r），v-（r）+K·σv（r）］。其中，K是降雨區(qū)域?qū)挾纫蜃?。?dāng)信雜比（signal-to-clutter ratio， SCR）較高時(shí)，由于地雜波及其旁瓣信號對氣象信號的影響較小，K要取相對較大的值，以便將氣象信號全部包含進(jìn)來;當(dāng)SCR較低時(shí)，雜波作為主要的成分，每個(gè)譜系數(shù)都有顯著的偏差，為了降低雜波對降雨信號譜矩估計(jì)的影響，K值會(huì)取得小一些［26］。本文將K取為2.4［26］，其中K的選取是根據(jù)雜信比的大小利用RD遞歸算法計(jì)算得到，其具體估計(jì)方法見文獻(xiàn)［25］。

2.2.2" 補(bǔ)償降雨目標(biāo)

降雨目標(biāo)回波的功率譜密度通常被認(rèn)為近似于高斯分布，即

S（v）=Pr2πσvexp－（v－v-）22σ2v+cc（13）

式中：Pr是該距離單元降雨目標(biāo)的峰值功率;v-為降雨目標(biāo)的平均速度;σv為降雨目標(biāo)的譜寬;cc是噪聲功率。

相對于降雨目標(biāo)，地物雜波的譜寬較窄，且其平均多普勒速度基本為零，但是降雨目標(biāo)的平均多普勒速度可以是奈奎斯特速度中的任意值。一般情況下，降雨目標(biāo)的譜寬比地雜波的譜寬大，故可以利用此特征重構(gòu)與地雜波混疊的降雨目標(biāo)。降雨目標(biāo)功率譜補(bǔ)償?shù)乃惴ㄈ鐖D4所示，當(dāng)降雨區(qū)域與地雜波重疊時(shí)，先用凹口濾波器去除該距離庫的地雜波，再利用高斯擬合重構(gòu)與地雜波重合的降雨。

由于地雜波和降雨的sρco相似，本文引入CPA來輔助識別地雜波。CPA的概率分布與速度的關(guān)系如圖5所示?？梢钥闯觯蟛糠终瓗Я闼倌繕?biāo)的CPA值大于0.88，本文將CPA的閾值設(shè)為0.88，可以識別大部分地雜波。

在去除地雜波之后，可以通過定義均方根誤差（root mean square error， RMSE）來確定該距離單元上的峰值功率。對于某個(gè)距離上降雨目標(biāo)的功率譜，如果我們有H個(gè)多普勒單元包含降雨，那么該距離上功率的RMSE可以表征為

δS（Pr）=1M∑vHv1（Stru（v）－Sest（Pr，v））（14）

式中：Stru（v）為該距離上速度為v處的測量真值;Sest（Pr，v）為當(dāng)峰值功率取Pr時(shí)速度為v處的估計(jì)值;Pr可?。跴rmin，Prmax］內(nèi)任意值;Prmin為噪聲功率;Prmax為此方位向該距離處的峰值功率。為了更好地恢復(fù)與地雜波混疊的降雨目標(biāo)，取使得δS（Pr）最小的峰值功率Pr為高斯擬合的峰值功率。

3" 實(shí)測數(shù)據(jù)驗(yàn)證

3.1" 雷達(dá)系統(tǒng)與雜波特性分析

將所提算法與OBSpol濾波方法以及業(yè)務(wù)中常用的CMD+GMAP算法進(jìn)行對比，并利用廣州站S波段氣象雷達(dá)的實(shí)測數(shù)據(jù)驗(yàn)證所提算法的性能。該雷達(dá)工作于同時(shí)發(fā)射同時(shí)接收（simultaneous transmission and simultaneous reception， STSR）模式，雷達(dá)的具體參數(shù)如表1所示，其中掃描模式為體積掃描（volume coverage pattern， VCP）21模式。

為了驗(yàn)證所提算法的雜波抑制與降雨保留性能，使用該雷達(dá)在0.5°仰角上測量得到的含降雨目標(biāo)和雜波的實(shí)測數(shù)據(jù)進(jìn)行性能驗(yàn)證，該數(shù)據(jù)測量于2016年5月20日，其原始雷達(dá)PPI如圖6所示。

由圖6可知，0.5°仰角時(shí)地雜波對降雨的影響較大。圖6中第56個(gè)徑向的譜極化參量如圖7所示。

由圖7（a）可見，該方位向近地區(qū)除了含有降雨目標(biāo)外還有地雜波及地雜波的旁瓣。地雜波的能量主要集中在零多普勒處，但是地雜波的旁瓣主要有以下幾個(gè)特點(diǎn)：① 在其出現(xiàn)的距離庫遍布整個(gè)奈奎斯特速度區(qū)間;② 在不同方位向地雜波及其旁瓣出現(xiàn)的距離庫不同;③ 地雜波旁瓣與降雨目標(biāo)的強(qiáng)度相當(dāng)，且有可能會(huì)貫穿降雨目標(biāo)。地雜波及其旁瓣的特征增加了雜波抑制的難度。

信噪比和SCR對雷達(dá)雙偏振參量的影響較大［30］。由圖7（b）可見，降雨目標(biāo)和地雜波旁瓣的sρco差別不大。這意味著僅用sρco參數(shù)無法提取較完整的降雨目標(biāo)，故在濾波之后需要利用降雨目標(biāo)的空間連續(xù)性進(jìn)行信號恢復(fù)。

3.2" 譜極化處理

為了評估本文所提算法的濾波及降雨保留性能，利用不同的處理方法處理圖7中S波段氣象雷達(dá)的數(shù)據(jù)，結(jié)果如圖8所示。

圖8（a）是OBSpol濾波器濾波后得到的結(jié)果，圖8（b）是采用CMD+GMAP方法處理后得到的結(jié)果，圖8（c）是使用本文所提算法處理得到的結(jié)果。對比圖8（a）和圖8（c），兩種濾波方法都可以有效地區(qū)別地雜波與降雨目標(biāo)，但是由于地雜波的sρco和降雨目標(biāo)相當(dāng)，OBSpol濾波器在選擇降雨區(qū)域的同時(shí)也會(huì)保留一部分地雜波的旁瓣和其他雜波（見圖8（a）中矩形框）。OBSpol濾波器在零多普勒速度附近使用一個(gè)凹口來去除地雜波，故當(dāng)降雨目標(biāo)和地雜波混疊時(shí)，去除地雜波的同時(shí)也會(huì)損失降雨目標(biāo)（見圖8（a）中圓圈所示）。80 km以內(nèi)的地雜波及其旁瓣較強(qiáng)且會(huì)覆蓋降雨，OBSpol濾波器濾波時(shí)會(huì)將與地雜波旁瓣重疊的降雨目標(biāo)去除（見圖8（a）中圓圈所示）。本文所提的算法處理效果較好，且在50～70 km處能夠基本恢復(fù)與地雜波旁瓣混疊的降雨目標(biāo)。

對比圖8（b）和圖8（c），業(yè)務(wù)中常用的CMD+GMAP算法對地雜波有一定的抑制效果，在地雜波主瓣和降雨混疊的情況下，該算法能在抑制地雜波的同時(shí)恢復(fù)降雨，但是當(dāng)?shù)仉s波旁瓣較高時(shí)，該方法會(huì)保留一部分地雜波旁瓣，這也會(huì)使得降雨估計(jì)出現(xiàn)偏差。相比而言，本文所提方法不僅能將與地雜波混疊的降雨信號恢復(fù)，還能去除地雜波旁瓣，具有較好的雜波抑制效果。

為了定量衡量不同方法在不同SCR情況下的有效性，在RD圖上定義檢測概率Pd和虛警概率Pfa為

Pd=TPTP+FN

Pfa=FPFP+TN（15）

式中：TP表示降雨被判斷為降雨的點(diǎn)數(shù);FN表示降雨被判斷為雜波的點(diǎn)數(shù);FP表示雜波被判斷為降雨的點(diǎn)數(shù);TN表示雜波被判斷為雜波的點(diǎn)數(shù)［9］。第56個(gè)方位向不同SCR情況下的Pd和Pfa結(jié)果如圖9所示，考慮到業(yè)務(wù)中常用的CMD+GMAP算法沒有選擇降雨并去除噪聲的能力，故這里只對OBSpol濾波器和所提方法進(jìn)行了對比。

由圖9可知，所提算法的Pd在任何SCR情況下均高于OBSpol方法，且Pd趨近于1，即所提方法可以將降雨完整地保留下來，但是OBSpol方法在SCR較低的情況下，Pd會(huì)降低至0.2以下。由Pfa可知，相對于OBSpol方法，所提方法的Pfa低于0.1。

從實(shí)時(shí)性上來說，將所有的算法在Windows 10的電腦上運(yùn)行，其中電腦的處理器為i7-1165G7，2.8 GHz的主頻，16 GB的內(nèi)存，OBSpol濾波器用時(shí)0.2 s，CMD+GMAP算法用時(shí)46.4 s，所提算法用時(shí)0.7 s。故綜合考慮處理效果和處理時(shí)間，所提算法有較好的雜波抑制與降雨保留性能，同時(shí)耗時(shí)較短，實(shí)時(shí)性相對較強(qiáng)。

3.3" PPI處理與結(jié)果分析

將所提算法應(yīng)用于全掃描角度下的每個(gè)RD譜，其雷達(dá)反射率以及徑向速度的PPI顯示結(jié)果如圖10所示。OBSpol濾波器處理后得到的雷達(dá)反射率和徑向速度如圖10（a）和圖10（c）所示，圖10（b）和圖10（d）是利用文中所提算法得到的雷達(dá)反射率和徑向速度。對比圖10（a）和圖10（b），OBSpol濾波器在去除地雜波的同時(shí)會(huì)損失大量降雨目標(biāo)。相比而言，本文所提算法得到的降雨區(qū)域較完整。對比圖10（c）和圖10（d），OBSpol濾波得到的徑向速度不連續(xù)且存在大量野值點(diǎn)，本文得到的徑向速度野值點(diǎn)較少且速度較為連續(xù)。

利用該站點(diǎn)1.5°仰角的實(shí)測數(shù)據(jù)驗(yàn)證所提算法的性能。圖11為經(jīng)過OBSpol濾波以及所提方法處理后的雷達(dá)反射率和徑向速度。對比圖11（b）與圖11（c），在保留降雨方面，本文所提算法性能優(yōu)于OBSpol濾波器。對比圖11（e）與圖11（f），相比于OBSpol濾波器，本文所提算法得到的徑向速度更加連續(xù)且野值點(diǎn)較少。

4" 結(jié)" 論

針對氣象雷達(dá)中氣象目標(biāo)和雜波混疊情況下的雜波抑制問題，本文提出一種基于譜連續(xù)的極化氣象雷達(dá)信號恢復(fù)方法。氣象雷達(dá)實(shí)測數(shù)據(jù)的驗(yàn)證表明，本文所提算法在雜波抑制和降雨目標(biāo)保留方面，相較于OBSpol濾波器以及業(yè)務(wù)中常用的CMD+GMAP，具有更好的性能，且算法相對簡單，有望應(yīng)用于業(yè)務(wù)氣象雷達(dá)。

參考文獻(xiàn)

［1］ DOVIAK R J， ZRNIC D S. Doppler radar and weather observations［M］. Mineola： Dover Publications， 2006.

［2］ BRINGI V N， CHANDRASEKAR V. Polarimetric Doppler weather radar： principles and applications［M］. Cambridge： Cambridge University Press， 2001.

［3］ ZHANG G F. Weather radar polarimetry［M］. Boca Raton： CRC Press， 2016.

［4］ MATROSOV S Y， CLARK K A， MARTNER B E， et al. X-band polarimetric radar measurements of rainfall［J］. Journal of Applied Meteorology， 2002， 41（9）： 941952.

［5］ RYZHKOV A V， SCHUUR T J， BURGESS D W， et al. The joint polarization experiment： polarimetric rainfall measurements and hydrometeor classification［J］. Bulletin of the American Meteorological Society， 2005， 86（6）： 809824.

［6］ CHEN H， CHANDRASEKAR V. The quantitative precipitation estimation system for dallas-fort worth （DFW） urban remote sensing network［J］. Journal of Hydrology， 2015， 531（2）： 259271.

［7］ DIXON M， WIENER G. TITAN： thunderstorm identification， tracking， analysis， and nowcasting—a radar-based methodology［J］. Journal of Atmospheric and Oceanic Technology， 1993， 10（6）： 785797.

［8］ FUKAO S， HAMAZU K， DOVIAK R J. Radar for meteorological and atmospheric observations［M］. Tokyo： Springer， 2014.

［9］ 殷加鵬， 李健兵， 龐晨， 等. 一種極化多普勒氣象雷達(dá)的射頻干擾濾波方法［J］. 雷達(dá)學(xué)報(bào)， 2021， 10（6）： 905918.

YIN J P， LI J B， PANG C， et al. A radio frequency interference mitigation method for polarimetric Doppler weather radar［J］. Journal of Radars， 2021， 10（6）： 905918.

［10］ YIN J P， HOOGEBOOM P， UNAL C， et al. Radio frequency interference characterization and mitigation for polarimetric weather radar： a study case［J］. IEEE Trans.on Geoscience and Remote Sensing， 2021， 60： 5102616.

［11］ MEISCHNER P. Weather radar： principles and advance applications［M］. Berlin： Springer， 2004.

［12］ GROGINSKY H L， GLOVER K M. Weather radar canceller design［C］∥Proc.of the 19th Conference on Radar Meteorology， 1980： 192198.

［13］ SIGGIA A D， PASSARELLI R. Gaussian model adaptive processing （GMAP） for improved ground clutter cancellation and moment calculation［C］∥Proc.of the 3rd European Conference on Radar Meteorology， 2004： 6773.

［14］ WU C， LIU L P， CHEN C， et al. Challenges of the polarimetric update on operational radars in China—ground clutter contamination of weather radar observations［J］. Remote Sensing， 2021， 13（2）： 217235.

［15］ LI Y G， ZHANG G F， DOVIAK R J， et al. A new approach to detect ground clutter mixed with weather signals［J］. IEEE Trans.on Geoscience and Remote Sensing， 2013， 51（4）： 23732387.

［16］ HUBBERT J， DIXON M， ELLIS S， et al. Weather radar ground clutter. Part I： identification， modeling， and simulation［J］. Journal of Atmospheric and Oceanic Technology， 2009， 26（7）： 11651180.

［17］ HUBBERT J， DIXON M， ELLIS S， et al. Weather radar ground clutter. Part II： real-time identification and filtering［J］. Journal of Atmospheric and Oceanic Technology， 2009， 26（7）： 11811197.

［18］ LI N， WANG Z H， SUN K Y， et al. A quality control method of ground-based weather radar data based on statistics［J］. IEEE Trans.on Geoscience and Remote Sensing， 2017， 56（4）： 22112219.

［19］ WARDE D A， TORRES S M. The autocorrelation spectral density for Doppler-weather-radar signal analysis［J］. IEEE Trans.on Geoscience and Remote Sensing， 2014， 52（1）： 508518.

［20］ HUBBERT J， MEYMARIS G， ROMATSCHKE U， et al. Improving signal statistics using a regression ground clutter filter. Part 1： theory and simulations［J］. Journal of Atmospheric and Oceanic Technology， 2021， 38（8）： 13531375.

［21］ UNAL C. Spectral polarimetric radar clutter suppression to enhance atmospheric echoes［J］. Journal of Atmospheric and Oceanic Technology， 2009， 26（9）： 17811797.

［22］ YIN J P， UNAL C M H， RUSSCHENBERG H W J. Narrow-band clutter mitigation in spectral polarimetric weather radar［J］. IEEE Trans.on Geoscience and Remote Sensing， 2017， 55（8）： 46554667.

［23］ YIN J P， UNAL C， HOOGEBOOM P. Object-orientated filter design in spectral domain for polarimetric weather radar［J］. IEEE Trans.on Geoscience and Remote Sensing， 2018， 57（5）： 27252740.

［24］ 安孟昀， 殷加鵬， 黃建開， 等. 一種雙極化氣象雷達(dá)自適應(yīng)譜極化濾波方法［J］. 雷達(dá)學(xué)報(bào)， 2022， 11（3）： 408417.

AN M Y， YIN J P， HUANG J K， et al. Adaptive spectral polarization filter design for dual-polarization weather radar［J］. Journal of Radars， 2022， 11（3）： 408417.

［25］ YIN J P， SCHLEISS M， WANG X S. Clutter-contaminated signal recovery in spectral domain for polarimetric weather radar［J］. IEEE Geoscience and Remote Sensing Letters， 2021， 19： 3504505.

［26］ NAI F， TORRES S， PALMER R. On the mitigation of wind turbine clutter for weather radars using range-Doppler spectral processing［J］. IET Radar， Sonar and Navigation， 2013， 7（2）： 178190.

［27］ IVIC I R， CURTIS C， TORRES S M. Radial-based noise po-wer estimation for weather radars［J］. Journal of Atmospheric and Oceanic Technology， 2021， 38（6）： 12311232.

［28］ MORSE C S， GOODRICH R K， CORNMAN L B. The NIMA method for improved moment estimation from Doppler spectra［J］. Journal of Atmospheric and Oceanic Technology， 2002， 19（3）： 274295.

［29］ HUBBERT J， MEYMARIS G， ROMATSCHKE U， et al. Using a regression ground clutter filter to improve weather radar signal statistics： theory and simulations［J］. Journal of Atmospheric and Oceanic Technology， 2021， 38（8）： 13531375.

［30］ 張帥. 雙偏振天氣雷達(dá)雜波識別方法研究［D］. 南京： 南京信息工程大學(xué)， 2022.

ZHANG S. Study on clutter identification method based on dual-polarization weather radar［D］. Nanjing： Nanjing University of Information Science and Technology， 2022.

作者簡介

安孟昀（1997—），女，博士研究生，主要研究方向?yàn)闃O化雷達(dá)信號處理。

殷加鵬（1990—），男，副研究員，博士，主要研究方向?yàn)闃O化雷達(dá)信號處理。

黃建開（1994—），男，博士研究生，主要研究方向?yàn)闃O化雷達(dá)信號處理。

李永禎（1977—），男，研究員，博士研究生導(dǎo)師，博士，主要研究方向?yàn)槔走_(dá)極化信息處理、空間電子對抗、目標(biāo)檢測與識別。

王雪松（1972—），男，教授，博士研究生導(dǎo)師，博士，主要研究方向?yàn)闃O化雷達(dá)信號處理、雷達(dá)目標(biāo)識別、雷達(dá)電子對抗。