微波濕度計(jì)條帶噪聲的分析與抑制方法研究

劉明緒，張升偉，何杰穎

1.中國科學(xué)院國家空間科學(xué)中心 微波遙感重點(diǎn)實(shí)驗(yàn)室,北京 100190;

2.中國科學(xué)院大學(xué),北京 100049

1 引言

衛(wèi)星氣象數(shù)據(jù)的降噪一直是遙感氣象數(shù)據(jù)處理領(lǐng)域研究的重點(diǎn)問題。Zou 等（2012）在處理“風(fēng)云三號(hào)”A星、B星（FY-3A、FY-3B）上的微波濕度計(jì)（MWHS）和美國NOAA系列衛(wèi)星上的微波濕度計(jì)（Microwave Humidity Sounder，MHS）的數(shù)據(jù)時(shí)，觀測(cè)到了一種沿軌道方向隨掃描角變化的直線型噪聲，且噪聲不可通過平均消除。因此Zou 等（2012）提出了一種通過主成分分析（PCA）與五點(diǎn)平滑法抑制該噪聲的方法。Qin 等（2013）分析了SNPP（Suomi National Polar-orbiting Partnership）衛(wèi)星上ATMS（Advanced Technology Microwave Sounder）水汽探測(cè)通道觀測(cè)數(shù)據(jù)中的條帶噪聲，并通過主成分分析聯(lián)合集合經(jīng)驗(yàn)?zāi)B(tài)分解的方法抑制了該噪聲。Li 和Liu（2016）在對(duì)MWTS-2 輻射同化的質(zhì)量控制程序中引入了該條帶抑制方法。對(duì)比發(fā)現(xiàn)，去除帶噪聲后分析誤差小于不修正帶噪聲吸收數(shù)據(jù)，這表明該降噪方法的實(shí)際應(yīng)用價(jià)值。Zou 等（2017）在處理ATMS 的條帶噪聲時(shí)，對(duì)原PCA/EEMD 方法進(jìn)行了修改，避免了在海陸邊界亮溫急劇變化導(dǎo)致的降噪后的誤差。同時(shí)，Zou 和Tian（2019）利用MWTS-2 掃描周期變化前后數(shù)據(jù)的差異對(duì)條帶噪聲進(jìn)行了分析，表明噪聲產(chǎn)生的根本原因不在于掃描速度，并給出了一種條帶噪聲的分析指標(biāo)。金旭等（2019）通過計(jì)算噪聲分布與方差，定量評(píng)估了風(fēng)云三號(hào)微波溫度計(jì)條帶噪聲的抑制效果。

目前，對(duì)于PCA/EEMD 降噪方法的各類改進(jìn)、分析與應(yīng)用均基于EEMD 的基礎(chǔ)上進(jìn)行。EEMD 的理論基礎(chǔ)之一是噪聲輔助數(shù)據(jù)分析NADA（Noise-Assisted Data Analysis）（Wu 和Huang，2009），因而該方法可能受到輔助噪聲的影響，而導(dǎo)致殘余噪聲與虛假分量等問題。模態(tài)分解作為時(shí)頻分析領(lǐng)域一種常見方法，已有諸多學(xué)者提出了如變分模態(tài)分解VMD（Variational Mode Decomposition）、多元經(jīng)驗(yàn)?zāi)B(tài)分解MEMD（Multivariate Empirical Mode Decomposition）、改進(jìn)的自適應(yīng)噪聲總體集合經(jīng)驗(yàn)?zāi)B(tài)分解ICEEMDAN（Improved Complete Ensemble Empirical Mode Decomposition With Adaptive Noise）等變體或改進(jìn)算法。其中，互補(bǔ)集合經(jīng)驗(yàn)?zāi)B(tài)分解CEEMD（Complementary Ensemble Empirical Mode Decomposition）、自適應(yīng)噪聲總體集合經(jīng)驗(yàn)?zāi)B(tài)分解CEEMDAN（Complete Ensemble Empirical Mode Decomposition With Adaptive Noise）及ICEEMDAN 等方法是一類試圖改進(jìn)并降低輔助噪聲影響的方法，可有效降低重構(gòu)信號(hào)后的噪聲。本文通過對(duì)各類方法進(jìn)行比較，試圖改進(jìn)原有的EEMD 方法，并應(yīng)用于MWHS-2 的數(shù)據(jù)處理中。同時(shí)，本文針對(duì)風(fēng)云三號(hào)微波濕度計(jì)探測(cè)數(shù)據(jù)與模態(tài)分解的特有性質(zhì)，給出了一種能有效判別噪聲和信號(hào)模態(tài)的新方法，可用于風(fēng)云三號(hào)微波濕度計(jì)探測(cè)數(shù)據(jù)噪聲特征的分析，并有助于提升分解的準(zhǔn)確性，提升數(shù)據(jù)質(zhì)量。

