于 龍 ，孫 佳，熊學(xué)軍*，郭延良

(1.自然資源部 第一海洋研究所,山東 青島 266061；2.自然資源部 海洋環(huán)境科學(xué)與數(shù)值模擬重點(diǎn)實(shí)驗(yàn)室,山東 青島 266061；3.山東省海洋環(huán)境科學(xué)與數(shù)值模擬重點(diǎn)實(shí)驗(yàn)室,山東 青島 266061；4.青島海洋科學(xué)與技術(shù)試點(diǎn)國(guó)家實(shí)驗(yàn)室 區(qū)域海洋動(dòng)力學(xué)與數(shù)值模擬功能實(shí)驗(yàn)室,山東 青島 266237)

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

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

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

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

1 資料與方法

1.1 資料介紹

原始數(shù)據(jù)是以船載觀測(cè)方式觀測(cè)得到的多個(gè)站點(diǎn)的溫、鹽、深單點(diǎn)剖面資料,其中實(shí)例數(shù)據(jù)D1如圖1所示。觀測(cè)儀器為SBE9plus CTD 溫鹽深測(cè)量系統(tǒng),采樣頻率24 Hz,其溫度傳感器分辨率為0.000 2 ℃,準(zhǔn)確度為0.001 ℃；電導(dǎo)率傳感器分辨率為0.000 4 mS/cm,準(zhǔn)確度為0.003 mS/cm；壓力傳感器分辨率為0.068 95 dbar(0.001%FS),準(zhǔn)確度為1.034 dbar(0.015%FS)。由于SBE9plus CTD 溫、鹽、深測(cè)量系統(tǒng)的壓力傳感器輸出單位為dbar,為方便表述以及與儀器自帶軟件處理結(jié)果進(jìn)行比較,我們?cè)诤罄m(xù)分析中仍然用dbar作為壓力單位,而不使用壓力的國(guó)際標(biāo)準(zhǔn)單位Pa,其換算關(guān)系為1 dabr＝104Pa。

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

圖1 CTD 原始數(shù)據(jù)D1Fig.1 The CTD raw data at D1

圖2 消除逆壓后CTD 數(shù)據(jù)Fig.2 The CTD data after Loop Edit

1.2 研究方法

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

針對(duì)以上區(qū)別,首先通過(guò)EMD(Empirical Mode Decomposition)方法[12]求取原始數(shù)據(jù)的最優(yōu)趨勢(shì)項(xiàng),進(jìn)而與原始數(shù)據(jù)做差求得距平,以距平為切入點(diǎn)進(jìn)行奇異值的判別與剔除。具體步驟如下:

1)求取最優(yōu)趨勢(shì)項(xiàng)與距平

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

2)識(shí)別奇異值區(qū)

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

3)剔除奇異值

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

2 方法應(yīng)用

以數(shù)據(jù)D1(圖1)進(jìn)行方法說(shuō)明,在應(yīng)用本方法之前對(duì)原始數(shù)據(jù)進(jìn)行初步質(zhì)量控制,包括去除感溫?cái)?shù)據(jù)、逆壓數(shù)據(jù)及明顯不符合常識(shí)的數(shù)據(jù)等,上述初步質(zhì)量控制步驟完全可以通過(guò)編程實(shí)現(xiàn)自動(dòng)化處理[15]。

第一步,采用EMD 方法求取最優(yōu)趨勢(shì)項(xiàng)。如圖3a所示,當(dāng)模態(tài)分量個(gè)數(shù)為9時(shí),方差比率最小。即后9項(xiàng)模態(tài)分量的和為最優(yōu)趨勢(shì)項(xiàng)(圖3b中紅色實(shí)線)。將最優(yōu)趨勢(shì)項(xiàng)與原始數(shù)據(jù)做差得到距平(圖3c)。

圖3 方差比率、最優(yōu)趨勢(shì)項(xiàng)及距平Fig.3 Variance ratio,optimal trend and anomalies

