北京MST雷達(dá)功率譜密度處理算法改進(jìn)

陳 澤 田玉芳* 呂達(dá)仁

1)(中國科學(xué)院大氣物理研究所中層大氣和全球環(huán)境探測重點(diǎn)實(shí)驗(yàn)室，北京 100029) 2)(中國科學(xué)院大氣物理研究所香河大氣綜合觀測試驗(yàn)站，香河 065400) 3)(中國科學(xué)院大學(xué)，北京 100049)

引 言

MST(mesosphere-stratosphere-troposphere，中間層-平流層-對流層)雷達(dá)根據(jù)大氣折射指數(shù)不規(guī)則體對雷達(dá)發(fā)射電磁波產(chǎn)生的后向散射，獲取對流層、下平流層以及中間層的大氣三維風(fēng)場及大氣折射率結(jié)構(gòu)常數(shù)等信息。可對3～25 km，60～90 km高度及以上大氣層進(jìn)行高時(shí)空分辨率探測，是研究大氣多層區(qū)動力學(xué)特征與過程的獨(dú)特手段。MST雷達(dá)的回波機(jī)制主要有3種：湍流散射、鏡面反射以及熱散射。大氣湍流活動會造成大氣折射率分布隨機(jī)起伏，從而對電磁波產(chǎn)生散射，稱為湍流散射，主要發(fā)生在對流層及以下大氣層。大氣中折射率水平分布較均勻，但在垂直方向上存在很大梯度，使用米波段或分米波段雷達(dá)探測該結(jié)構(gòu)時(shí)，存在鏡面反射機(jī)制，如對流層頂區(qū)。高層大氣中離子、電子熱運(yùn)動形成的自由電子密度隨機(jī)起伏，造成大氣折射率不均勻分布，雷達(dá)電磁波產(chǎn)生的散射稱為熱散射，在中間層-低熱層起作用。大氣湍流散射是MST雷達(dá)回波的主要成因，有效散射體的空間尺度與入射電磁波半波長相當(dāng)。MST雷達(dá)回波非常微弱，同時(shí)MST雷達(dá)目標(biāo)信號也會受到本機(jī)噪聲、地雜波以及間歇性雜波等影響。因此，從有干擾的弱信號中提取大氣信號，需經(jīng)過信號處理與數(shù)據(jù)處理。

國內(nèi)外有關(guān)提高獲取目標(biāo)信號能力和處理信號的譜分析方法已取得一些進(jìn)展[1-16]。目前北京MST雷達(dá)存儲的原始數(shù)據(jù)是經(jīng)信號處理后的雷達(dá)功率譜密度，在此基礎(chǔ)上改進(jìn)數(shù)據(jù)處理算法既可提高實(shí)時(shí)觀測數(shù)據(jù)質(zhì)量，也可提高歷史數(shù)據(jù)質(zhì)量。

功率譜密度中存在的噪聲與其他雜波，對目標(biāo)信號(晴空大氣湍流信號)提取有很大影響。數(shù)據(jù)處理算法的理論基礎(chǔ)可概括為大氣風(fēng)場在空間與時(shí)間上具有一定連續(xù)性和不同類型的信號在功率譜密度圖上性質(zhì)不同，如最大鏈檢測法[17]、高斯擬合判斷[14]、NIMA方法[18]與綜合識別法等[19-26]。此外，還有研究通過抑制降水[27-28]、地雜波[29]等干擾信號以提高水平風(fēng)場數(shù)據(jù)質(zhì)量[30]。由于技術(shù)限制和算法本身的復(fù)雜性，有關(guān)數(shù)據(jù)處理的研究有限。值得注意的是，數(shù)據(jù)處理結(jié)果對噪聲電平估算[31-34]準(zhǔn)確性十分敏感，噪聲電平估算是雷達(dá)數(shù)據(jù)處理過程中的一個(gè)非常關(guān)鍵步驟。