2 降噪算法原理

2.1 主成分分析原理

主成分分析PCA（Principal Components Analysis）是一種常見的數(shù)據(jù)分解與降維的方法。該方法基于最大化數(shù)據(jù)方差原則，將數(shù)據(jù)進(jìn)行投影到低維超平面上，實(shí)現(xiàn)對(duì)數(shù)據(jù)的降維。Lorenz（1956）將PCA 方法引入氣象數(shù)據(jù)分析領(lǐng)域，稱為經(jīng)驗(yàn)正交函數(shù)法EOF（Empirical Orthogonal Function）。該方法將觀測(cè)數(shù)據(jù)分解為空間模態(tài)和時(shí)間系數(shù)兩部分，從空間和時(shí)間域?qū)?shù)據(jù)進(jìn)行分析（金旭等，2019）。

設(shè)MWHS觀測(cè)亮溫矩陣為

式中，M表示每條掃描線的視場（FOV）點(diǎn)，N為掃描線個(gè)數(shù)，則Tij為第j條掃描線中第i個(gè)視場的觀測(cè)亮溫值。

根據(jù)PCA 的原理，對(duì)亮溫矩陣A的協(xié)方差矩陣AAT做特征分解，即

λi為特征值，為對(duì)應(yīng)特征向量。將上式寫作矩陣形式，有

E為特征向量組成矩陣，表示了使掃描線方差最大的投影方向。此外，定義

該基下，原矩陣A可表示為

表示A在第i主成分中的分量信息，稱為第i模態(tài)。

2.2 模態(tài)分解原理

經(jīng)驗(yàn)?zāi)?態(tài)分解 EMD（Empirical Mode Decomposition）是由Huang 等（1998）提出的一種針對(duì)非線性非平穩(wěn)信號(hào)的自適應(yīng)的時(shí)頻分析方法。該方法將信號(hào)分解為若干個(gè)本征模態(tài)函數(shù)IMF（Intrinsic Mode Function），每一IMF 均包含有信號(hào)的不同時(shí)間分辨率的信息。

EMD大體步驟為

（1）根據(jù)原始信號(hào)x(t)的所有極大值極小值點(diǎn)，利用三次樣條插值擬合出上下包絡(luò)線。

（2）求上下包絡(luò)線的均值m1(t)，用原信號(hào)減去該均值，得到一新序列

（3）判斷中間信號(hào)h1(t)是否滿足IMF 的條件，若滿足，則作為第一個(gè)IMF 分量IMF1(t)，否則將h1(t)作為待分解信號(hào)，重復(fù)步驟（1）—（3），直至滿足為止。

（4）用上面的方法得到第一個(gè)IMF 后，將該IMF從原始信號(hào)中分離出來，得到

（5）將r1(t)作為新的信號(hào)重復(fù)步驟（1）—（4），直到最終分解得到的rn(t)為一單調(diào)函數(shù)為止。此時(shí)分解得到的IMF1(t)，…，IMFn(t)即為所求的n個(gè)模態(tài)。

此后，為解決EMD 在處理跳變或間斷時(shí)可能產(chǎn)生的模態(tài)混疊問題，Wu和Huang等（2009）提出了EEMD（Ensemble Empirical Mode Decomposition）方法。該方法通過在分解前添加白噪聲解決了EMD 中的模態(tài)混疊問題，之后通過多次計(jì)算取平均去除白噪聲的影響。

EEMD步驟如下：

已知信號(hào)為x(t)，信號(hào)中分別加入白噪聲n(t)，得到

對(duì)s(t)進(jìn)行EMD 分解，得到所求IMF 分量IMF1(t)，…，IMFn(t)。

重復(fù)分解m次，將每一分量所得m個(gè)IMF進(jìn)行平均，最終每一IMF 分量即為m組IMF 平均的結(jié)果。

