基于趨勢項距平值迭進式方差控制的SBE9plus CTD奇異值剔除方法

于 龍 ，孫 佳，熊學軍*，郭延良

(1.自然資源部 第一海洋研究所,山東 青島 266061；2.自然資源部 海洋環境科學與數值模擬重點實驗室,山東 青島 266061；3.山東省海洋環境科學與數值模擬重點實驗室,山東 青島 266061；4.青島海洋科學與技術試點國家實驗室 區域海洋動力學與數值模擬功能實驗室,山東 青島 266237)

現場調查是認知和把握海洋環境最直接、最有效的方法,為了能準確地揭示海洋環境的時空分布特征和變化規律,保證所獲取調查數據的準確性和可靠性就顯得尤為重要。目前觀測方式和觀測儀器多種多樣,數據處理日趨自動化。然而,在調查現場,由于儀器或人為因素所產生的測量誤差在所難免,主要包括隨機誤差和過失誤差[1-3]。這就要求數據的后期處理人員對數據獲取情況有足夠了解,包括數據來源、現場調查環境等,在進行數據處理之前,根據現場情況對數據進行質量審查和控制,做出切合實際的判斷,保留下那些能反應海洋環境異常變化的數據,去掉因誤差造成的奇異值,盡最大可能還原現場海洋環境。

目前,對SBE9plus CTD(美國SEA-BIRD 公司生產的溫、鹽、深測量系統)剖面觀測資料的處理,基本按照儀器自帶的SBE Data Processing處理軟件的流程進行[4-6],主要步驟包括低通濾波(Filter)、滯后訂正(Align CTD)、熱效應訂正(Cell Thermal Mass)、消除逆壓(Loop Edit)、奇異值剔除(Wild Edit與Window Filter)及衍生計算(Derive)等。其中,低通濾波對壓力數據進行平滑,為后續消除逆壓數據做準備；滯后訂正目前多通過SBE11甲板單元實現[2,4],并在原始數據頭文件中給出訂正說明；對于熱效應訂正,可以采用Lueck[7]和Morison等[8]提供的校正系數對電導率進行熱效應校正。通過低通濾波-滯后訂正-熱效應訂正-消除逆壓這4個步驟,可有效去除因各傳感器不同步造成的數據“尖峰”及船舶起伏造成的數據“打結”。但在進行奇異值剔除時筆者發現,SBE Data Processing軟件中的Wild Edit與Window Filter程序需要數據處理人員根據經驗進行多次嘗試[4,9],處理結果受人為因素影響較大,或條件過松無法盡除奇異值,或條件過嚴造成正常數據失真。此外,亦可采用中位數濾波結合最大偏差(Max Difference,MDIF)[10]或中位數絕對偏差(Mean Absolute Deviation,MAD)[11]進行奇異值剔除。其中,中位數濾波結合最大偏差法剔除效果較好,但其剔除結果取決于處理者對MDIF或最大梯度選擇的合適程度,而MDIF隨調查海區的不同而變化。若無法有效剔除溫度、電導率等直接觀測要素中的奇異值,則鹽度、密度等導出量的數據質量更無從談起。為此,本文對溫度、電導率等直接觀測量的奇異值剔除進行討論。

首先明確奇異值的定義。每一個物理量均客觀存在,其測量值分為真值、誤差限內值和超出誤差限值。絕對真實值無法通過測量得到,只能通過提高測量儀器的測量精度、減少環境因素和人為因素所造成的誤差等,將測量值控制在誤差限內,進而無限逼近其絕對真實值；而超出誤差限值的數據則被認為無效數據,即奇異值。其特征如下:①物理上違背了觀測值分布的基本規律和趨勢,超出了觀測對象的物理性范疇；②由于海水連續性特征,溫度、鹽度和密度等大多數物理海洋要素在時間與空間分布上的變化趨勢,表現在曲線上應基本呈平滑狀態,而奇異值在數學上違背了數據連續性法則,表現在曲線上即與曲線平滑性顯著沖突；③儀器性能上表現在傳感器的物理偏差、反應時間滯后及各傳感器之間未能達到同步性調整等,從而造成“尖峰”等突變現象。