當(dāng)前北京MST雷達(dá)數(shù)據(jù)除了風(fēng)場數(shù)據(jù)，其他數(shù)據(jù)產(chǎn)品未得到充分應(yīng)用，主要因?yàn)槔走_(dá)數(shù)據(jù)處理系統(tǒng)所使用的算法仍需要改進(jìn)。如噪聲電平估算明顯偏小，未能對目標(biāo)回波進(jìn)行精確識別，雜波干擾未被有效抑制，雖然對一階矩(多普勒速度)影響較小，但會導(dǎo)致零階矩(回波功率)與二階中心矩(譜寬)的值偏大。因此，本文將從噪聲電平估算與目標(biāo)回波準(zhǔn)確識別與提取兩個(gè)方面改進(jìn)，在噪聲電平估算過程采用對數(shù)-線性擬合方案實(shí)現(xiàn)客觀分析法，以提高實(shí)時(shí)觀測與歷史觀測數(shù)據(jù)的質(zhì)量。

1 數(shù)據(jù)及北京MST雷達(dá)簡介

1.1 數(shù) 據(jù)

本文使用北京MST雷達(dá)數(shù)據(jù)、北京探空數(shù)據(jù)和ERA5再分析數(shù)據(jù)。

北京MST雷達(dá)數(shù)據(jù)：選用2012年1月1日00:00—12月31日23:59(世界時(shí)，下同)的中模式數(shù)據(jù)，每30 min探測1次，垂直分辨率為600 m。

北京探空數(shù)據(jù)：選用2012年1月1日—12月31日北京市觀象臺探空的水平風(fēng)場數(shù)據(jù)，每日兩次探測，探空球施放時(shí)間為11:15與23:15，每1～2 s記錄1次，垂直分辨率約為10 m。為方便計(jì)算，將探空數(shù)據(jù)插值成垂直分辨率為50 m的數(shù)據(jù)集。

ERA5再分析數(shù)據(jù)：選用2012年1月1日00：00—12月31日23：00單點(diǎn)(39.75°N，117°E)模式高度層的垂直速度數(shù)據(jù)，時(shí)間分辨率為1 h。將垂直速度線性插值為垂直分辨率為600 m的數(shù)據(jù)集，與北京MST雷達(dá)數(shù)據(jù)高度相同。

1.2 北京MST雷達(dá)

2 MST雷達(dá)功率譜密度處理算法改進(jìn)

MST雷達(dá)信號并非只有目標(biāo)信號，還存在非目標(biāo)信號，因此對功率譜密度進(jìn)行處理，精確提取目標(biāo)信號是一項(xiàng)必需的工作。MST雷達(dá)信號可分為兩大類[37]，即大氣要素與非大氣要素信號。大氣要素信號包括晴空大氣湍流、降水等，非大氣要素則包括地雜波(如地面的建筑、山、樹等造成的信號)、間歇性雜波(如飛機(jī)、鳥類的飛行等所造成的信號)和系統(tǒng)噪聲等。當(dāng)以大氣湍流為探測目標(biāo)，則噪聲、間歇性雜波等非目標(biāo)信號常對基礎(chǔ)數(shù)據(jù)估算有很強(qiáng)干擾。

2.1 噪聲電平估算

在功率譜密度圖中，一般認(rèn)為高于噪聲電平的每個(gè)包絡(luò)均對應(yīng)著信號。因此噪聲電平估算的準(zhǔn)確性不僅對SNR估算有重大影響，同時(shí)在信號提取中也有重要作用。

為抑制直流與窄而強(qiáng)的回波信號的影響，在處理功率譜密度前需進(jìn)行預(yù)處理，包括零頻去直流以及三點(diǎn)平滑處理(如圖1所示)。

圖1 2012年1月4日00:10北京MST雷達(dá)原始功率譜密度及數(shù)據(jù)預(yù)處理疊加經(jīng)客觀分析法和八段平均法處理的噪聲電平(15 km高度東波束)

估算噪聲電平時(shí)，最常用的方法是Hildebrand等[31]提出的客觀分析法?？陀^分析法的基本原理是假定雷達(dá)測量的噪聲信號屬于高斯白噪聲，將功率譜密度在某個(gè)待定噪聲閾值以下的點(diǎn)組成該距離庫的噪聲功率譜密度序列記為{Sm}，頻率范圍為F，第m點(diǎn)對應(yīng)的頻率為fm(m=1,2,…,M)，通過調(diào)整預(yù)置的噪聲閾值，當(dāng)噪聲功率譜密度序列計(jì)算出的方差與高斯白噪聲理論給出的方差相等時(shí)，該預(yù)設(shè)噪聲閾值即為所求噪聲閾值，該序列{Sm}的平均值即為噪聲電平。Hildebrand等[31]根據(jù)高斯白噪聲的均勻分布與高斯分布兩個(gè)特征，給出如下兩個(gè)方程組:

(1)

(2)

其中，fm為第m點(diǎn)對應(yīng)的頻率，共有M個(gè)點(diǎn)，{Sm}為噪聲功率譜密度序列，F(xiàn)為頻率范圍，p是數(shù)據(jù)預(yù)處理時(shí)滑動平均的點(diǎn)數(shù)，如果不進(jìn)行滑動平均處理，則p=1。

通過調(diào)整預(yù)設(shè)噪聲閾值，當(dāng)R1=R2≈1時(shí)，即得到所需噪聲閾值。實(shí)際應(yīng)用中，常規(guī)方案的噪聲閾值預(yù)設(shè)從該距離庫中功率譜密度最大值Smax開始向下搜索，將Smax與預(yù)設(shè)噪聲閾值之差(dP)作為自變量計(jì)算R1和R2。R1從大于1的值開始以不規(guī)則方式逐漸接近1，而R2從小于1的值開始逐漸增大到1以上(詳見文獻(xiàn)[31]中圖2b)，可見R2=1可作為確定噪聲閾值的循環(huán)終止條件。以上客觀分析法的基本原理與實(shí)現(xiàn)方法可參考Hildebrand等[31]和胡明寶[20]的相關(guān)工作。

由于客觀分析法的實(shí)現(xiàn)存在一定難度，實(shí)際應(yīng)用時(shí)多采用分段平均法。該方法是基于噪聲服從自由度為2N/k的χ2分布的假設(shè)，將FFT(fast Fourier transform)點(diǎn)數(shù)為N的頻譜均分為k段，每一段數(shù)據(jù)計(jì)算其均值，將最小的平均值視為整個(gè)頻譜的噪聲電平值。根據(jù)Petitdidier等[32]的研究，當(dāng)k=8或k=16時(shí)效果最佳，實(shí)際使用一般取k=8，稱為八段平均法。該方法算法復(fù)雜度低，易實(shí)現(xiàn)，且較以低于峰值10 dB或15 dB 作為噪聲閾值估算噪聲電平值的方法更客觀[33]。

經(jīng)過比對客觀分析法與八段平均法的估算結(jié)果，發(fā)現(xiàn)八段平均法計(jì)算的噪聲電平較低(圖1)。對八段平均法，在功率譜密度圖上，高于噪聲電平的某些包絡(luò)會出現(xiàn)多個(gè)極大值，其中部分代表噪聲信號的包絡(luò)譜寬值異常偏大，不利于后續(xù)處理，而客觀分析法在此方面明顯優(yōu)于八段平均法。

客觀分析法通過從功率譜密度最大值不斷向下調(diào)整噪聲閾值估算噪聲電平，算法實(shí)現(xiàn)耗時(shí)長，選取合適的調(diào)整步長比較困難。為了快速實(shí)現(xiàn)客觀分析法，本文使用兩種方案，即對數(shù)-線性擬合方案與二分法方案。

二分法方案基本思路：對于某距離庫的功率譜密度，當(dāng)預(yù)設(shè)的噪聲閾值取功率譜密度最大值Smax時(shí)計(jì)算的R2<1，即目標(biāo)信號未淹沒于噪聲；接下來若預(yù)設(shè)的噪聲閾值取八段平均法計(jì)算出的噪聲電平Nk=8，通常計(jì)算出的R2>1，即實(shí)際噪聲閾值在Smax與Nk=8之間。因此可將Nk=8和Smax作為預(yù)設(shè)噪聲閾值的初始上下邊界值，通過二分法(也稱折半法)縮小噪聲閾值的搜索區(qū)域，達(dá)到快速實(shí)現(xiàn)客觀分析法的目的，該方案與客觀分析法常規(guī)方案的原理完全一致，是一種常規(guī)方案。

