基于時序相關性分析方法的浮標異常數據識別

張宇 ，周燕，陶邦一,3*，顧吉星，趙傳高，郝增周，張藝蔚,5，黃海清，毛志華

(1.自然資源部第二海洋研究所 衛星海洋環境動力學國家重點實驗室，浙江 杭州 310012；2.浙江省海洋科學院，浙江 杭州 310007；3.南方海洋科學與工程廣東實驗室，廣東 廣州 511458；4.國家海洋局煙臺海洋環境監測中心站，山東 煙臺 264006；5.中國科學院上海技術物理研究所，上海 200083)

1 引言

海上浮標是獲得長時序、高精度海洋環境參數最主要的手段，確保浮標數據的質量可靠性是開展數據應用所面臨的首要問題，因此開展浮標異常數據的檢測識別是其中一項重要工作[1]。異常數據一般指超過正常合理數值范圍的以及偏離由海洋環境引起的變化規律的數據[2]。將臺州大陳（TZ01）和溫州南麂島（NJ01）的兩個浮標葉綠素數據與Aqua/MODIS、VIIRS和GOCI 海洋水色衛星反演的葉綠素產品進行比對（圖1），研究發現浮標反演的葉綠素產品存在兩種異常類型：（1）浮標數據在時序分布上連續且與衛星數據有較好的一致性（圖1a），但由于海洋隨機過程產生如圖1a 中紅色方框標記的跳變數據，屬于傳統意義上的跳變異常數據類型；（2）紅色條帶標記的一段浮標葉綠素測量異常數據呈現出：在時序變化過程中連續平穩，但隨時間逐漸推移，最后整體偏離正常數據的分布特征（圖1b），這種異常數據屬于一種新的漸變異常數據類型。基于海洋環境要素時序數據分布平穩的假設[3-6]，傳統的異常數據統計識別方法僅對跳變異常數據類型的數據檢測有效，而對漸變的質量異常數據類型無法識別[5-12]。主要原因在于異常發生的初始階段，其變化特征與由海洋環境變化引起的變化趨勢很難在沒有先驗知識的條件下進行區分，只有利用后續正常數據分布特征等后驗證知識進行識別。這類漸變異常數據可能與傳感器探頭受污、供電異常等因素有關。漸變的異常數據類型在長時間觀測的浮標數據中時有發生，因此如何在這一類型異常數據發生的初始階段進行有效識別，對于浮標異常的實時監測、及時維護、保證數據的可靠和連續性具有現實意義。

圖1 浮標數據與衛星數據葉綠素a 濃度對比Fig.1 Comparison of chlorophyll a concentration between buoy data and satellite data

國內外都已開展海上浮標觀測應用工作多年，但實際上實現各類型異常數據的自動檢測識別仍有較大難度[3-4]，國內對海洋數據的檢測主要依賴專家經驗、歷史資料以及常識形成的海洋環境監測數據檢驗標準庫[5]。目前已有的異常檢測方法主要有極值檢驗、一致性判斷、遞增性判斷、格拉布斯檢驗、狄克遜檢驗、拉依達檢驗、過度梯度檢測、尖峰檢測和無梯度檢測等[1,6-19]，這些異常檢測方法主要是針對單一參數在某一時間尺度的平穩隨機過程中進行統計學的分析處理，在傳統跳變異常數據類型識別中取得了較好的檢測效果，但對于漸變異常數據類型的自動檢測識別研究較少。

