張新文，林冠英，劉同木，周保成

（1.自然資源部南海調(diào)查中心，廣東 廣州 510300；2.自然資源部海洋環(huán)境探測技術(shù)與應(yīng)用重點實驗室，廣東 廣州 510300）

海洋浮標是利用無動力漂浮載體獲取海洋環(huán)境信息的自動化監(jiān)測系統(tǒng)，具有“海洋上的地球同步衛(wèi)星”和“海上天氣偵察兵”等稱號，是海洋環(huán)境監(jiān)測、海洋資源調(diào)查評估和海洋科學(xué)研究的重要平臺[1-2]。受海洋環(huán)境噪聲、海洋生物附著和傳感器元件老化等因素的影響，浮標搭載的傳感器測量值不可避免地存在環(huán)境噪聲，直接影響數(shù)據(jù)應(yīng)用效果[3]。當前，海洋浮標受到衛(wèi)星通信帶寬限制，一般在浮標端需要對原始數(shù)據(jù)按照固定時間將觀測數(shù)據(jù)融合后回傳。因此，在浮標端實時開展數(shù)據(jù)異常值檢測和測量噪聲預(yù)處理，是提高數(shù)據(jù)融合精度的可靠手段之一。

國內(nèi)外學(xué)者在環(huán)境監(jiān)測數(shù)據(jù)預(yù)處理方面進行了相關(guān)研究，主要分為以下三類預(yù)處理方法：（1）基于數(shù)值統(tǒng)計規(guī)律的數(shù)據(jù)預(yù)處理方法，該類方法通過計算數(shù)據(jù)組的分布規(guī)律識別離散異常值[4-5]，但通常閾值設(shè)定需要人工干預(yù)，在延時數(shù)據(jù)預(yù)處理方面應(yīng)用較好[6-7]。（2）基于時空特征的數(shù)據(jù)預(yù)處理方法，該類方法將觀測要素的基于時間/空間的變化特征進行相關(guān)性分析，使得預(yù)處理結(jié)果符合時間/空間變化規(guī)律，主要包括卡爾曼濾波算法[8-10]、自回歸移動平均模型（ARMA）[11]、基于時序相關(guān)性分析方法[12]等。（3）基于人工智能的數(shù)據(jù)預(yù)處理方法，該類方法以大數(shù)據(jù)量的歷史數(shù)據(jù)挖掘分析為主，具有數(shù)據(jù)量大、算法復(fù)雜度高特點，適用于監(jiān)測數(shù)據(jù)的延時預(yù)處理過程[13-14]。綜上，基于數(shù)值統(tǒng)計規(guī)律的數(shù)據(jù)預(yù)處理方法在數(shù)據(jù)異常值剔除方面具有良好的表現(xiàn)，但不能實現(xiàn)浮標原始數(shù)據(jù)的實時去噪聲處理，而浮標數(shù)據(jù)采集器的計算能力無法滿足基于人工智能方法的高復(fù)雜度、大計算量的需求。基于時空特征的數(shù)據(jù)預(yù)處理方法計算開銷較低且符合浮標測量要素隨時間序列變化特征，但在工程應(yīng)用上，浮標端尚未實現(xiàn)原始數(shù)據(jù)的自動化預(yù)處理，國內(nèi)外的相關(guān)研究也以數(shù)據(jù)異常值檢測或延時數(shù)據(jù)的批量質(zhì)控處理為主，對于在浮標端開展原始測量數(shù)據(jù)實時自動化預(yù)處理方面的研究較少。因此，本研究結(jié)合海洋浮標的數(shù)據(jù)特征，提出一種基于改進的自適應(yīng)卡爾曼濾波自動化數(shù)據(jù)預(yù)處理方法，該方法采用箱形圖法對原始數(shù)據(jù)中的異常值進行檢測和修正，以期通過改進的卡爾曼濾波方法實現(xiàn)浮標觀測數(shù)據(jù)的噪聲濾波處理，并在廣東某海域進行實際的海洋工程應(yīng)用實驗，以驗證該方法的有效性和可行性。

1 數(shù)據(jù)預(yù)處理方法

1.1 浮標數(shù)據(jù)特征