在實(shí)際計(jì)算過程中，對大量距離庫功率譜密度處理時(shí)發(fā)現(xiàn)，R2隨dP增長呈類似指數(shù)變化趨勢，若將lgR2與dP做線性擬合，則當(dāng)lgR2=0時(shí)，對應(yīng)(Smax-dPlgR2=0)為所求噪聲閾值(圖2)。由于圖2b中l(wèi)gR2=0時(shí)對應(yīng)的dP與圖2a中R2=1所對應(yīng)的dP值相同，因此可以考慮應(yīng)用先取對數(shù)再進(jìn)行線性擬合的方案實(shí)現(xiàn)客觀分析法，確定噪聲閾值并估算噪聲電平，具體步驟：①dP選取計(jì)算序列，計(jì)算對應(yīng)的R2值；②對lgR2與dP進(jìn)行線性擬合，當(dāng)lgR2=0時(shí)所對應(yīng)的dP為目標(biāo)值，(Smax-dPlgR2=0)為噪聲閾值；③根據(jù)所得噪聲閾值，選取噪聲序列，計(jì)算此序列的平均值即為噪聲電平。

圖2 2012年3月16日15:10北京MST雷達(dá)觀測數(shù)據(jù)的R2及l(fā)gR2隨dP變化

為進(jìn)一步檢驗(yàn)對數(shù)-線性擬合方案的有效性和可靠性，利用北京MST雷達(dá)2012年5月7.8～12 km 高度共1715個(gè)距離庫的功率譜密度，分別采用對數(shù)-線性擬合方案及二分法方案估算噪聲電平值并進(jìn)行比較(如圖3所示)。結(jié)果表明：對數(shù)-線性擬合方案估算的噪聲電平與二分法方案估算值有很好一致性，二者差值的標(biāo)準(zhǔn)差、平均絕對差值與相關(guān)系數(shù)分別是0.43 dB，0.21 dB與0.998，前者的平均值為168.6 dB，后者的平均值為168.5 dB，表明對數(shù)-線性擬合方案可行有效。

圖3 2012年5月北京MST雷達(dá)7.8～12 km高度內(nèi)(1715個(gè)距離庫)功率譜數(shù)據(jù)應(yīng)用對數(shù)-線性擬合方案與二分法方案估算噪聲電平對比

相比二分法方案，對數(shù)-線性擬合方案具有明顯優(yōu)勢：①能夠快速且準(zhǔn)確估算噪聲電平值。對數(shù)-線性擬合方案不需要預(yù)設(shè)參數(shù)，而二分法方案或Hildebrand等[31]方案，需要提前設(shè)定參數(shù)，且計(jì)算結(jié)果與計(jì)算速度對參數(shù)設(shè)定十分敏感；②對數(shù)-線性擬合方案可判斷各距離庫能否提取出噪聲以外的目標(biāo)回波信號，即當(dāng)擬合直線的斜截距，即預(yù)設(shè)噪聲閾值選取功率譜密度最大值時(shí)對應(yīng)的lgR2>0時(shí)，目標(biāo)信號淹沒于噪聲中，該距離庫將不進(jìn)行下一步的目標(biāo)識別處理。而分段平均法則無法進(jìn)行該判斷。因此本文提出的對數(shù)-線性擬合方案快速實(shí)現(xiàn)客觀分析法估算噪聲電平是雷達(dá)功率譜數(shù)據(jù)處理過程中的重要改進(jìn)，具有實(shí)際應(yīng)用價(jià)值。

2.2 目標(biāo)回波識別與提取

功率譜密度進(jìn)行預(yù)處理后，需要進(jìn)一步識別并提取目標(biāo)回波信號。根據(jù)目標(biāo)信號在功率譜密度圖中的特點(diǎn)，將峰值最大、譜寬最大與功率最大的峰區(qū)稱為三峰[22-23]，目標(biāo)信號通常存在于三峰，視此三峰為對應(yīng)距離庫上的準(zhǔn)目標(biāo)信號。本文實(shí)現(xiàn)目標(biāo)信號的識別與提取的步驟：①提取所有回波；②在所有信號中確定三峰為準(zhǔn)目標(biāo)信號；③以三峰為主鏈進(jìn)行鏈?zhǔn)綑z驗(yàn)，對整個(gè)距離庫進(jìn)行目標(biāo)信號的識別；④對稱性檢驗(yàn)，進(jìn)一步抑制干擾信號；⑤根據(jù)地雜波對應(yīng)的峰區(qū)特點(diǎn)對其進(jìn)行抑制，最終準(zhǔn)確識別和提取出目標(biāo)回波。

2.2.1 信號提取

為提取目標(biāo)信號，應(yīng)當(dāng)識別功率譜中所有信號。以噪聲電平為閾值，將功率譜中大于閾值的每個(gè)峰區(qū)視為有效信號(Si)。首先識別每個(gè)峰區(qū)的極大值(Ci),再識別每個(gè)極大值左右兩側(cè)的極小值或噪聲電平作為峰區(qū)的左右邊界(li,ri)，則Si=Si(li,Ci,ri)。

