基于向量自回歸模型的高頻地波雷達(dá) 海流異常值的識別算法研究

方 今，劉京城

（1.上海市海洋監(jiān)測預(yù)報(bào)中心，上海 200062；2. 湖北中南鵬力海洋探測信息系統(tǒng)工程有限公司，湖北 宜昌 443005)

海洋作為社會(huì)經(jīng)濟(jì)發(fā)展的重要依托和載體，是國家未來發(fā)展的新方向。監(jiān)測海洋表面流場、風(fēng)場和浪場等重要的海洋動(dòng)力學(xué)要素，有助于保障海上和沿岸人類活動(dòng)和安全以及海防事業(yè)，有利于發(fā)展海洋經(jīng)濟(jì)、預(yù)防海洋自然災(zāi)害和維護(hù)國家利益。高頻地波雷達(dá)[1]可以全天候地對雷達(dá)波束內(nèi)的中遠(yuǎn)距離（如200 km）海面進(jìn)行實(shí)時(shí)監(jiān)測，在海態(tài)參數(shù)監(jiān)測中具有不可替代的作用。隨著高頻地波雷達(dá)在探測海洋表面動(dòng)力學(xué)要素方面的應(yīng)用和發(fā)展，海洋表面流場探測技術(shù)在海態(tài)檢測領(lǐng)域已經(jīng)趨于成熟，各方面的對比驗(yàn)證[2-3]表明其具有優(yōu)秀的反演性能。目前，國外成型的高頻地波雷達(dá)產(chǎn)品主要有美國SeaSoned系統(tǒng)、德國WERA系統(tǒng)、英國OSCR系統(tǒng)、俄羅斯“向日葵”系統(tǒng)和加拿大SWR-503系統(tǒng)等，國內(nèi)成型的海態(tài)高頻地波雷達(dá)產(chǎn)品主要分為陣列式高頻地波雷達(dá)OSMAR[4]和便攜式高頻地波雷達(dá)OSMAR-S[5]。其中陣列式高頻地波雷達(dá)分辨率更高、探測范圍更遠(yuǎn)，而便攜式高頻地波雷達(dá)采用單極子/交叉環(huán)天線，占地面積更小。本文中所選用的數(shù)據(jù)為OSMAR-S100系統(tǒng)反演得到的海流結(jié)果。

由于雷達(dá)站所處電磁環(huán)境一般比較復(fù)雜，電離層的干擾、廣播電臺的干擾、瞬態(tài)干擾以及其他干擾，都會(huì)影響雷達(dá)回波的質(zhì)量。雖然已有相應(yīng)的干擾抑制算法[6-7]，但這些算法不能一直奏效，因此，海流反演結(jié)果中經(jīng)常會(huì)出現(xiàn)時(shí)間或空間上的異常值，異常值的存在會(huì)影響海流流速和流向的可靠性。一般的海流異常值識別方法[8-9]是在單一維度中進(jìn)行的，能夠識別部分異常值，由于忽略了另一個(gè)維度，存在誤識別或漏識別的情況。

向量自回歸（VAR）模型[10]是一種多變量預(yù)測模型，從縱向比較來看，它是單變量時(shí)間序列自回歸（AR）模型在多元時(shí)間序列上的衍生；從橫向比較來看，它是用線性關(guān)系刻畫一個(gè)平穩(wěn)的系統(tǒng)。VAR模型的特性使其在平穩(wěn)的條件下，參數(shù)估計(jì)與最小二乘估計(jì)一致，并且統(tǒng)計(jì)上具有優(yōu)良性質(zhì)，便于區(qū)間估計(jì)、誤差分析和模型診斷等，可用于刻畫多變量之間相互影響的動(dòng)態(tài)線性相關(guān)關(guān)系，該模型應(yīng)用在風(fēng)速預(yù)測中[11-12]，取得了較好的結(jié)果。

本文基于海流流速在時(shí)間維度與空間維度是一個(gè)緩慢變化的過程，同一點(diǎn)的海流流速在時(shí)間維度與空間維度具有相關(guān)性，提出了基于VAR模型的海流異常值的識別算法，利用OSMAR-S100便攜式高頻地波雷達(dá)數(shù)據(jù)，結(jié)合時(shí)間域與空間域，通過VAR模型識別高頻地波雷達(dá)海流異常值，從而提高數(shù)據(jù)準(zhǔn)確度。