在確定數(shù)據(jù)預(yù)處理方法前，對海洋浮標搭載的傳感器采樣原理、數(shù)據(jù)特征和采樣/融合頻率等進行分析，以保證數(shù)據(jù)處理方法的合理應(yīng)用，提高數(shù)據(jù)預(yù)處理效果。

依據(jù)GB/T 14914.3—2021《海洋觀測規(guī)范 第3部分：浮標潛標觀測》[15]和HY/T 147.6—2013《海洋監(jiān)測技術(shù)規(guī)程第6 部分：海洋水文、氣象與海冰》[16]等海洋環(huán)境監(jiān)測相關(guān)標準，各監(jiān)測要素在到達融合和回傳時間后，浮標端將監(jiān)測期間的原始數(shù)據(jù)進行提取并融合成該時段的測量值，浮標數(shù)據(jù)具有采樣與融合頻率固定、單參數(shù)融合、隨時間序列演變等特征。因此，必須在融合前對該時段內(nèi)所有監(jiān)測數(shù)據(jù)進行統(tǒng)一的去噪聲處理。同時，浮標監(jiān)測參數(shù)在時間序列上符合環(huán)境變化規(guī)律，在設(shè)計自動化數(shù)據(jù)處理程序時可以根據(jù)參數(shù)的融合頻率，采用基于預(yù)測與修正的實時濾波去噪方法。

1.2 卡爾曼濾波算法

卡爾曼濾波算法是一種具有較強容錯能力，以最小均方誤差為準則的遞推無偏估計算法，可有效去除系統(tǒng)中的噪聲干擾，在無人導(dǎo)航、電氣控制、傳感器數(shù)據(jù)融合等多個領(lǐng)域得到廣泛應(yīng)用[17-18]。傳統(tǒng)的卡爾曼濾波是在標準條件下獲得的，但海洋浮標是搭載多傳感器的實時環(huán)境監(jiān)測平臺，在工程應(yīng)用中由于噪聲統(tǒng)計特性不準確而產(chǎn)生濾波精度降低甚至產(chǎn)生濾波發(fā)散問題。為解決卡爾曼濾波發(fā)散問題，有學(xué)者提出Sage-Husa 自適應(yīng)卡爾曼濾波方法[19]，該方法通過時變噪聲估計器，在實時處理數(shù)據(jù)過程中，更新系統(tǒng)噪聲和測量噪聲，從而抑制濾波發(fā)散。其狀態(tài)空間模型如下：

式(1)中，xk為系統(tǒng)k時刻的n維狀態(tài)量；Φk,k-1為系統(tǒng)從k-1 時刻到k時刻的n×n維狀態(tài)轉(zhuǎn)移矩陣；zk為系統(tǒng)k時刻的m維觀測量；Hk為系統(tǒng)k時刻的m×n維測量矩陣；wk-1為系統(tǒng)k-1 時刻的過程噪聲；vk為m維測量噪聲，滿足如下統(tǒng)計特性：

式(2)中，j,k=1,2,3,...，當k=j時，δkj=1，否 則δkj=0；qk和qj為wk和wj的數(shù)學(xué)期望；rk為vk的數(shù)學(xué)期望；Qk和Rk為系統(tǒng)噪聲和測量噪聲方差矩陣。

Sage-Husa自適應(yīng)卡爾曼濾波算法可描述為

1.3 改進的自適應(yīng)卡爾曼濾波

受海洋環(huán)境噪聲、海洋生物附著和傳感器元件老化等因素的影響，浮標數(shù)據(jù)不可避免的存在異常測量值，這些異常值將導(dǎo)致算法在式(3)求狀態(tài)值時被放大并引入到狀態(tài)估計中，特別是當異常值偏離較大時，殘差無法正確修正，會嚴重影響濾波精度，因而在使用卡爾曼濾波器進行濾波時，必須對測量異常數(shù)據(jù)進行剔除并修正。同時，測量噪聲方差在計算過程中，當初始噪聲設(shè)置過大或?qū)嶋H海洋環(huán)境噪聲與理論建模值相比較小時，卡爾曼濾波算法將出現(xiàn)濾波發(fā)散問題。為提高算法的穩(wěn)定性，需要解決異常值對濾波的干擾和抑制濾波發(fā)散問題。