2.2.2 鏈?zhǔn)綑z驗(yàn)

Clothiaux等[17]以鏈最長以及功率譜值之和最大為原則，提出最大鏈檢測法進(jìn)行目標(biāo)回波信號的識別與提取。基于大氣水平風(fēng)場具有時(shí)空連續(xù)性，且雷達(dá)以大氣折射指數(shù)不規(guī)則體，即大氣湍渦為探測目標(biāo)，作為大氣風(fēng)場探測的示蹤物，使雷達(dá)能夠探測到這種時(shí)空連續(xù)性的特點(diǎn)，結(jié)合最大鏈檢測法以及三峰的概念，本文應(yīng)用鏈?zhǔn)綑z驗(yàn)的方法進(jìn)行目標(biāo)回波識別。

某距離庫中干擾信號比目標(biāo)信號更強(qiáng)的情況在實(shí)際觀測中較為常見，為了盡可能識別弱的目標(biāo)信號需進(jìn)行鏈?zhǔn)綑z驗(yàn)增點(diǎn)處理。經(jīng)此處理后仍可能存在具有一定連續(xù)性的干擾信號，為抑制此類干擾信號需進(jìn)行鏈?zhǔn)綑z驗(yàn)減點(diǎn)，步驟如下：

第1步，以判斷徑向速度梯度以及信號譜寬比值大小的方式實(shí)現(xiàn)鏈?zhǔn)綑z驗(yàn)增點(diǎn)?；诶走_(dá)某一波束所有距離庫，將所有三峰視作準(zhǔn)目標(biāo)信號。計(jì)算某距離庫中每個(gè)信號與其相鄰兩個(gè)距離庫中所有準(zhǔn)目標(biāo)信號的徑向速度之差以及譜寬之比，即該距離庫中每個(gè)信號對應(yīng)多個(gè)徑向速度差值以及譜寬比值。將每個(gè)信號對應(yīng)的每個(gè)徑向速度差值與當(dāng)?shù)靥卣髦颠M(jìn)行比較，當(dāng)小于當(dāng)?shù)靥卣髦禂?shù)量大于1時(shí)，繼續(xù)下一步的處理，否則視為非目標(biāo)信號。徑向速度之差的特征值由北京探空數(shù)據(jù)確定，且徑向速度梯度特征值的精確度對此步驟的影響很小。判斷譜寬比值的大小，小于某一閾值時(shí)則視此信號為準(zhǔn)目標(biāo)信號。在處理北京MST雷達(dá)數(shù)據(jù)時(shí)，譜寬比的閾值取0.5時(shí)效果好。進(jìn)行譜寬大小的判斷是因?yàn)槿魞H以徑向速度差值為判斷標(biāo)準(zhǔn)，噪聲信號有很大干擾；由于噪聲信號譜寬比其他信號譜寬偏窄，因此可通過增加譜寬大小的判斷實(shí)現(xiàn)鏈?zhǔn)皆鳇c(diǎn)的同時(shí)達(dá)到噪聲信號抑制的目的。經(jīng)過鏈?zhǔn)綑z驗(yàn)增點(diǎn)處理后，得到的準(zhǔn)目標(biāo)信號則主要包含目標(biāo)信號以及有一定空間連續(xù)性的雜波和降水信號等。