隨著技術發展，目前浮標平臺上搭載的傳感器數目和測量的參數越來越多，而在這些測量參數中存在某些參數相互關聯的特征。多元時間序列數據分析方法（如建立矢量自回歸（VAR）、多元譜分析，廣義自回歸條件異方差模型（GARCH）等）被廣泛地應用到質量異常數據的檢測和識別上[20-23]。Tsay[24]將4 種類型單參數時間序列異常數據識別方法拓展到了多元序列數據。此外，異常質量數據檢測中也應用到了矢量自相關性系數、差分整合移動平均自回歸模型、遺傳算法等方法[25-29]。然而，上述方法主要適用于跳變異常數據的識別，而且對數據平穩性要求較高、計算流程復雜，并未在本文發現的漸變異常數據類型上得到應用。其中，在海洋多元長時序數據異常識別方面，Schuckmann 等[12]提出相關性分析方法成功地識別了葉綠素濃度高而溶解氧濃度低的錯誤數據類型。在浮標數據質量控制中的應用，僅給出了白天葉綠素濃度高而溶解氧濃度低的錯誤數據類型的識別案例。竇文潔等[18]則根據海洋碳酸鹽系統中海水CO2分壓本身于水體溫度鹽度存在定量相關性關系的特點，在假設觀測參數變化在非常小的時間尺度內為一平穩過程的基礎上，提出了基于多參數觀測序列差分統計特征的異常點識別方法。雖然該方法由于僅基于參數平穩性假設而無法進一步有效識別漸變異常數據類型，但相比于單一參數分析方法，利用多參數強關聯性對異常數據進行檢測，會對數據的處理有更加全面、深入的把控。

目前，我國常規生態浮標通常會同時觀測酸堿度（pH）、溶解氧（DO）濃度以及葉綠素a（Chla）濃度等數據。大量的研究表明，它們之間雖然具有較緊密的相關性[30-31]，特別是在海水藻類生長暴發期間[31]，但它們并不存在穩定的相關關系，如謝群等[32]在流沙灣得出溶解氧濃度與葉綠素a濃度成正比例關系，尤其是冬季，海水中的溶解氧濃度與葉綠素a濃度具有極顯著正相關，春季次之，夏秋兩季兩者之間不存在相關性的結論。可見在不同海域、不同季節及不同海洋過程中的參數之間相關性特征具有明顯的差異性，并不能類似于Schuckmann 等[12]采用事先設定的相關性特征進行多年長時序數據的處理。Hollinger 和Richardson[33]在海洋數據不確定分析時提出了“單塔日變化法”，其基本假設是相鄰日期間在相同或相似的環境條件下數據變化過程相似。因此本研究認為，浮標觀測到的正常多參數數據不僅單一參數在時序變化上平穩連續，并且兩兩參數間的相關性在一定時序上穩定甚至一致。

本文基于上述假設，希望通過對浙江沿岸海域浮標多年的pH、DO 濃度、Chla濃度數據相關性進行分析，了解當某一參數出現漸變異常時，與其他參數的相關性特征的變化規律，基于多參數相關性變化提出一種簡單、適用的漸變異常數據檢測識別方法，并且探討該方法在該海域的適用性。

2 數據來源

浮標數據采用浙江省沿岸的溫州南麂大沙岙（NJ01）、臺州大陳（TZ01）、寧波南韭山（NB01）、寧波漁山（NB03）、舟山嵊泗綠華（ZS03）和舟山普陀東極（ZS04）6 處浮標數據，浮標分布如圖2 所示。

觀測時間在2012 年8 月至2017 年5 月之間，數據每15 min 或1 h 傳輸1 次，以同一浮標同一時間獲取的數據為一組，共計183 967 組數據，其中狀態顯示異常、故障或維護的數據有8 662 組，儀器運行正常狀態的原始浮標數據（DO 濃度、Chla濃度和pH）有175 305 組，占總數據量的95%，數據情況如表1 所示。對儀器運行正常狀態的175 305 組數據進行分析處理，對其他狀態的數據不予處理。

3 基于相關性的漸變異常數據識別方法

本文根據pH 與DO 濃度具有正相關關系，Chla濃度與pH 和DO 濃度的關系因藻類生長、季節變換等因素呈現顯著正相關或不相關關系等特點，利用最基本相關性統計學方法來計算pH 與DO 濃度、pH與Chla濃度、DO 與Chla濃度兩兩相關性系數。在對生態浮標數據進行多參數協同分析后發現，異常判定方法的關鍵是相關性計算時所要選取的時間窗口以及基于相關平穩性異常的判定方法。

圖2 研究區域浮標分布Fig.2 Distribution of buoys in the study area

表1 浮標數據統計Table 1 Statistical buoys data

3.1 時間窗口的確定