1 基于向量自回歸模型的海流異常值識別算法

1.1 海態(tài)流場數(shù)據(jù)來源

本文海流異常值識別數(shù)據(jù)來源于基于單極子交叉環(huán)緊湊型天線陣OSMAR-S100便攜式高頻地波雷達(dá)，該雷達(dá)部署于上海雞骨礁，坐標(biāo)位于（122.388°E，31.1730°N），工作頻率13 MHz，根據(jù)雷達(dá)指標(biāo)，海流最大探測距離80 km，波束覆蓋扇角范圍120°，海流探測的徑向角度和距離分辨率分別為3°和2.5 km，海流測量速度分辨率3 cm，該雷達(dá)系統(tǒng)每6.5 min輸出一場海流數(shù)據(jù)。

下面以雞骨礁站2020年7月5日數(shù)據(jù)為例，對雷達(dá)實(shí)際應(yīng)用時(shí)的判斷流程進(jìn)行介紹。該日一共產(chǎn)生了221場數(shù)據(jù)。經(jīng)過統(tǒng)計(jì)，數(shù)據(jù)點(diǎn)中存在大量缺失情況，圖1展示了數(shù)據(jù)完整性統(tǒng)計(jì)的結(jié)果，橫縱坐標(biāo)代表經(jīng)緯度，數(shù)值代表一天中有效數(shù)據(jù)的數(shù)量。圖1結(jié)果表明，越靠近雷達(dá)站，數(shù)據(jù)完整性越高，越遠(yuǎn)離雷達(dá)站，數(shù)據(jù)完整性越低，在20個(gè)距離元之后的數(shù)據(jù)基本可以判定不存在或者數(shù)據(jù)可靠性較差。在進(jìn)行數(shù)據(jù)分析時(shí)，選擇的位置均位于20個(gè)距離元以內(nèi)，盡量避免數(shù)據(jù)缺失所帶來的不利影響。

圖1 有效數(shù)據(jù)分布圖

1.2 自回歸模型

OSMAR-S100雷達(dá)系統(tǒng)一天數(shù)據(jù)可視為一個(gè)時(shí)間序列，海流流速具有時(shí)間相關(guān)性，某一時(shí)間海流流速與前后兩個(gè)時(shí)間點(diǎn)流速具有時(shí)間相關(guān)性。圖2展示了位置處于數(shù)據(jù)點(diǎn)（122.367°E，30.866°N）的海流觀測值對于不同滯后階數(shù)的自相關(guān)系數(shù)，海流流速自相關(guān)系數(shù)隨著滯后階數(shù)的增加而緩慢減小，這種自相關(guān)系數(shù)隨滯后階數(shù)增加而緩慢減小的特征被稱為“ACF圖拖尾”，是AR模型的顯著特征，即AR模型適用于海流觀測值的處理。

圖2 海流流速不同滯后階數(shù)的自相關(guān)系數(shù)圖

AR模型是一種分析具有時(shí)間序列特性數(shù)據(jù)的重要模型，可以表示為:

式中：c為常數(shù)；φi為回歸系數(shù)；p為AR模型的階數(shù)；εt為隨機(jī)誤差，通常認(rèn)為是白噪聲序列，階數(shù)p可從偏自相關(guān)系數(shù)隨滯后階數(shù)的變化情況估計(jì)得到，即利用AR模型的偏相關(guān)圖（PACF圖）截尾特征進(jìn)行估計(jì)，其中偏相關(guān)性（PACF）是兩個(gè)隨機(jī)變量在排除了其余部分或全部隨機(jī)變量影響情形下的凈相關(guān)性或純相關(guān)性，是兩個(gè)隨機(jī)變量在處于同一體系的其余部分或全部隨機(jī)變量取給定值的情形下的條件相關(guān)性，偏相關(guān)性由偏相關(guān)系數(shù)決定。圖3展示了上述地理位置的海流流速AR模型PACF圖。平穩(wěn)序列的偏相關(guān)圖表現(xiàn)為截尾，當(dāng)階數(shù)為1的時(shí)候，偏相關(guān)系數(shù)值還是很大，為0.814；二階長的時(shí)候突然就變成了0.125，后面的值都很小，認(rèn)為是趨于0，這種狀況就是截尾，可作為海流流速AR模型數(shù)據(jù)序列平穩(wěn)的判斷。