第2步，以判斷連續(xù)信號鏈上下兩端距離庫與各自相鄰距離庫準(zhǔn)目標(biāo)回波信號的徑向速度差值大小的方式實(shí)現(xiàn)鏈?zhǔn)綑z驗(yàn)減點(diǎn)。經(jīng)過鏈?zhǔn)綑z驗(yàn)增點(diǎn)后，將準(zhǔn)目標(biāo)信號連接成鏈，在某些高度范圍可能存在多條鏈。對于具有空間連續(xù)性的非目標(biāo)信號所連接成的鏈，通常僅在鏈的兩端存在極大的多普勒速度變化，即鏈兩端距離庫與各自臨近距離庫的徑向速度差值有異常大值，而目標(biāo)信號鏈在整個(gè)波束探測的距離庫范圍內(nèi)連續(xù)，不存在這樣的特點(diǎn)，據(jù)此可以將非目標(biāo)信號的鏈剔除。實(shí)際應(yīng)用時(shí)，計(jì)算某一距離庫中某一準(zhǔn)目標(biāo)回波信號的徑向速度與上下各3層(共6層)距離庫中所有準(zhǔn)目標(biāo)回波信號的徑向速度差值，并記錄徑向速度差值絕對值的最小值，共有6個(gè)正值的序列。當(dāng)此序列的方差過大或值異常偏大的點(diǎn)數(shù)大于2，則認(rèn)為此信號為非目標(biāo)信號。徑向速度差值異常偏大的判斷標(biāo)準(zhǔn)是與當(dāng)?shù)靥卣髦颠M(jìn)行比較，又因異常偏大值遠(yuǎn)大于當(dāng)?shù)靥卣髦?，因此徑向速度梯度特征值精度對此步驟影響很小。本文算法以大氣晴空湍流為目標(biāo)信號，在后續(xù)研究中，也可以對其他信號鏈進(jìn)行提取，如降水造成的信號[38-39]。

2.2.3 對稱性檢驗(yàn)

對于具有五波束的MST雷達(dá)，其中東西波束和南北波束的方位角夾角均為180°，假設(shè)大氣水平風(fēng)場均勻，則認(rèn)為對稱波束在同高度距離庫上測得的徑向速度具有對稱性。探測數(shù)據(jù)也證明各對稱波束的目標(biāo)信號幾乎均關(guān)于原點(diǎn)對稱。

經(jīng)過鏈?zhǔn)綑z驗(yàn)后，仍可能有噪聲與地雜波等非目標(biāo)信號殘留,可利用對稱波束回波信號具有對稱性的特點(diǎn)，通過對稱性檢驗(yàn)繼續(xù)抑制剩余干擾信號。

同鏈?zhǔn)綑z驗(yàn)相比，雖然對稱性檢驗(yàn)也能抑制間歇性雜波以及降水信號等不具有對稱性的回波信號，但實(shí)際情況中對稱性檢驗(yàn)有一定局限性，如鏈?zhǔn)綑z驗(yàn)?zāi)芎芎靡种圃肼曅盘?，而直接進(jìn)行對稱性檢驗(yàn)則無法抑制噪聲信號;當(dāng)對稱波束中一個(gè)方向缺測時(shí)，無法進(jìn)行對稱檢驗(yàn);垂直波束無法使用對稱性檢驗(yàn)。因而鏈?zhǔn)綑z驗(yàn)是必不可少的前提工作。

2.2.4 地雜波抑制

對稱性檢驗(yàn)對地雜波具有一定抑制作用，但仍可能有地雜波殘留，如對稱波束均有地雜波干擾時(shí)，無法通過對稱性檢驗(yàn)抑制地雜波干擾?；诟咚箶M合判斷的思路，根據(jù)地雜波在零頻附近以及譜寬較窄的特點(diǎn)，可對殘余的地雜波再做抑制處理。

經(jīng)數(shù)據(jù)處理后的東西波束5個(gè)距離庫的功率譜密度如圖4所示，改進(jìn)算法能夠很好地抑制干擾信號，無需進(jìn)行時(shí)間連續(xù)性檢驗(yàn)。

圖4 2012年1月4日00:10北京MST雷達(dá)東波束與西波束經(jīng)改進(jìn)算法處理后5個(gè)距離庫高度功率譜密度(加粗曲線表示識別出的目標(biāo)回波信號)

3 基礎(chǔ)數(shù)據(jù)估算及算法有效性與可靠性檢驗(yàn)

各階矩參數(shù)計(jì)算方法主要有兩種：一種是高斯擬合法(Gaussian fitting method)[41]，另一種為加權(quán)矩法(the weighted moment method)。根據(jù)Yamamoto等[42]與Wilson等[43]的研究，高SNR時(shí)加權(quán)矩法估算結(jié)果更好，低SNR時(shí)則相反。經(jīng)改進(jìn)算法處理后，目標(biāo)信號具有高SNR，因此本文使用加權(quán)矩法。Woodman[44]給出了加權(quán)矩法的各階矩計(jì)算公式。