由于浙江沿岸海域生化參數日變化動態范圍較大，如以太短或者太長的時間段內的兩兩相關性來建立成段異常數據方法，則存在較大的隨機與不確定性，不利于對長時序浮標數據的穩定性研究與漸變異常數據的早期識別。因此，選擇合適的時間窗口，對于建立相關性分析處理異常數據模型至關重要。本文將浙江沿岸6 處儀器運行正常狀態的175 305 組浮標數據經過不可能出現的范圍和5S 方法剔除異常數據等預處理后，剩余156 305 組浮標數據參與多參數協同分析。其中，選出13 620 余組各參數質量較好的浮標數據對其進行兩兩相關性分析。部分正確的浮標數據序列如圖3 所示。

將圖3 的pH、DO 濃度和Chla濃度數據的兩兩相關系數（Rnd）計算的時間窗口逐天擴大，從1 d 擴大到16 d，結果如圖4a 所示。由圖可見，隨著時間窗口的擴大，相關性逐步提升，并且當擴大到8 d 時Rnd都處于穩定狀態，即當時間窗口大于8 d 后相關性并未明顯增強。同時以8 d 為時間窗口對圖3 中的多組浮標長時序數據進行8 d 兩兩相關系數（R8d）的計算（圖4b），可以看出正常原始浮標數據的R8d在一定時期內同樣非常穩定，因此時間窗口選定為8 d，同時將R8d作為檢測漸變異常數據的核心參數。

3.2 基于相關性的異常判定方法

圖3 部分正確的浮標數據序列Fig.3 Partially correct buoy data sequence

圖4 不同時間窗口的相關系數（a）和基于8 d 時間窗口的相關系數（b）Fig.4 Correlation coefficient for different time windows (a),and correlation coefficient for 8 d time window (b)

首先，利用多參數之間相關性程度來進行異常數據判定。如圖4b 所示，正常數據的R8d變化平穩，狀態穩定。為定量化正常R8d變化的范圍，利用6 處浮標多年中的正常數據集，統計了浙江海域R8d的數值分布情況，如圖5 和表2 所示。統計結果表明：（1）pH與DO 濃度之間正相關性最強，幾乎所有正常數據的R8d(pH-DO)都大于0；（2）DO 濃度與Chla濃度之間相關性次之，其R8d(DO-Chla) 大于-0.3，其中大于0 的數據近95%；（3）pH 與Chla濃度之間相關性變化較大，但有92% 數據的R8d(pH-Chla) 大于0，另6.8%的R8d(pH-Chla)數據在-0.3～0 之間，并且僅有1.2%的R8d(pH-Chla) 小于-0.3。因此，（1）正常數據pH 與DO 濃度、pH 與Chla濃度、DO 濃度與Chla濃度明顯存在較高的正相關關系，其判定原則較為簡單，即兩項以上的相關性R8d都大于0.5 可以作為數據正常有效標志；（2）因pH 與DO 濃度之間不存在負相關關系，明顯錯誤數據的判定原則為當R8d(pHDO)＜0 時為異常值；另外當R8d(pH-DO)＞0 時，如果Chla濃度與DO 濃度、pH 之間不存在較強的負相關關系，即R8d(pH-Chla)＜-0.3 時，或R8d(DO-Chla)＜-0.3，可識別為可疑數據。實際上，單一的R8d只能用于識別相對明確的正確及異常數據，而對于其中某項相關性R8d小于0.5 的浮標數據需要進一步采用其他相關性時序穩定性指標來進行識別。

圖5 R8d（a）和ΔR（b）的分布情況Fig.5 Distribution of R8d (a) and ΔR(b)

表2 R8d 的分布情況Table 2 Distribution of R8d