圖3 海流流速AR模型PACF圖

圖中紅線為顯著性閾值，超出該值表示該階滯后階數(shù)顯著不為0，反之則表示該階系數(shù)不顯著，可以理解為近似為0。海流觀測值的偏自相關(guān)系數(shù)從二階開始下降到紅線以下，因此海流觀測值符合一階自回歸模型。

1.3 VAR模型

AR模型同樣僅考慮了時(shí)間維度，為了同時(shí)考慮時(shí)間維度和空間維度，我們選取了空間上5個(gè)較近位置的海流觀測序列進(jìn)行了分析，所選擇5個(gè)位置的空間相對位置示意圖如圖4所示。

圖4 選取海流的空間相對位置示意圖

統(tǒng)計(jì)了數(shù)據(jù)有效范圍內(nèi)一周的所有數(shù)據(jù)點(diǎn)（黑色圓）的海流序列與其一階滯后序列以及周圍較近的4個(gè)位置（白色圓）的海流序列間的相關(guān)系數(shù)，統(tǒng)計(jì)結(jié)果的直方圖如圖5所示。其中數(shù)據(jù)點(diǎn)與其一階滯后序列的相關(guān)系數(shù)代表了時(shí)間維度上的相關(guān)性，與周圍的海流觀測序列的相關(guān)系數(shù)代表了空間維度上的相關(guān)性。圖中，80%以上數(shù)據(jù)點(diǎn)的前述5個(gè)相關(guān)系數(shù)都大于0.9，表明海流的觀測值在時(shí)間維度和空間維度中都具有較強(qiáng)的相關(guān)性。因此，我們對上述5個(gè)序列建立了VAR模型，階數(shù)仍沿用AR模型中確定的階數(shù)，即建立的模型為一階VAR模型，其表達(dá)式如下。

圖5 所有數(shù)據(jù)點(diǎn)時(shí)空相關(guān)系數(shù)分布圖

式中：

T為轉(zhuǎn)置運(yùn)算符；Xi為對應(yīng)的海流序列；φi為對應(yīng)的回歸系數(shù)，其中，φf為前數(shù)據(jù)點(diǎn)的海流序列回歸系數(shù)；φb為后數(shù)據(jù)點(diǎn)的海流序列回歸系數(shù)；φs為數(shù)據(jù)點(diǎn)自身海流序列回歸系數(shù)；φl為左數(shù)據(jù)點(diǎn)的海流序列回歸系數(shù)；φr為右數(shù)據(jù)點(diǎn)的海流序列回歸系數(shù)；Xt-1為數(shù)據(jù)點(diǎn)前一時(shí)刻海流數(shù)據(jù)，數(shù)據(jù)點(diǎn)X通過該時(shí)刻其前、后、左、右與前一時(shí)刻自身海流數(shù)據(jù)對現(xiàn)海流數(shù)據(jù)驗(yàn)證。

對式（2）稍做修改，得到式（5）。

式中：

對于式（5），可以利用最小二乘法求解其系數(shù)向量θ。

2 海流異常值識別

根據(jù)海流的時(shí)間維度與空間維度相關(guān)性提出的一階VAR模型是一種線性預(yù)測模型，根據(jù)式（2），通過某位置前一時(shí)刻的海流以及臨近4個(gè)位置的海流，可以預(yù)測該位置當(dāng)前時(shí)刻的海流。預(yù)測海流與真實(shí)海流之間的差值被定義為殘差，如式（9）所示。

式中：r為殘差序列；X為真實(shí)海流序列；為預(yù)測海流序列，理論上，認(rèn)為r是均值為0、方差不隨時(shí)間改變的白噪聲序列。

殘差序列的標(biāo)準(zhǔn)差反映了殘差整體的性質(zhì)，殘差標(biāo)準(zhǔn)差越大說明模型擬合效果越差，殘差中還可能殘留系統(tǒng)部分的信息。對于單個(gè)殘差點(diǎn)而言，如果它的絕對值超過了3倍殘差標(biāo)準(zhǔn)差，則說明該點(diǎn)在99%的顯著性水平上不等于0，不滿足白噪聲序列的假設(shè)，則可以標(biāo)記為海流異常值，對于標(biāo)記出來的海流異常值，用VAR模型的估計(jì)值代替原始的異常值，從而實(shí)現(xiàn)對海流異常值的識別及修正。