根據(jù)MST雷達(dá)五波束探測方式，若假設(shè)觀測區(qū)域內(nèi)的水平風(fēng)場均勻，則可以用多波束估算法進(jìn)行三維風(fēng)場計(jì)算。多波束估算法是將三波束估算法與五波束估算法進(jìn)行聯(lián)合使用的方法，何平[45]詳細(xì)介紹了對相關(guān)計(jì)算方法。

應(yīng)用改進(jìn)算法處理得到的估算結(jié)果、雷達(dá)產(chǎn)品數(shù)據(jù)以及探空數(shù)據(jù)廓線如圖5所示，改進(jìn)算法估算的水平風(fēng)速與風(fēng)向與探空數(shù)據(jù)有很好的一致性，改進(jìn)后算法估算的基礎(chǔ)數(shù)據(jù)，連續(xù)性明顯增強(qiáng)，異常值出現(xiàn)概率減小。由2012年1月1日—12月31日改進(jìn)算法與雷達(dá)原算法得到的SNR、譜寬(各11593組對比廓線)對比(圖6)可以看到，對于SNR，改進(jìn)算法得到的結(jié)果更符合實(shí)際，即隨高度增加SNR減小，直至為負(fù)數(shù)，而雷達(dá)原算法得到的結(jié)果全部為正值(圖6a)。雷達(dá)原算法由于未精確識別目標(biāo)回波，得到的譜寬明顯偏大，改進(jìn)算法得到的譜寬平均值為2.5 m·s-1，標(biāo)準(zhǔn)差也小于雷達(dá)原算法，所得譜寬更可信(圖6b)。

圖5 2012年3月16日11：40北京MST雷達(dá)原算法與改進(jìn)算法得到的水平風(fēng)速、水平風(fēng)向與探空(11:15)對比

圖6 2012年1月1日—12月31日信噪比(a)和譜寬(b)的年平均值(曲線)和標(biāo)準(zhǔn)差(陰影)廓線

為進(jìn)一步檢驗(yàn)改進(jìn)算法的有效性和可靠性，分別應(yīng)用改進(jìn)算法和雷達(dá)原算法得到2012年1月1日—12月31日的緯向風(fēng)、經(jīng)向風(fēng)與垂直速度。緯向風(fēng)和經(jīng)向風(fēng)同北京探空測值進(jìn)行對比(共504組比對廓線，見圖7)，垂直速度與ERA5再分析數(shù)據(jù)進(jìn)行比對(共5312組對比廓線，見圖8)。由圖7可知，在3～22 km高度范圍內(nèi)，改進(jìn)算法得到的緯向風(fēng)與探空測值的平均差值在-1～0 m·s-1范圍內(nèi)，與雷達(dá)原算法相比更接近探空測值。由圖8可知，在3～22 km高度范圍內(nèi)，改進(jìn)算法得到的垂直速度與ERA5再分析數(shù)據(jù)的平均差值和差值的標(biāo)準(zhǔn)差分別在-0.05～0.05 m·s-1和0.25～0.3 m·s-1范圍內(nèi)，明顯小于雷達(dá)原算法與ERA5數(shù)據(jù)的差異(平均差值和差值的標(biāo)準(zhǔn)差分別在-0.35～-0.25 m·s-1和0.4～0.5 m·s-1范圍內(nèi))。改進(jìn)算法得到的垂直速度與ERA5數(shù)據(jù)的差異明顯小于雷達(dá)產(chǎn)品與ERA5數(shù)據(jù)的差異。

圖7 2012年1月1日—12月31日北京MST雷達(dá)原算法與改進(jìn)算法得到的經(jīng)向風(fēng)、緯向風(fēng)與探空對比

圖8 2012年1月1日—12月31日北京MST雷達(dá)原算法與改進(jìn)算法得到的垂直速度與ERA5再分析數(shù)據(jù)對比

改進(jìn)算法可更準(zhǔn)確地識別大氣湍流信號，有效抑制雜波及降水信號的干擾。為檢驗(yàn)改進(jìn)算法在降水天氣條件下的有效性，將2012年共16組降水天氣條件下(根據(jù)雷達(dá)站所在地歷史天氣記錄全天有降水時(shí)段，7月25，26，28，31日，8月3日及12月13日，20日)的探空與北京MST雷達(dá)原算法及改進(jìn)算法得到的經(jīng)向風(fēng)和緯向風(fēng)廓線進(jìn)行統(tǒng)計(jì)，結(jié)果如圖9所示。由圖9可知，降水條件下改進(jìn)算法得到的經(jīng)向風(fēng)和緯向風(fēng)與探空測值的平均偏差及均方根誤差均小于雷達(dá)原算法與探空測值比對結(jié)果，驗(yàn)證了改進(jìn)算法在降水天氣條件下的有效性和可靠性。