如前文所述，本文認為正常浮標多參數數據之間的兩兩相關性在一定時序上是穩定甚至是一致的，因此需要建立一個指標來表征R8d本身的穩定性。本文利用前后兩天R8d之差的絕對值（ΔR）作為判斷相關性時序分布穩定與否的指標。通過統計正常浮標數據的ΔR分布情況（圖5b，表3）可見，其中約有60%～70%數據的ΔR＜0.06，且ΔR＜0.1 數據都達到了77.6%以上。從數據分布上可以看出，正常浮標數據的多參數之間相關性變化是穩定或緩變的過程，符合前文的穩定性假設。由于計算求得ΔR的標準差為0.068，同時從表3 的統計結果也可看出，各相關性中ΔR＞0.34的數據僅占1.0%左右，因此選取5 倍標準差[19]即0.34為判斷穩定性閾值，即當0＜R8d(pH-DO)＜0.5，-0.3＜R8d(DO-Chla，pH-Chla)＜0.5 時，有一項ΔR＞0.34則判定為異常值。

表3 ΔR 的分布情況Table 3 Distribution of ΔR

利用R8d和ΔR兩項指標進行漸變異常數據的判斷與識別流程如圖6 所示。第一步利用單一指標R8d來判定簡單的正確數據和異常數據；第二步則是利用ΔR作為R8d穩定性指標來進一步判定異常數據。

4 結果與討論

4.1 識別算法實例與結果對比分析

圖6 數據判斷流程圖Fig.6 Flow chart of buoy data processing

圖7 2015 年4-6 月（a-c）和2014 年6-7 月（d-f）TZ01 浮標的原始數據（a，d）、R8d（b，e）和ΔR（c，f）Fig.7 TZ01 buoy raw data (a,d),R8d (b,e),and ΔR(c,f) in April to June,2015 (a-c) and June to July,2014 (d-f)

本文利用浙江溫州、臺州及舟山海域NJ01 浮標、TZ01 浮標以及ZS04 處浮標典型數據，對漸變異常數據的判定方法進行了適用性驗證。首先，選取同樣位于臺州外海TZ01 浮標的2015 年4-6 月（圖7a-c）和2014 年6-7 月（圖7d-f）的兩組正常數據。第一組2015 年的原始數據是十分具有代表性的正確數據，Chla濃度、DO 濃度和pH 之間存在非常高的正相關性，兩兩R8d大于0.5，并且相關性的變化平穩，ΔR都小于0.34。而第二組2014 年的正確數據相比于第一組數據變化更加復雜，從6 月下旬，Chla濃度與DO濃度和pH 存在極弱的相關性，R8d(DO-Chla)和R8d(pHChla)接近于0，但隨著7 月初藻華事件的出現，上述R8d逐漸升高，并在整個藻華期間處在一個平穩的高相關性時期。雖然在這一過程中R8d的總體變化很大，但根據ΔR的計算結果都小于0.34，可以說明這個變化過程是穩定的漸變過程。那么根據圖6 的識別方法仍然可以準確判定為正確數據，因此證明了本文方法的適用性。

本研究同樣利用了浙江臺州海域TZ01 浮標2013年5-6 月（圖8a-c），以及溫州海域NJ01 浮標2014 年6 月（圖8d-f）和2015 年3-4 月（圖8h-j）的3 組存在漸變異常的數據集對本文識別方法進行了適用性驗證。

第一組案例是2013 年5-6 月TZ01 浮標數據（圖8a），其中，5 月初有一次藻華事件，3 個參數變化同步浮標數據正常，而漸變異常數據實際上出現在5 月24 日前后，pH 上升發生偏離，后續在5 月30 日前后恢復正常。圖8b 和圖8c 分別給出了對應的R8d和ΔR數據，圖中紅色為異常值區間（5 月24-29 日），灰色部分為相關性計算受異常值影響區間。可以看出5 月15 日出現ΔR＞0.34 的情況，但是根據圖6 判斷流程，5 月15 日3 組R8d都升高到0.5 以上，因而仍然判定為正確數據。而在5 月24 日（表4），R8d(pH-DO)下降到-0.42，R8d(Chla-pH)下降到-0.39，并且ΔR(pH-DO)為0.89，大于0.34，多項指標符合本文漸變異常數據的判定標準，因此成功判定為異常數據，實現了數據異常早期識別。另外DO-Chla的R8d雖然有所下降，但仍然保持在0.5 以上。如果利用排除法，本方法可以進一步判定是3 個參數中pH 值出了問題。