本文以SBE9plus CTD 剖面數據為例,提出了一種基于趨勢項距平值迭進式方差控制的奇異值剔除方法。此方法主要參數涉及觀測儀器采樣頻率與下放速度,可為其他型號CTD 剖面數據的奇異值剔除提供參考。

1 資料與方法

1.1 資料介紹

原始數據是以船載觀測方式觀測得到的多個站點的溫、鹽、深單點剖面資料,其中實例數據D1如圖1所示。觀測儀器為SBE9plus CTD 溫鹽深測量系統,采樣頻率24 Hz,其溫度傳感器分辨率為0.000 2 ℃,準確度為0.001 ℃；電導率傳感器分辨率為0.000 4 mS/cm,準確度為0.003 mS/cm；壓力傳感器分辨率為0.068 95 dbar(0.001%FS),準確度為1.034 dbar(0.015%FS)。由于SBE9plus CTD 溫、鹽、深測量系統的壓力傳感器輸出單位為dbar,為方便表述以及與儀器自帶軟件處理結果進行比較,我們在后續分析中仍然用dbar作為壓力單位,而不使用壓力的國際標準單位Pa,其換算關系為1 dabr＝104Pa。

如圖1所示,原始剖面數據中溫度和電導率均存在逆壓數據,經過消除逆壓(Loop Edit)后,僅電導率存在明顯奇異值(圖2)。觀測甲板單元在現場觀測時已對電導率傳感器進行了滯后訂正(Align CTD)0.073 s,因此在數據處理過程中對僅對電導率進行熱效應訂正(Cell Thermal Mass)。但在采用Wild Edit與Window Filter程序進行奇異值剔除時,需經多次嘗試,且最終剔除結果或因條件過松無法盡除奇異值,或因條件過嚴造成正常數據失真。

圖1 CTD 原始數據D1Fig.1 The CTD raw data at D1

圖2 消除逆壓后CTD 數據Fig.2 The CTD data after Loop Edit

1.2 研究方法

由于原始數據在整個剖面上變化范圍較大,特別是在溫鹽數據本身的變化幅度與奇異值導致的變化幅度之間存在交叉的躍層區或逆躍層區等海洋環境變化劇烈的水層,因此很難直接提出有效的奇異值判別條件,但奇異值具有震蕩性,而躍層區與逆躍層區正常數據具有單調性。

針對以上區別,首先通過EMD(Empirical Mode Decomposition)方法[12]求取原始數據的最優趨勢項,進而與原始數據做差求得距平,以距平為切入點進行奇異值的判別與剔除。具體步驟如下:

1)求取最優趨勢項與距平

在求取趨勢項時本文采用EMD 方法[12],該方法對非線性、非平穩過程數據進行距平化效果顯著。通過EMD 方法得到的殘余項本身就是趨勢項,并且,還可以根據物理背景的不同,調整殘余項與部分模態分量相加作為趨勢項[13]。為了使趨勢項盡可能地反映原始數據的所有信息,本文以王金良和李宗軍提出的最小方差比率[14](后n項模態分量和的方差與原始數據方差的比)作為最優趨勢項的判別條件。

2)識別奇異值區

根據數據采樣頻率F和觀測儀器移動速度V,確定每k＝min(F/V)個數據作為一個區塊,對距平A進行分割,得到一系列距平區塊A n(j)(其中j＝1,2,3,…,k,代表每個距平區塊內的數據量,max(j)＝k；n為距平區塊個數,n＝1,2,3,…,max(n)＝[max(i)/k],i為原始數據的下標),若無法整分,取A中的最后k個數據作為最后一個距平區塊,保證每個區塊包含k個數據；然后求出每個距平區塊的方差V n。由于這一步僅為尋找可能存在奇異值的區域,因此可以縮小方差控制區間,放寬判別條件,即對距平區塊方差通過1倍方差控制識別奇異值區。本文所用資料的觀測儀器采樣頻率為24 Hz,觀測儀器下放速度在1.0～2.0 m/s,故取k＝12,使得區塊的深度變化范圍在1 m 之內。對于1 Hz等低頻采樣頻率型號的CTD,可嘗試直接進行數學分割或采用其他方法,但識別效果有待進一步研究。