圖9 2012年16組降水天氣條件下北京MST雷達(dá)原算法與改進(jìn)算法得到的經(jīng)向風(fēng)、緯向風(fēng)與探空測值平均差值和均方根誤差廓線

綜上所述，基于2012年多種數(shù)據(jù)對比分析，使用改進(jìn)算法得到的經(jīng)向風(fēng)、緯向風(fēng)與探空測值更接近，垂直速度、SNR及譜寬數(shù)據(jù)可靠性更高，同時(shí)改進(jìn)算法在降水條件下依然有效。

4 結(jié)論與討論

本文從噪聲電平估算和目標(biāo)回波準(zhǔn)確識別與提取兩方面對北京MST雷達(dá)功率譜處理算法進(jìn)行改進(jìn)，并驗(yàn)證算法的有效性與可靠性。改進(jìn)算法可用于北京MST雷達(dá)歷史觀測數(shù)據(jù)，結(jié)果表明：

1)對數(shù)-線性擬合方案估算得到的噪聲電平與二分法方案估算值一致，二者差值的標(biāo)準(zhǔn)差、平均絕對差值與相關(guān)系數(shù)分別是0.43 dB，0.21 dB和0.998，前者的平均值為168.6 dB，后者的平均值為168.5 dB。對數(shù)-線性擬合方案可快速、準(zhǔn)確估算高斯白噪聲的噪聲電平值，利于后續(xù)目標(biāo)回波信號的準(zhǔn)確識別且業(yè)務(wù)實(shí)用性強(qiáng)。

2)對數(shù)-線性擬合方案具有兩個(gè)明顯優(yōu)勢：對數(shù)-線性擬合方案不需要預(yù)設(shè)參數(shù)，使算法更加簡便；對數(shù)-線性擬合方案能夠判斷功率譜密度中能否提取出目標(biāo)回波信號，是否需要進(jìn)行后續(xù)數(shù)據(jù)處理步驟，提高計(jì)算速度與效率。

3)北京MST雷達(dá)使用改進(jìn)算法后，非目標(biāo)回波被有效抑制，目標(biāo)回波可被準(zhǔn)確識別與提取，可在功率譜密度圖中更精確地顯示大氣湍流信號。

4)在3～22 km高度范圍內(nèi)，改進(jìn)算法處理結(jié)果與探空緯向風(fēng)均方根誤差為2～3 m·s-1，而雷達(dá)產(chǎn)品與探空緯向風(fēng)均方根誤差為3～4 m·s-1。此外，改進(jìn)算法得出的譜寬平均值為2.5 m·s-1，小于雷達(dá)產(chǎn)品平均值，改進(jìn)算法得到的垂直速度與ERA5數(shù)據(jù)差異明顯小于雷達(dá)產(chǎn)品與ERA5數(shù)據(jù)的差異。

5)在降水天氣條件下，改進(jìn)算法處理結(jié)果和探空測值在不同高度的水平風(fēng)速平均偏差和均方根誤差均小于雷達(dá)產(chǎn)品和探空差值。

從本文分析可知，改進(jìn)算法有效可靠，明顯改善雷達(dá)數(shù)據(jù)質(zhì)量且易于實(shí)現(xiàn)。該算法的原理與實(shí)現(xiàn)方法具有普適性，可用于有相同探測機(jī)制的各類MST雷達(dá)、ST雷達(dá)和風(fēng)廓線雷達(dá)的數(shù)據(jù)處理，但算法的有效性與可靠性需進(jìn)一步檢驗(yàn)。本文是基于信號處理方式固定條件下改進(jìn)了數(shù)據(jù)處理算法，若將信號處理與數(shù)據(jù)處理的改進(jìn)相結(jié)合，則可更有效地提高雷達(dá)數(shù)據(jù)質(zhì)量。

致 謝:本文使用了國家重大科技基礎(chǔ)設(shè)施“子午工程”科學(xué)數(shù)據(jù)。感謝中國科學(xué)院大氣物理研究所香河大氣綜合觀測試驗(yàn)站工作人員對MST雷達(dá)運(yùn)維所做的工作。