1.3.1 異常值修正 基于統(tǒng)計規(guī)律的方法在異常值檢測方面應(yīng)用最為廣泛，如拉依達準則和狄克遜準則等基于測量值正態(tài)分布規(guī)律的方法。但基于正態(tài)分布的疏失誤差檢測方法是以假定數(shù)據(jù)服從正態(tài)分布為前提的，其判斷異常值的標準是以計算數(shù)據(jù)組的均值和標準差為基礎(chǔ)的，異常值本身會對它們產(chǎn)生較大影響，進而導(dǎo)致方法的耐抗性較小。同時，海洋浮標數(shù)據(jù)往往難以服從正態(tài)分布規(guī)律，在使用基于正態(tài)分布的異常值檢測方法時存在漏剔除和誤剔除現(xiàn)象?；诮y(tǒng)計規(guī)律的還有箱形圖法，如圖1所示。

圖1 箱形圖結(jié)構(gòu)Fig.1 Box plot Structure

該方法利用數(shù)據(jù)序列中的5個統(tǒng)計量：下限、下四分位數(shù)（Q1）、中位數(shù)（Q2）、上四分位數(shù)（Q3）與上限來描述數(shù)據(jù)離散程度的方法。假設(shè)預(yù)處理數(shù)據(jù)集含n個測量值，則箱形圖Q1取1×(n+1)/4 位置的數(shù)；若n為奇數(shù)，則Q2取(n+1)/2 位置的數(shù)，否則取n/2 和n/ 2+1 位置數(shù)的均值；Q3取3 × (n+1)/ 4 位置的數(shù)。因此，四分位距IQR 為Q3—Q1，上限Q3+1.5 × IQR 和下限Q1—1.5 × IQR 作為異常值截斷點，測量值分布在該區(qū)間內(nèi)的判定為有效測量值，否則判定為異常測量值。箱形圖判斷異常值的標準以四分位數(shù)和四分位距為基礎(chǔ)，具有一定的耐抗性，能夠客觀地識別數(shù)據(jù)中的異常值。在識別到異常值后，為使得數(shù)據(jù)維數(shù)保持一致，本研究使用臨近數(shù)據(jù)求取的平均值進行插補。

1.3.2 方差受限抑制濾波發(fā)散 為避免初始噪聲設(shè)置過大或?qū)嶋H系統(tǒng)噪聲與理論建模值相比較小，將造成測量噪聲方差失去正定性而導(dǎo)致濾波發(fā)散問題。本研究設(shè)計采用方差受限的方法解決測量噪聲導(dǎo)致的濾波發(fā)散問題，這種改變雖然降低了一定的濾波精度，但提高了濾波的穩(wěn)定性。

式(4)中，βk作為方差加權(quán)系數(shù)，其取值為βk=βk-1/(βk-1+b)，其中β0=1，β∞=1 -b，b為漸消因子，取值區(qū)間為[0.90,0.99]；Rmin為設(shè)置的方差下限，保證為正，Rmax為設(shè)置的方差上限，能夠快速降低測量值的可信度并判定是否異常，若異常則放棄本次更新。

2 數(shù)據(jù)處理流程

結(jié)合海洋浮標數(shù)據(jù)特征和處理程序，提出基于改進擴展卡爾曼濾波的自動化數(shù)據(jù)處理方法，如圖2所示。

圖2 數(shù)據(jù)預(yù)處理流程Fig.2 Data preprocessing process

在到達參數(shù)數(shù)據(jù)融合/回傳時間后，浮標端將監(jiān)測時段內(nèi)的原始數(shù)據(jù)集提取作為數(shù)據(jù)預(yù)處理方法的輸入，對數(shù)據(jù)集依次進行常規(guī)檢驗、箱型圖異常值檢測以及卡爾曼濾波去噪處理。主要步驟如下：

（1）判斷是否到達參數(shù)數(shù)據(jù)融合與回傳時間。若到達則提取該監(jiān)測時段內(nèi)參數(shù)的所有采樣數(shù)據(jù)，作為預(yù)處理數(shù)據(jù)集D；若未到達則繼續(xù)進行采樣監(jiān)測。

