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

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

引 言

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

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

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

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

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

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探測(cè)1次，垂直分辨率為600 m。

北京探空數(shù)據(jù)：選用2012年1月1日—12月31日北京市觀象臺(tái)探空的水平風(fēng)場(chǎng)數(shù)據(jù)，每日兩次探測(cè)，探空球施放時(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á)信號(hào)并非只有目標(biāo)信號(hào)，還存在非目標(biāo)信號(hào)，因此對(duì)功率譜密度進(jìn)行處理，精確提取目標(biāo)信號(hào)是一項(xiàng)必需的工作。MST雷達(dá)信號(hào)可分為兩大類[37]，即大氣要素與非大氣要素信號(hào)。大氣要素信號(hào)包括晴空大氣湍流、降水等，非大氣要素則包括地雜波(如地面的建筑、山、樹等造成的信號(hào))、間歇性雜波(如飛機(jī)、鳥類的飛行等所造成的信號(hào))和系統(tǒng)噪聲等。當(dāng)以大氣湍流為探測(cè)目標(biāo)，則噪聲、間歇性雜波等非目標(biāo)信號(hào)常對(duì)基礎(chǔ)數(shù)據(jù)估算有很強(qiáng)干擾。

2.1 噪聲電平估算

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

為抑制直流與窄而強(qiáng)的回波信號(hào)的影響，在處理功率譜密度前需進(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á)測(cè)量的噪聲信號(hào)屬于高斯白噪聲，將功率譜密度在某個(gè)待定噪聲閾值以下的點(diǎn)組成該距離庫的噪聲功率譜密度序列記為{Sm}，頻率范圍為F，第m點(diǎn)對(duì)應(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)對(duì)應(yīng)的頻率，共有M個(gè)點(diǎn)，{Sm}為噪聲功率譜密度序列，F(xiàn)為頻率范圍，p是數(shù)據(jù)預(yù)處理時(shí)滑動(dòng)平均的點(diǎn)數(shù)，如果不進(jìn)行滑動(dòng)平均處理，則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)過比對(duì)客觀分析法與八段平均法的估算結(jié)果，發(fā)現(xiàn)八段平均法計(jì)算的噪聲電平較低(圖1)。對(duì)八段平均法，在功率譜密度圖上，高于噪聲電平的某些包絡(luò)會(huì)出現(xiàn)多個(gè)極大值，其中部分代表噪聲信號(hào)的包絡(luò)譜寬值異常偏大，不利于后續(xù)處理，而客觀分析法在此方面明顯優(yōu)于八段平均法。

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

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

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

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

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

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

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

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

2.2.1 信號(hào)提取

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

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

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

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

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

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

2.2.3 對(duì)稱性檢驗(yàn)

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

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

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

2.2.4 地雜波抑制

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

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

圖4 2012年1月4日00:10北京MST雷達(dá)東波束與西波束經(jīng)改進(jìn)算法處理后5個(gè)距離庫高度功率譜密度(加粗曲線表示識(shí)別出的目標(biāo)回波信號(hà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)信號(hào)具有高SNR，因此本文使用加權(quán)矩法。Woodman[44]給出了加權(quán)矩法的各階矩計(jì)算公式。

根據(jù)MST雷達(dá)五波束探測(cè)方式，若假設(shè)觀測(cè)區(qū)域內(nèi)的水平風(fēng)場(chǎng)均勻，則可以用多波束估算法進(jìn)行三維風(fēng)場(chǎng)計(jì)算。多波束估算法是將三波束估算法與五波束估算法進(jìn)行聯(lián)合使用的方法，何平[45]詳細(xì)介紹了對(duì)相關(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組對(duì)比廓線)對(duì)比(圖6)可以看到，對(duì)于SNR，改進(jìn)算法得到的結(jié)果更符合實(shí)際，即隨高度增加SNR減小，直至為負(fù)數(shù)，而雷達(dá)原算法得到的結(jié)果全部為正值(圖6a)。雷達(dá)原算法由于未精確識(shí)別目標(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)對(duì)比

圖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)同北京探空測(cè)值進(jìn)行對(duì)比(共504組比對(duì)廓線，見圖7)，垂直速度與ERA5再分析數(shù)據(jù)進(jìn)行比對(duì)(共5312組對(duì)比廓線，見圖8)。由圖7可知，在3～22 km高度范圍內(nèi)，改進(jìn)算法得到的緯向風(fēng)與探空測(cè)值的平均差值在-1～0 m·s-1范圍內(nèi)，與雷達(dá)原算法相比更接近探空測(cè)值。由圖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)與探空對(duì)比

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

改進(jìn)算法可更準(zhǔn)確地識(shí)別大氣湍流信號(hào)，有效抑制雜波及降水信號(hào)的干擾。為檢驗(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)與探空測(cè)值的平均偏差及均方根誤差均小于雷達(dá)原算法與探空測(cè)值比對(duì)結(jié)果，驗(yàn)證了改進(jìn)算法在降水天氣條件下的有效性和可靠性。

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

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

4 結(jié)論與討論

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

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

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

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

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é)果和探空測(cè)值在不同高度的水平風(fēng)速平均偏差和均方根誤差均小于雷達(dá)產(chǎn)品和探空差值。

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

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