圖8 2013 年5-6 月TZ01 浮標原始數據（a）、R8d（b）和ΔR（c），2014 年6 月（d-f）和2015 年3-4 月（g-i）NJ01 浮標原始數據（d，g）、R8d（e，h）和ΔR（f，i）Fig.8 TZ01 buoy raw data (a),R8d (b),and ΔR(c) in May to June,2013,and NJ01 buoy raw data (d，g),R8d (e，h),and ΔR(f，i) in June 2014 (d-f) and March to April 2015 (g-i)

第二組案例是2014 年6 月NJ01 浮標的數據，如圖8d 所示漸變異常數據出現在6 月9 日前后。根據圖8e 和圖8f 對應的R8d和ΔR結果（表4），在6 月9 日雖然R8d（DO-Chla）和R8d（pH-Chla）的下降并未超過-0.3，但是DO-Chla和pH-Chla的ΔR分別高達0.41、0.44，皆大于0.34，故可判斷出在6 月9 日數據出現了異常。根據后續Chla濃度的變化也可看出6 月9 日是異常數據出現較早時期，證明了本文方法早期識別的有效性。

表4 浮標出錯日期的R8d 和ΔR 情況Table 4 R8d and ΔR of buoy error date

第三組案例是2015 年3 月和4 月NJ01 浮標數據（圖8g-i），同樣可以看出漸變異常數據開始出現在4月7 日的前后。然而與前兩組數據不同的是，本組不同的R8d出現相反的變化趨勢（表4），其中R8d（DO-Chla）降到了-0.41，而R8d（pH-Chla）卻上升到了0.45，但兩者的ΔR都大于0.34，最終4 月7 日被判定為異常數據起始點。通過上述多組案例的驗證，本文的識別方法能夠適用于浙江沿海多參數浮標數據的漸變異常類型識別。

4.2 水體生化參數變化帶來滯后現象的影響

圖9 2015 年9 月ZS04 浮標原始數據（a）、R8d（b）和ΔR（c）Fig.9 ZS04 buoy raw data (a),R8d (b),and ΔR(c) in September,2015

海洋水體生化特性變化并不完全同步，一些參數相比于其他參數具有滯后現象。部分大型赤潮發生時，如圖9a 所示的舟山2015 年9 月25 日前后的一次赤潮事件，藻類暴發導致葉綠素峰值的出現往往早于DO 濃度或pH。這就導致在初期，葉綠素濃度與DO 濃度或pH 的R8d數值相對較小（圖9b），同時隨著赤潮的發展，R8d迅速升高，導致ΔR可能大于閾值0.34（圖9c）。此種情況下，如果R8d都高于0.5 呈現極高正相關性亦可判定為正確數據，但是如果小于0.5 則很容易被錯誤識別為異常數據。在這種劇烈而快速的海洋過程中，如何有效識別正常海洋規律現象與漸變異常數據是十分重要但又極具挑戰性的問題。實際上，本文采用的方法是時間同步相關性分析。如果利用時間延遲模式或許可以有效避免此類異常數據的錯誤識別。然而，目前對于該類滯后現象產生的海洋學機制并不十分明了。因此，在識別方法中如何具體引入時間延遲模式（如延遲區間等）還需進一步研究。

4.3 季節性（冬季）相關特征差異的影響

本文數據主要集中在上半年，而季節性變化（特別是冬季）對近海海域的海洋現象有著重要影響。由于浙江沿岸受河流沖淡水、季風和各類水團影響較大，冬季浙江沿海海域受浙閩沿岸流影響，水中懸浮泥沙含量高，限制了藻類生長。在這一時期，水體Chla濃度、pH 和DO 濃度的相關性比春、夏和秋3 個季節要弱很多。圖10 為2014 年冬季（2014 年11 月至2015 年2 月）TZ01 浮標的觀測結果，從圖10a 可看出，3 個參數的時序變化依然連續平穩，ΔR也基本在0.34 以內（圖10c），也說明了數據的平穩性，但是相關性系數R8d時高時低（圖10b），甚至出現極強的負相關。主要原因是在冬季水溫低，藻類豐富度較小，導致葉綠素濃度變化很小，對pH 和DO 濃度的影響作用有限。反而在這一時期，DO 濃度的變化受溫度影響較大，有個緩慢上升過程。因此，冬季浙江海域的3 個參數在機理上并不存在明確的相關性，而本文方法也僅基于同一年份前序時間數據的相關性進行漸變異常數據識別，所以在冬季可能會失效。