3)剔除奇異值

奇異值區內的數據并不一定全部為奇異值,因此并不對上一步識別出的奇異值區進行整體方差控制,而是針對識別出的每個奇異值區Qm進行逐點判別。在判別第一個奇異值區內的第一個數據點時,由于其前面的數據可以確定為非奇異值,其后面的數據則暫時無法判定,因此取與其相鄰的至少前2k－1個正常數據共同組成至少2k個數據量的被判別區,保證被判別區內的數據量至少為每個區塊數據量的2倍。由于這一步需要識別出真正的奇異值對其進行剔除,從而擴大方差控制區間,加嚴判別條件,謹防誤判,因此進行3倍方差控制。針對識別出的奇異值,在保證上升過程數據趨勢與下降過程數據趨勢基本一致且對上升過程數據質量有把控的前提下,取上升過程對應水層正常數據予以平移替代。在第一個數據點被糾正后可作為下一個被判別點的正常引入數據點,依次逐點推進,最終剔除整個剖面數據奇異值。

2 方法應用

以數據D1(圖1)進行方法說明,在應用本方法之前對原始數據進行初步質量控制,包括去除感溫數據、逆壓數據及明顯不符合常識的數據等,上述初步質量控制步驟完全可以通過編程實現自動化處理[15]。

第一步,采用EMD 方法求取最優趨勢項。如圖3a所示,當模態分量個數為9時,方差比率最小。即后9項模態分量的和為最優趨勢項(圖3b中紅色實線)。將最優趨勢項與原始數據做差得到距平(圖3c)。

圖3 方差比率、最優趨勢項及距平Fig.3 Variance ratio,optimal trend and anomalies

第二步,識別奇異值區。將第一步得到的距平按照每12個數據一個區塊進行分割,并求取每個區塊的方差,然后對距平區塊方差通過1倍方差控制識別奇異值區。數據的4個奇異值區如圖4所示,依次與原始數據(圖1紅線)中從上到下的4個“尖峰”相對應。可以看到,原始數據中的4個奇異值區全部被識別。

第三步,剔除奇異值。針對識別出的奇異值(圖4黑色圓圈),在保證上升過程數據趨勢與下降過程數據趨勢基本一致,且對上升過程數據質量有把控的前提下,取上升過程對應水層正常數據予以平移替代(圖4紅色實點)。在第一個數據點被糾正后可作為下一個被判別點的正常引入數據點,依次逐點推進,最終剔除整個剖面數據奇異值。最終結果如圖5所示,在不改變其他正常數據的前提下,對奇異值進行了精準識別與靶向剔除。

圖4 奇異值區、奇異值及使用同水層上升數據平移替代結果Fig.4 The identification results of abnormal block,abnormal data and their replacement with upcast data

圖5 電導率原始數據、奇異值剔除結果、趨勢項距平值及距平區塊方差Fig.5 The raw data,elimination result,trend anomaly and variances of blocks

3 方法對比

在數據D1(圖1)的基礎上,增加另外6個同樣存在奇異值的CTD 剖面數據實例(圖6),分別使用本文提出的距平迭進方差法、SBE Data Processing中的Wild Edit與Window Filter方法以及中位數濾波結合最大偏差法共四種方法進行處理和對比分析,結果如下:

①采用本文提出的距平迭進方差法,處理過程無需針對每個實例進行個性參數設置,可對所有實例進行批處理,均可對所有實例中的奇異值進行靶向剔除,原始數據如圖6黑點及圖7a所示,剔除結果如圖6紅點及圖7b所示。