第二步,識(shí)別奇異值區(qū)。將第一步得到的距平按照每12個(gè)數(shù)據(jù)一個(gè)區(qū)塊進(jìn)行分割,并求取每個(gè)區(qū)塊的方差,然后對(duì)距平區(qū)塊方差通過(guò)1倍方差控制識(shí)別奇異值區(qū)。數(shù)據(jù)的4個(gè)奇異值區(qū)如圖4所示,依次與原始數(shù)據(jù)(圖1紅線)中從上到下的4個(gè)“尖峰”相對(duì)應(yīng)?？梢钥吹?原始數(shù)據(jù)中的4個(gè)奇異值區(qū)全部被識(shí)別。

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

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

圖5 電導(dǎo)率原始數(shù)據(jù)、奇異值剔除結(jié)果、趨勢(shì)項(xiàng)距平值及距平區(qū)塊方差Fig.5 The raw data,elimination result,trend anomaly and variances of blocks

3 方法對(duì)比

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

①采用本文提出的距平迭進(jìn)方差法,處理過(guò)程無(wú)需針對(duì)每個(gè)實(shí)例進(jìn)行個(gè)性參數(shù)設(shè)置,可對(duì)所有實(shí)例進(jìn)行批處理,均可對(duì)所有實(shí)例中的奇異值進(jìn)行靶向剔除,原始數(shù)據(jù)如圖6黑點(diǎn)及圖7a所示,剔除結(jié)果如圖6紅點(diǎn)及圖7b所示。

圖6 CTD 原始數(shù)據(jù)(黑點(diǎn))及本方法剔除結(jié)果(紅點(diǎn))Fig.6 The CTD raw data(black dot)and eliminate results(red dot)

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

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

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

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

圖8 電導(dǎo)率偏差Fig.8 Conductivity differences

為了進(jìn)一步比較4種方法對(duì)正常數(shù)據(jù)的影響,將原始數(shù)據(jù)中的正常數(shù)據(jù)與結(jié)果中的相應(yīng)數(shù)據(jù)做差,結(jié)果如表1所示,距平跌進(jìn)方差法對(duì)所有實(shí)例的正常數(shù)據(jù)均無(wú)改變；中位數(shù)濾波結(jié)合MDIF法對(duì)數(shù)據(jù)D1、數(shù)據(jù)D2、數(shù)據(jù)D5、數(shù)據(jù)D7 無(wú)改變,但對(duì)數(shù)據(jù)D3、數(shù)據(jù)D4、數(shù)據(jù)D6 正常數(shù)據(jù)的改變較大,最大改變量達(dá)0.33 mS·cm－1；SBE Data Processing處理軟件中的Wild Edit程序和Window Filter程序?qū)λ袑?shí)例的原始數(shù)據(jù)均有改變,最大改變量同樣達(dá)到0.33 mS·cm－1。

表1 正常數(shù)據(jù)改變量統(tǒng)計(jì)(mS·cm-1)Table 1 Statistics of the change of normal data(mS·cm－1)

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

4 結(jié)語(yǔ)

本文針對(duì)SBE9plus CTD 剖面數(shù)據(jù)的奇異值剔除問(wèn)題,基于多個(gè)CTD 剖面數(shù)據(jù)實(shí)例,提出了一種利用趨勢(shì)項(xiàng)距平值迭進(jìn)式方差識(shí)別CTD 剖面數(shù)據(jù)奇異值的方法。首先通過(guò)EMD 方法求取最優(yōu)趨勢(shì)項(xiàng)進(jìn)而得到距平,通過(guò)方差控制識(shí)別出奇異值區(qū),在奇異值區(qū)內(nèi)逐點(diǎn)推進(jìn),再次通過(guò)方差控制識(shí)別奇異值。并與SBE Data Processing軟件中的Wild Edit程序和Window Filter程序,以及中位數(shù)濾波結(jié)合最大偏差法進(jìn)行了比較分析。實(shí)例應(yīng)用與方法對(duì)比結(jié)果表明,本文提出的基于趨勢(shì)項(xiàng)距平值迭進(jìn)式方差控制的奇異值剔除方法具有以下特點(diǎn):

①可以對(duì)奇異值進(jìn)行精準(zhǔn)識(shí)別與靶向剔除,對(duì)奇異值的剔除更有針對(duì)性,且不會(huì)改變?cè)紨?shù)據(jù)中的正常數(shù)據(jù),最大程度保留數(shù)據(jù)原始特性；

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

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