EEMD 方法雖然有效的解決了EMD 中的模態(tài)混疊，但也產(chǎn)生了一些新的問題。例如，在平均去噪時(shí)，由于平均次數(shù)有限，噪聲不可能被完全中和，導(dǎo)致重構(gòu)后的信號(hào)存在殘余噪聲。此外，由于每次添加的噪聲隨機(jī)，最終分解出的IMF 數(shù)量的可能會(huì)受到噪聲的影響，因此平均的各IMF 并不一定對(duì)齊，而導(dǎo)致了虛假分量的產(chǎn)生。

CEEMD（Complementary Ensemble Empirical Mode Decomposition）方法（Yeh 等，2010）將EEMD 中添加的單個(gè)白噪聲替換為一組符號(hào)相反的白噪聲，即

s1(t)、s2(t)分別為加入正負(fù)噪聲后的信號(hào)。

對(duì)s1(t)、s2(t)分別進(jìn)行EMD 分解，最終將正負(fù)兩組結(jié)果取平均，完成分解。該方法有效減少了EEMD分解后可能存在的殘余噪聲。

CEEMDAN（Complete Ensemble Empirical Mode Decomposition With Adaptive Noise）方法由Torres等（2011）在EEMD 的基礎(chǔ)上，在EMD 計(jì)算出每個(gè)IMF 后不使用原加噪信號(hào)s(t)，而使用噪聲對(duì)應(yīng)IMF 分量x(t) +Ej(n(t))（Ej(·)表示EMD 分解后的第j模態(tài)）求殘差，從而降低了輔助噪聲的影響。同時(shí)，該方法在分解得到IMF 后立即進(jìn)行平均得 到，用平均后的IMF 求殘差rj(t)，即

相較于先前先求所有IMF 再平均的方法，該方法徹底避免了模態(tài)數(shù)不同難以平均的問題。

之后，Colominas等（2014）又基于CEEMDAN提出了ICEEMDAN（Improved Complete Ensemble Empirical Mode Decomposition with Adaptive Noise）。先前方法提取IMF 的過程均是對(duì)信號(hào)添加噪聲并進(jìn)行多次迭代最終得到IMF，即

M(·)為迭代完成后移除的包絡(luò)線均值總和，為均值算子。

可以看出，最終得到的IMF中包含有加噪信號(hào)和包絡(luò)均值兩部分平均后的殘余噪聲。ICEEMDAN在求IMF時(shí)使用原始信號(hào)去除包絡(luò)均值

改進(jìn)后的方法噪聲只存在于包絡(luò)均值中，去掉了原信號(hào)中的噪聲，從而有效減少了最終平均后殘余的噪聲。此外，ICEEMDAN 對(duì)初始信號(hào)不再添加白噪聲而直接使用噪聲模態(tài)，從而避免了CEEMDAN中偽模式的問題。

各方法關(guān)系如圖1所示。

圖1 模態(tài)分解算法關(guān)系Fig.1 The relationship between modal decomposition algorithms

3 降噪方法與結(jié)果分析

3.1 降噪方法

測(cè)試所用觀測(cè)亮溫?cái)?shù)據(jù)選取2016 年10 月8 日與2019 年1 月10 日濕度計(jì)采集的各14 組亮溫?cái)?shù)據(jù)進(jìn)行測(cè)試。其中，118.75 GHz 附近的第4 通道（118.75±0.3 GHz）有較為明顯的條帶噪聲。模擬亮溫?cái)?shù)據(jù)由對(duì)應(yīng)時(shí)間的ERA5每小時(shí)全球再分析分層網(wǎng)格數(shù)據(jù)與ERA5-land每小時(shí)地面網(wǎng)格數(shù)據(jù)導(dǎo)入RTTOV-SCATT 模塊仿真生成。網(wǎng)格分辨率0.25°×0.25°，氣壓層選取50 hPa—1000 hPa間共29層。

部分區(qū)域的觀測(cè)亮溫如圖2（a）所示?？梢钥闯?，亮溫圖像中有明顯的掃描線方向的條紋狀噪聲。

圖2 條帶噪聲去除前后亮溫圖及（a）、（b）間差異Fig.2 Brightness temperature before and after removing the striping noise and difference between（a）and（b）

對(duì)數(shù)據(jù)進(jìn)行PCA 分解，得到其各個(gè)PC 分量（圖3），可以看出第一PC 模態(tài)中有明顯的水平條紋。根據(jù)PCA 的原理，第一主成分方向?yàn)樾盘?hào)差異最大的方向，而觀測(cè)亮溫中的差異主要來源于衛(wèi)星軌跡方向宏觀的天氣變化，因此第一主成分分量包含了最多的沿衛(wèi)星軌跡方向包括條帶噪聲在內(nèi)的天氣信息。通過計(jì)算方差貢獻(xiàn)率