3 結(jié)果與分析

本文通過對位于上海雞骨礁的OSMAR-S100便攜式高頻地波雷達(dá)站海流反演數(shù)據(jù)進(jìn)行分析與處理，驗(yàn)證了基于VAR模型的高頻地波雷達(dá)海流異常值識別算法的有效性，并統(tǒng)計(jì)了海流異常值在空間和時(shí)間上的分布情況。

在對海流進(jìn)行處理時(shí)，以天為間隔進(jìn)行處理，共處理了一周的數(shù)據(jù)。針對檢測點(diǎn)每一天的海流時(shí)間序列，利用最小二乘法估計(jì)出模型參數(shù)，建立一階VAR模型，圖6對比了一個(gè)位置海流時(shí)間序列VAR模型估計(jì)的海流與雷達(dá)輸出的海流，二者結(jié)果非常接近，表明利用VAR模型預(yù)測某一時(shí)刻的海流切實(shí)可行。兩條曲線之間的差值即為殘差，根據(jù)殘差的結(jié)果判斷并標(biāo)記出所有的海流異常值。圖7展示了部分被標(biāo)記為異常海流值的殘差序列以及相鄰的5個(gè)海流序列，海流曲線中紅色實(shí)線代表檢測點(diǎn)的海流序列，虛線部分代表周圍4個(gè)點(diǎn)的海流序列，殘差曲線中為判決閾值，超出紅線范圍的殘差即被標(biāo)記為異常值，對應(yīng)位置的海流被標(biāo)記為異常海流值。圖7結(jié)果表明，當(dāng)檢測點(diǎn)的變化趨勢與周圍點(diǎn)的變化趨勢明顯不同，或者檢測點(diǎn)本身出現(xiàn)較劇烈的變化時(shí)，檢測點(diǎn)會(huì)被標(biāo)記成為異常點(diǎn)，這與人們對于海流的認(rèn)識（海流在空間和時(shí)間上是緩慢變化的）基本一致，即在時(shí)間維度和空間維度上變化劇烈的海流為異常海流值。

圖6 VAR模型預(yù)測海流與雷達(dá)輸出海流對比

圖7 VAR模型識別的部分異常值及其殘差

本文對一周內(nèi)的海流數(shù)據(jù)進(jìn)行了異常識別，結(jié)果表明對于孤立的異常值，該算法具有較好的識別效果，但是對于連續(xù)的異常值，可能由于后續(xù)變化趨勢一致而無法被識別出來，見圖7（a）和圖7（b），這種比較復(fù)雜的情況需要進(jìn)一步改進(jìn)預(yù)測模型予以判斷。同時(shí)，本文統(tǒng)計(jì)了海流異常值在時(shí)間和空間上的分布情況，統(tǒng)計(jì)結(jié)果如圖8所示。受雷達(dá)波束和距離衰減的影響，在空間上海流異常值的數(shù)量隨距離的增加而增加，且偏離陣列法線方向越遠(yuǎn)，異常值數(shù)量越多。而其在時(shí)間上的分布，近似于均勻分布，僅在少數(shù)時(shí)刻異常值數(shù)量偏大，表明該場次數(shù)回波受到了較強(qiáng)的外部干擾。

圖8 海流異常值在空間維度和時(shí)間維度上的分布

4 結(jié) 論

高頻地波雷達(dá)在海洋表面環(huán)境觀測方面具有超視距、大范圍、全天候的特性。針對高頻地波雷達(dá)海流結(jié)果中的異常值，本文提出了基于VAR模型的高頻地波雷達(dá)海流異常值的識別算法。由于本文算法通過空間與時(shí)間維度對海流流速進(jìn)行趨勢分析，在區(qū)域流速較大和較小時(shí)，這種異常點(diǎn)識別方法依舊能正常工作。通過對區(qū)域海流流速進(jìn)行基于VAR模型的時(shí)間維度分析，可提高高頻地波雷達(dá)海流數(shù)據(jù)的準(zhǔn)確性。整體上看，基于VAR模型的高頻地波雷達(dá)海流異常值的識別算法能夠?yàn)楹Ａ鳟a(chǎn)品和流場時(shí)空分布分析提供準(zhǔn)確的數(shù)據(jù)。