（2）數(shù)據(jù)常規(guī)校驗。對數(shù)據(jù)集D內(nèi)的各測量值進行遍歷處理，判斷數(shù)據(jù)是否符合格式要求，對不符合要求的數(shù)據(jù)進行修正或剔除，同時判斷測量值是否在理論區(qū)間，對于超過傳感器理論量程的數(shù)據(jù)予以剔除，輸出數(shù)據(jù)集D1。

（3）異常值檢測。完成常規(guī)校驗后，以數(shù)據(jù)集D1作為箱形圖輸入，統(tǒng)計該時段內(nèi)的數(shù)據(jù)分布規(guī)律，對于異常值進行剔除并使用臨近數(shù)據(jù)的平均值進行插補，保證數(shù)據(jù)序列隨時間的變化特征，輸出數(shù)據(jù)集D2。

（4）初始化卡爾曼濾波算法的參數(shù)?？紤]到浮標監(jiān)測系統(tǒng)噪聲由傳感器元件自身的特性所決定的，如風(fēng)速采樣頻率為3 s，浮標觀測系統(tǒng)相對穩(wěn)定，相鄰間隔的測量值變化不大，因此根據(jù)k—1 時刻的最優(yōu)估計值來預(yù)測k時刻的預(yù)測值時，我們由經(jīng)驗認為，測量噪聲和系統(tǒng)噪聲近似為高斯噪聲，變量Φ和H初始設(shè)置為1，方差初始設(shè)置為q=0.05和r=0.50。

（5）去噪聲處理。采用改進的自適應(yīng)卡爾曼濾波方法進行去噪聲處理，以數(shù)據(jù)集D2各測量值按時間序列作為算法輸入，輸出預(yù)處理后的數(shù)據(jù)集D3。

3 應(yīng)用與分析

3.1 海洋工程應(yīng)用

2021年6月28日至2021年8月16日，本研究數(shù)據(jù)預(yù)處理方法的海洋浮標在廣東汕尾海域進行實際工程應(yīng)用，該浮標搭載2 臺風(fēng)速風(fēng)向儀、1 臺聲學(xué)多普勒流速剖面儀（ADCP）、1 臺波浪儀進行海洋風(fēng)、水溫、海流、波浪的監(jiān)測，實時回傳數(shù)據(jù)為廣東某地的系統(tǒng)工程建設(shè)提供數(shù)據(jù)支撐。站點位置與傳感器配置情況詳見表1。

表1 浮標系統(tǒng)部署與配置情況Table 1 Deployment and configuration of buoy system

3.2 應(yīng)用結(jié)果與分析

3.2.1 系統(tǒng)運行穩(wěn)定性 傳感器完成數(shù)據(jù)采集后實時發(fā)送至浮標端的數(shù)據(jù)采集與控制系統(tǒng)，由該系統(tǒng)按照數(shù)據(jù)融合/回傳時間，提取原始數(shù)據(jù)進行異常值修正和濾波降噪處理，并同步存儲在浮標本地。海上實驗結(jié)果表明，在整個監(jiān)測周期內(nèi)，部署本研究自動化數(shù)據(jù)預(yù)處理程序的浮標系統(tǒng)，數(shù)據(jù)采集率達100%，云數(shù)據(jù)中心數(shù)據(jù)接收率達97%。因此，本研究提出的自動化數(shù)據(jù)預(yù)處理程序計算開銷較低且不影響浮標系統(tǒng)的正常采集和數(shù)據(jù)融合工作，系統(tǒng)整體運行穩(wěn)定。

3.2.2 常規(guī)檢驗及異常值檢測情況 2021 年8 月16日回收浮標后，將原始數(shù)據(jù)和處理后的數(shù)據(jù)下載后進行統(tǒng)計分析。2021 年7 月18 日臺風(fēng)“查帕卡”期間浮標監(jiān)測的風(fēng)、海流流速數(shù)據(jù)異常值檢測前后對比情況如圖3所示。

圖3 數(shù)據(jù)異常值檢測Fig.3 Data outlier detection