圖3 第1—4 PC模態(tài)空間分布圖Fig.3 Spatial distribution of 1-4 PC modes

第一主成分方差貢獻(xiàn)比例超過99.98%，表明絕大部分的沿軌方向數(shù)據(jù)差異都成功保留在了第一主成分分量中。后續(xù)模態(tài)包含了沿掃描線變化及各類小尺度波動(dòng)信息。因此，僅對(duì)進(jìn)行模態(tài)分解，也可較為準(zhǔn)確且高效的提取出絕大部分條帶噪聲，且可有效降低處理的數(shù)據(jù)量，提高計(jì)算效率。

分解后所得前7 個(gè)IMF 如圖4?？梢钥闯?，前幾個(gè)IMF主要包含了信號(hào)中的高頻噪聲，后續(xù)模態(tài)保留了各類天氣信息。分解后的PC分量可表示為

IMFn和IMFs分別為噪聲主導(dǎo)和信號(hào)主導(dǎo)的IMF分量，r(t)為提取完所有IMF后的直流分量。

因此，只需選擇合適的劃分位置，分離出噪聲主導(dǎo)的前幾個(gè)IMF分量并去除，即可得到重構(gòu)信號(hào)

定義信號(hào)能量密度為

式中，N為IMF 長度，IMFj(i)表示第j個(gè)IMF 序列的第i個(gè)元素。

由于IMF 極值點(diǎn)與零點(diǎn)相同（或差一）的定義，IMF可近似看作一周期振蕩的函數(shù)。因此，可利用其極值點(diǎn)個(gè)數(shù)nl定義IMF 的平均周期T（Wu和Huang，2004）

同時(shí)，Wu 和Huang（2004）指出，對(duì)于白噪聲，IMF函數(shù)S(lnT，j)在對(duì)數(shù)坐標(biāo)下積分后的數(shù)值相同（c為常數(shù)）。

該式等價(jià)于

對(duì)于本文所述條帶噪聲，雖無法證明其為白噪聲，但觀察各IMF的頻譜（圖5）可以看出，前5個(gè)IMF 峰值明顯低于后續(xù)IMF，且除包含殘余噪聲的第一IMF外，其余幾個(gè)IMF具有相似的波形，這與Wu和Huang（2009）的實(shí)驗(yàn)結(jié)果近似。因此，不妨嘗試計(jì)算各IMF能量密度與平均周期的乘積EnTn并計(jì)算相鄰IMF間能量之比。在噪聲IMF之間，該比值較小，而在噪聲與信號(hào)間，比值會(huì)突變。可尋找出比值變化最大的間隔作為信號(hào)與噪聲的劃分P：

圖5 第一PC分量IMF1-10頻譜圖Fig.5 Fourier spectra of the first PC coefficient IMF1-10

同時(shí)，為避免對(duì)無噪聲通道的誤處理，在尋找出最大比值后還需對(duì)其進(jìn)行約束。比較有噪聲通道和無噪聲通道的頻譜發(fā)現(xiàn)，無噪聲通道并無明顯的能量突變。因此可定義僅當(dāng)最大突變超過某一閾值時(shí)才將其定義為噪聲并進(jìn)行劃分。對(duì)于本文，閾值定義為0.8，實(shí)際操作過程中，可以實(shí)現(xiàn)較好的劃分效果。

在模態(tài)分解方法的選取上，如2.2 節(jié)所述，EEMD 會(huì)因殘余噪聲和虛假分量而在重構(gòu)后產(chǎn)生誤差，使得天氣信號(hào)分量受到噪聲的影響。同時(shí)，EEMD平均時(shí)模態(tài)數(shù)不一定對(duì)齊，致使分解后的直流分量也混雜到了IMF 中，使噪聲選取方法失效。因此本文采用了ICEEMDAN 方法進(jìn)行處理，并最后重構(gòu)出亮溫信號(hào)

3.2 結(jié)果分析

