劉明緒，張升偉，何杰穎

1.中國科學院國家空間科學中心 微波遙感重點實驗室,北京 100190;

2.中國科學院大學,北京 100049

1 引言

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

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

2 降噪算法原理

2.1 主成分分析原理

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

設MWHS觀測亮溫矩陣為

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

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

λi為特征值，為對應特征向量。將上式寫作矩陣形式，有

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

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

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

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

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

EMD大體步驟為

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

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

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

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

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

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

EEMD步驟如下：

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

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

重復分解m次，將每一分量所得m個IMF進行平均，最終每一IMF 分量即為m組IMF 平均的結果。

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

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

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

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

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

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

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

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

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

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

各方法關系如圖1所示。

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

3 降噪方法與結果分析

3.1 降噪方法

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

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

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

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

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

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

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

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

因此，只需選擇合適的劃分位置，分離出噪聲主導的前幾個IMF分量并去除，即可得到重構信號

定義信號能量密度為

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

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

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

該式等價于

對于本文所述條帶噪聲，雖無法證明其為白噪聲，但觀察各IMF的頻譜（圖5）可以看出，前5個IMF 峰值明顯低于后續(xù)IMF，且除包含殘余噪聲的第一IMF外，其余幾個IMF具有相似的波形，這與Wu和Huang（2009）的實驗結果近似。因此，不妨嘗試計算各IMF能量密度與平均周期的乘積EnTn并計算相鄰IMF間能量之比。在噪聲IMF之間，該比值較小，而在噪聲與信號間，比值會突變?？蓪ふ页霰戎底兓畲蟮拈g隔作為信號與噪聲的劃分P：

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

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

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

3.2 結果分析

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

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

圖7 去除前1—5 IMF后重構出的第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

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

4 結論

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

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

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