圖3 顯示了數(shù)據(jù)異常值的檢測過程，由圖3（a、c、e）可知，風(fēng)、海流的原始數(shù)據(jù)中存在明顯的異常測量值，這是由于傳感器受海洋浮標搖晃、電信號不穩(wěn)定及環(huán)境噪聲干擾產(chǎn)生的。由于海洋浮標測量要素具有隨時間有序變化的特征，如圖3（b、d、f）所示，在到達數(shù)據(jù)融合時間后，通過常規(guī)檢驗和箱形圖法能夠很好地統(tǒng)計區(qū)間范圍內(nèi)的數(shù)據(jù)分布情況和識別離群異常值。2021 年6 月28 日至8 月16 日期間，浮標數(shù)據(jù)采集器未檢測到傳感器測量值超出理論閾值情況，但識別出部分數(shù)據(jù)存在格式錯誤和異常測量值。如表2 所示，浮標數(shù)據(jù)采集器檢測到各傳感器錯誤率均較低，最高的ADCP 錯誤率為8.90‰，主要原因是串口通信存在一定的通信誤碼率，而剖面海流的測量參數(shù)較多且單組數(shù)據(jù)量較大，單個參數(shù)格式有誤則直接判定為格式錯誤。浮標數(shù)據(jù)采集器對格式錯誤進行初步修正后存儲至本地，并對箱型圖檢測出的異常測量值使用臨近數(shù)據(jù)求取平均值進行插補。

表2 常規(guī)檢驗及異常值檢測情況Table 2 Routine inspection and detection of outliers

3.2.3 卡爾曼濾波數(shù)據(jù)預(yù)處理與應(yīng)用分析 在完成常規(guī)校驗和異常值檢測后，浮標端將在每個融合/數(shù)據(jù)回傳時刻，提取該時間段內(nèi)的原始數(shù)據(jù)并采用改進的自適應(yīng)卡爾曼濾波方法進行去噪聲處理。如圖4所示，浮標運行期間水溫、海流流速和風(fēng)速的原始數(shù)據(jù)經(jīng)過濾波處理前后的對比。

圖4 原始數(shù)據(jù)與卡爾曼濾波處理結(jié)果對比Fig.4 Comparison of raw data and Kalman filter processing results

根據(jù)圖4，原始數(shù)據(jù)經(jīng)卡爾曼濾波方法處理后，有效降低了水溫、海流流速和風(fēng)速的數(shù)據(jù)離散度，數(shù)據(jù)曲線更為平滑且符合隨時間序列變化的規(guī)律，濾波效果良好。由于卡爾曼濾波算法“預(yù)測+修正”的特征，原始數(shù)據(jù)經(jīng)濾波后，不僅能夠消除測量噪聲，也過濾了原始測量數(shù)據(jù)的高頻信號。但部分海上施工作業(yè)時甚至只關(guān)注極值、最大值和海流流速等高頻信號的變化，使得該方法在部分海洋工程領(lǐng)域存在應(yīng)用局限性。因此，各融合時段的平均風(fēng)速、平均風(fēng)向、水溫、氣壓、氣溫、平均波高等需要數(shù)據(jù)融合的參數(shù)，可直接采用濾波后的數(shù)據(jù)進行算數(shù)或矢量平均；極大風(fēng)、最大風(fēng)、最大波高等極值統(tǒng)計參數(shù)，則需要結(jié)合原始測量數(shù)據(jù)進行計算。針對大剖面、長序列的海流數(shù)據(jù)，研究人員一般更關(guān)注高頻信號的變化以開展中大尺度的海洋動力過程研究，如內(nèi)波實時監(jiān)測與預(yù)警，通過計算原始數(shù)據(jù)與濾波后數(shù)據(jù)的殘差，獲取剖面海流流速的高頻信號，進而通過海流高頻變化特征實現(xiàn)浮標端內(nèi)波的自動化識別[20-21]。因此，本研究提出的改進的自適應(yīng)卡爾曼方法在有效消除測量噪聲的同時，也能夠間接應(yīng)用于海流等參數(shù)高頻信號的提取，驗證了本研究方法的可行性。

3.3 與其他方法的仿真對比