降噪前后的亮溫圖像及提取出的噪聲如圖4（b）（c）?？梢钥闯觯惴ǔ晒μ崛〕隽诵盘?hào)中的條帶噪聲并重構(gòu)出了觀測(cè)亮溫。對(duì)比圖6可以看出，重構(gòu)后的信號(hào)O-B 亮溫并無明顯偏移與失真，該方法恢復(fù)后的數(shù)據(jù)中保留了信號(hào)中絕大多數(shù)的亮溫信息。其中單獨(dú)提取出第49 視場角并去除前n個(gè)IMF后繪制的亮溫曲線如圖7。可以看出，隨著去除的IMF 數(shù)增加，信號(hào)逐漸趨于平滑，且各級(jí)模態(tài)分解均主要抑制的是信號(hào)中的高頻振動(dòng)，對(duì)信號(hào)整體趨勢(shì)變化并無太大影響。同時(shí)，對(duì)提取出的噪聲進(jìn)行統(tǒng)計(jì)，其概率分布（圖8）整體呈現(xiàn)出近似于無偏的高斯噪聲的分布特征。

圖6 噪聲抑制前后O-B亮溫圖Fig.6 O-B before and after removing the striping noise

圖7 去除前1—5 IMF后重構(gòu)出的第49視場角亮溫隨掃描線變化圖Fig.7 Reconstruction of the 49th FOV brightness temperature changes with the scan line after remove the first 1—5 IMF

圖8 條帶噪聲概率分布直方圖Fig.8 Histogram of probability distribution of striping noise

表1 各模態(tài)分解方法降噪效果Table 1 Mitigation effect of each modal decomposition method

可以看出，各類方法相較于原始信號(hào)均在不同程度上對(duì)噪聲進(jìn)行了抑制。信噪比越高效果越好，方差和均方誤差相對(duì)越低越好。對(duì)比各類方法，可以看出相較于原始的EEMD 算法，改進(jìn)后的CEEMD 等方法在去噪效果上均有不同程度的提高。其中ICEEMDAN 在各方面效果最好。相較于EEMD，該方法使信噪比提高了0.031 dB，方差降低0.020 K2。且多次實(shí)驗(yàn)表明該方法具有較好的穩(wěn)定性，有助于提高算法的噪聲抑制效果。

4 結(jié)論

本文主要介紹了針對(duì)MWHS-2 中的條帶噪聲的一種有效的降噪方法，并給出了改進(jìn)方案。本文首先使用主成分分析對(duì)數(shù)據(jù)進(jìn)行了降維，隨后通過模態(tài)分解將信號(hào)分解包含噪聲和天氣信息的不同IMF分量，利用兩類IMF在能量密度上的差異進(jìn)行自動(dòng)劃分，重構(gòu)剩余模態(tài)，實(shí)現(xiàn)了噪聲的抑制。針對(duì)EEMD 殘余噪聲等問題，本文采用ICEEMDAN 進(jìn)行PC 分量的分解、降噪與重構(gòu)操作。在對(duì)MHWS-2 數(shù)據(jù)的處理過程中，證明了該方法的有效性。同時(shí)，本文針對(duì)各種不同的模態(tài)分解方法的降噪效果進(jìn)行了統(tǒng)計(jì)與對(duì)比，數(shù)據(jù)結(jié)果表明，各類改進(jìn)方法均可在不同程度上實(shí)現(xiàn)噪聲的抑制。其中，使用ICEEMDAN 進(jìn)行模態(tài)分解降噪并恢復(fù)后的數(shù)據(jù)使信噪比提高了0.054 dB，方差降低0.032 K2，具有良好的重建效果，可使重構(gòu)后的數(shù)據(jù)精度得到進(jìn)一步提高。

模態(tài)分解是經(jīng)驗(yàn)上近似的二進(jìn)濾波器，即以相鄰模態(tài)間周期大致為2 倍的關(guān)系對(duì)信號(hào)進(jìn)行分解，并不一定能將噪聲和信號(hào)準(zhǔn)確的劃分開。直接對(duì)模態(tài)進(jìn)行刪除只能濾除絕大多數(shù)噪聲，并不能保證剩余模態(tài)中沒有噪聲殘余。同時(shí)，本文以能量差異進(jìn)行劃分的方法，只針對(duì)所用微波濕度計(jì)數(shù)據(jù)進(jìn)行了測(cè)試，僅作為一種思路以供參考，對(duì)于其適用范圍的普遍性和對(duì)各類數(shù)據(jù)劃分的穩(wěn)定性，還有待進(jìn)一步的探索。

志 謝文中所用濕度計(jì)相關(guān)數(shù)據(jù)來自國家衛(wèi)星氣象中心，在此表示衷心的感謝。