圖6 CTD 原始數據(黑點)及本方法剔除結果(紅點)Fig.6 The CTD raw data(black dot)and eliminate results(red dot)

②利用SBE Data Processing處理軟件,首先對原始數據進行初步質量控制,去除感溫數據及明顯異常數據,然后依次執行Filter-Cell Thermal Mass-Loop Edit,采用2倍方差與4倍方差控制,從100～500不斷改變區塊大小多次運行Wild Edit[4]。但處理后數據D3～D7五個實例中仍存在明顯奇異值,且繼續多次改變區塊大小運行Wild Edit后,剔除效果再無增益,數據D1剔除結果如圖7c所示。

③執行SBE Data Processing-Window Filter-Median程序,剔除效果隨濾波窗口大小的不同而存在差異,對于數據D1,與中位數濾波結合MDIF方法窗口保持一致取25時,無法剔除奇異值,如圖7d所示。若窗口取為60,剔除效果稍有改善。此方法對其他正常數據同樣進行了濾波處理,對正常數據改變較大。

④采用中位數濾波結合MDIF進行奇異值剔除時,首先計算7個實例原始數據的偏差情況,對整體偏差進行評估。由于7個CTD 剖面資料所處海域不同,最大偏差的選取也有所差異。對于數據D1,選取最大偏差MDIF＝0.05 mS·cm－1(圖8中灰色實線所示)進行奇異值剔除,結果如圖7e所示。

圖7 奇異值剔除方法比較Fig.7 Comparison of the methods of eliminating abnormal data

圖8 電導率偏差Fig.8 Conductivity differences

為了進一步比較4種方法對正常數據的影響,將原始數據中的正常數據與結果中的相應數據做差,結果如表1所示,距平跌進方差法對所有實例的正常數據均無改變；中位數濾波結合MDIF法對數據D1、數據D2、數據D5、數據D7 無改變,但對數據D3、數據D4、數據D6 正常數據的改變較大,最大改變量達0.33 mS·cm－1；SBE Data Processing處理軟件中的Wild Edit程序和Window Filter程序對所有實例的原始數據均有改變,最大改變量同樣達到0.33 mS·cm－1。

表1 正常數據改變量統計(mS·cm-1)Table 1 Statistics of the change of normal data(mS·cm－1)

綜上所述,SBE Data Drocessing處理軟件中的Wild Edit程序與Window Filter程序在處理過程中需要多次嘗試,處理效果具有不確定性,且對其他正常數據的改變較大；中位數濾波結合MDIF法剔除效果優于前面2種方法,但同樣存在改變其他正常數據的情況,且需根據不同剖面數據選取合適的MDIF；本文提出的距平迭進方差法可實現對異常值的精確識別與靶向剔除,避免多次嘗試、誤判以及處理結果的不確定性等,對高頻采樣頻率的CTD 剖面數據的異常值剔除具有普適性。

4 結語

本文針對SBE9plus CTD 剖面數據的奇異值剔除問題,基于多個CTD 剖面數據實例,提出了一種利用趨勢項距平值迭進式方差識別CTD 剖面數據奇異值的方法。首先通過EMD 方法求取最優趨勢項進而得到距平,通過方差控制識別出奇異值區,在奇異值區內逐點推進,再次通過方差控制識別奇異值。并與SBE Data Processing軟件中的Wild Edit程序和Window Filter程序,以及中位數濾波結合最大偏差法進行了比較分析。實例應用與方法對比結果表明,本文提出的基于趨勢項距平值迭進式方差控制的奇異值剔除方法具有以下特點:

①可以對奇異值進行精準識別與靶向剔除,對奇異值的剔除更有針對性,且不會改變原始數據中的正常數據,最大程度保留數據原始特性；

②可以一次性盡除奇異值,避免不斷改變判別閾值進行多次嘗試,避免因人為因素選擇不同的判別閾值而導致不同的剔除結果；

③該方法適用于其他高頻采樣頻率的CTD 剖面數據,對于1 Hz等低頻采樣頻率的CTD 剖面數據,其適用性有待進一步探討。