為進一步驗證算法的準確性、可靠性和抑制濾波發(fā)散能力，將本研究算法與標準卡爾曼濾波方法[8]、自適應(yīng)卡爾曼濾波方法[9]和改進的Sage-Husa自適應(yīng)卡爾曼濾波方法[22]的數(shù)據(jù)處理效果進行對比分析，運行硬件為具有一定計算能力的低功耗數(shù)據(jù)采集設(shè)備（Raspberry Pi，運行內(nèi)存2 GB）。如圖5（a、b）所示，該實驗數(shù)據(jù)采用已知正弦信號作為模擬的實測值，添加高斯噪聲（sigma=0.5）和突變數(shù)據(jù)的正弦信號作為傳感器的測量值作為仿真實驗的原始數(shù)據(jù)，因為真實值已知，通過將預(yù)處理后的數(shù)據(jù)與真實測量值進行誤差比對分析，能夠很好地評價數(shù)據(jù)預(yù)處理效果和驗證算法的穩(wěn)定性。3 種算法的數(shù)據(jù)預(yù)處理效果對比如圖5（c）所示。

圖5 不同預(yù)處理方法仿真對比Fig.5 Simulation comparison of different preprocessing methods

通過計算處理后數(shù)據(jù)與正弦信號真實值之間的均方根誤差（RMSE）、平均絕對誤差（MAE）、相關(guān)系數(shù)（R）及算法運行耗時等指標，對比分析不同方法的處理效果。如表3 所示，文獻[8]方法受到數(shù)據(jù)噪聲及突變數(shù)據(jù)的影響，即噪聲過大及噪聲與理論建模值相比較小時，隨著誤差累積出現(xiàn)了濾波發(fā)散現(xiàn)象。文獻[9]、文獻[22]及本研究方法能夠有效的降低異常值對卡爾曼濾波發(fā)散的影響，濾波處理后的數(shù)據(jù)與真實數(shù)據(jù)曲線擬合效果較好，且本文方法處理后的數(shù)據(jù)均方根誤差、平均絕對誤差最低，與原始數(shù)據(jù)集的相關(guān)系數(shù)達0.98，數(shù)據(jù)預(yù)處理結(jié)果的準確性最優(yōu)。另外，本研究方法利用箱型圖法處理突變數(shù)據(jù)和限定方差上限，在有效降低異常值對濾波精度影響的同時，能夠快速降低測量值的可信度，降低了算法的計算開銷，其計算耗時相比文獻[9]和文獻[22]降低了50%，計算耗時最少。

表3 對添加噪聲及突變數(shù)據(jù)的預(yù)處理效果評價對比Table 3 Comparison of preprocessing performance evaluation on adding noise and mutation data

4 結(jié)論

結(jié)合海洋浮標數(shù)據(jù)采集特征，本文提出了基于改進自適應(yīng)卡爾曼濾波的自動化數(shù)據(jù)預(yù)處理方法，主要結(jié)論如下：

（1）該方法通過箱形圖進行異常測量值的檢測，并采用方差受限的方法解決測量噪聲導(dǎo)致的濾波發(fā)散問題。

（2）海上工程應(yīng)用結(jié)果表明，浮標端常規(guī)數(shù)據(jù)采集器能夠滿足該方法的計算開銷，數(shù)據(jù)采集及自動化數(shù)據(jù)預(yù)處理程序運行穩(wěn)定，浮標系統(tǒng)有效數(shù)據(jù)接收率達97%。根據(jù)統(tǒng)計，常規(guī)檢驗和箱形圖法能夠有效地識別格式錯誤和離群異常值，且各時段測量數(shù)據(jù)經(jīng)過異常值修正和濾波降噪處理后，數(shù)據(jù)曲線更為平滑且符合隨時間序列變化的規(guī)律，其濾波處理后的數(shù)據(jù)可應(yīng)用于各測量要素的數(shù)據(jù)融合與高頻信號的間接提取。

（3）與同類方法的仿真對比發(fā)現(xiàn)，本研究方法能夠有效的降低異常值對卡爾曼濾波發(fā)散的影響，濾波處理后的數(shù)據(jù)與真實數(shù)據(jù)曲線擬合效果較好，且計算耗時最小。

本研究成果豐富了我國海洋浮標數(shù)據(jù)預(yù)處理與質(zhì)控手段，可為其他領(lǐng)域的基于時間序列變化的數(shù)據(jù)處理提供參考和借鑒。