圖10 2014 年冬季TZ01 浮標原始數據（a）、R8d（b）和ΔR（c）Fig.10 TZ01 buoy raw data (a),R8d (b),and ΔR (c) in the winter of 2014

在上述不適用的情況下，我們需對時序相關性的概念進一步拓展，可利用同一海域季節性數據存在物候等現象，依靠歷史同一時期觀測數據集等，對浮標漸變異常數據進行有效地識別。如劉增宏等[34]采用歷史水文觀測資料集得到的溫-鹽度關系對Argo 剖面浮標鹽度資料進行校正，王輝贊等[35]也同樣通過尋找Argo 浮標不同剖面位置與其“最佳匹配”歷史剖面資料對比判別的途徑，對Argo 浮標鹽度偏移現象進行有效甄別。上述方法雖然用的是連續深度剖面數據，但是替換成連續時間序列數據同樣適用。如圖11所示，為臺州大陳浮標在2014-2015 年冬季與2015-2016 年冬季的pH、DO 濃度和Chla濃度數據對比結果，可以看出這兩年冬季同一時期的pH、DO 濃度和Chla濃度數據變化有較好的一致性趨勢。因此利用與多年歷史數據的相關性，可對pH、DO 濃度和Chla濃度數據進行異常識別。這或許是本文識別方法在冬季失效問題的一種有效解決方式，但需要大量歷史數據的積累，目前還無法實現，需要進一步研究。

5 結論

本研究通過對浙江沿岸6 處浮標多年多參數觀測數據進行分析，發現了與傳統跳變異常數據不同的漸變異常數據類型。該異常數據類型呈現出在時序變化過程連續平穩，但隨時間逐漸偏移，最后整體偏離正常數據的分布特征；并且在異常發生的初始階段，其變化特征與由海洋環境變化引起的變化趨勢很難在沒有先驗知識的條件下進行區分。因此本文提出了一種假設：浮標觀測到的正常多參數數據不僅單一參數在一定時序上的變化是平穩連續的，并且兩兩參數間的相關性在一定時序上是穩定甚至是一致的。根據上述假設，本文建立了基于pH、DO 濃度、Chla濃度數據兩兩相關性的漸變異常數據類型自動識別方法，確定了以8 d 時間窗口的兩兩相關系數（R8d）作為核心相關性表征指標，并將前后兩天R8d之差的絕對值（ΔR）作為判斷相關性時序分布穩定性指標，形成了利用R8d和ΔR兩項指標進行漸變異常數據判斷與識別的流程。

圖11 TZ01 浮標冬季原始數據Fig.11 TZ01 buoy raw data in winter

本文提出的方法重點突出了多元參數間相關性系數時間序列上的變化特征，各判別指數計算過程簡單、直觀，易于實際浮標監測工作人員的理解和掌握。通過浙江沿海浮標實際測量數據案例檢驗，證明了該方法可以用于漸變異常數據類型的實時監測，對浮標的傳感器漸變異常做到早期識別，特別是由生物污垢導致傳感器測量值持續增加而引起的假赤潮現象，有較好的識別效果，可解決由此帶來的赤潮預報虛警等問題。因此，在指導浮標日常檢查與維護、確保數據的準確性和完整性方面有實際意義。本文根據單參數自校方法無法識別漸變異常數據類型，提出了一種簡單，實用的有效解決方法。此方法為漸變異常值的自動識別及處理提供了新的思路。由于所處海域的不同，可能相關性穩定的時間窗口有所不同，需因地制宜，考慮季節性差異等因素的影響。因此在后續研究中應當針對在多參數變化不同步、冬季數據和非高斯分布數據等情況下，識別精度不高等局限性，可利用多年歷史數據對其進行物候特征分析，提高相關性識別方法精度。