鄭 旭，馮 晨，李騰躍，程 倫，何 波

(中國海洋大學信息科學與工程學院水下機器人實驗室，山東 青島 266100)

海洋資源的開發(fā)及全球自然環(huán)境保護的需要，極大程度推動了人們對于自主式水下機器人(Autonomous Underwater Vehicle,AUV)的研究[1-2]。AUV自身結(jié)構(gòu)復雜，系統(tǒng)的非線性與不確定性強，所處的海洋環(huán)境變化莫測[3]，這些因素會對AUV狀態(tài)監(jiān)測數(shù)據(jù)產(chǎn)生很大的干擾，因此，對AUV數(shù)據(jù)的波動模式進行準確描述有較大的難度，國內(nèi)外的相關(guān)研究也一直在不斷地進行[4-7]。

近年來，復雜系統(tǒng)分析方法—復雜網(wǎng)絡(luò)，被廣泛地應用于復雜系統(tǒng)研究中，基于復雜網(wǎng)絡(luò)的非線性時間序列特征的挖掘是其中一個重要的研究方向[8-11]。這是因為，復雜網(wǎng)絡(luò)在刻畫整體波動情況的同時，也可以對時間片段的波動模式進行精確的描述，有助于更清晰地了解時間序列的波動模式及規(guī)律。但實際監(jiān)測數(shù)據(jù)中往往包含許多噪聲，直接使用原始時間序列構(gòu)建網(wǎng)絡(luò)會對分析結(jié)果產(chǎn)生很大干擾[12-13]，因此，需要對數(shù)據(jù)進行預處理，才能使構(gòu)建出的復雜網(wǎng)絡(luò)中的節(jié)點更加準確、有效的描述時間序列中關(guān)鍵的波動模式。

近年來，隨著國內(nèi)外學者對時間序列預處理研究的深入，涌現(xiàn)了許多的優(yōu)秀算法[14-15]。分段累積近似算法(PAA)[16]成為一種常用的數(shù)據(jù)預處理方法，在此基礎(chǔ)上，Eammon又提出了符號化的相似性度量(SAX)方法[17]。該方法可以在一定程度上去除干擾，并以符號化的方式表明時間序列的大小特性，但缺乏對于序列模式之間轉(zhuǎn)化的描述，且需要數(shù)據(jù)服從高斯分布，不具有普適性。另外PAA算法在去除噪聲的同時也平滑掉了一些關(guān)鍵的峰值點，會使整體的波動模式發(fā)生變化，對分析結(jié)果產(chǎn)生不利影響。定性趨勢分析方法(Qualitative Trend Analysis，QTA)[18-19]被用于時間序列的描述中，該方法可以用符號化的方式準確表現(xiàn)出時間序列的波動模式，但無法排除冗雜波動的干擾，會產(chǎn)生許多不必要的波動模式。因此，需要探尋到一種新的時間序列預處理算法，既可以保留關(guān)鍵的峰值信息，又可以去除冗雜波動的干擾，并能適用大多數(shù)的數(shù)據(jù)。

針對AUV這類復雜的非線性系統(tǒng)的航行狀態(tài)描述問題，本文構(gòu)建AUV航行監(jiān)測時間序列復雜網(wǎng)絡(luò)，分析AUV時間序列波動模式。在構(gòu)建復雜網(wǎng)絡(luò)過程中，針對實際數(shù)據(jù)中噪聲干擾嚴重的問題，本文使用改進的預處理方法去除噪聲干擾，能夠優(yōu)化復雜網(wǎng)絡(luò)構(gòu)建，使復雜網(wǎng)絡(luò)更加準確地描述AUV系統(tǒng)狀態(tài)。

1 預處理算法與復雜網(wǎng)絡(luò)構(gòu)建方法

1.1 改進的預處理算法(Pro-DPCA)

為更好的實現(xiàn)數(shù)據(jù)的預處理，既保留關(guān)鍵的轉(zhuǎn)折點，又能去除多余的波動，本文提出了一種改進的數(shù)據(jù)預處理算法，算法分為映射、聚類、復原三步。首先，將長時間序列看作若干個長度為3的短序列，根據(jù)映射規(guī)則將短序列投影到二維平面內(nèi)，用二維平面的坐標表示段時間序列的波動模式。然后，利用密度峰值聚類算法對平面內(nèi)的點進行聚類，聚類以后，波動模式相似的點會被分到同一類。最后，根據(jù)復原函數(shù)重新將二維平面的點轉(zhuǎn)化為時間序列。

1.1.1 映射 設(shè)需要處理長度為n的時間序列{f1,f2,…,fn}，將其劃分為多個長度為3的段時間序列fi-1,fi,fi+1(若長度n不是3的倍數(shù)，可以適當增加或減少以一個數(shù)據(jù)點，一般來說n遠大于3)。映射后的參數(shù)作為二維平面的橫縱坐標X,Y，每個短時間序列映射到二維平面的結(jié)果如圖1所示。

圖1 映射過程示意圖Fig.1 Mapping process diagram

選取映射參數(shù)：將映射后的橫坐標定義為中間點到首尾兩點連線的距離，縱坐標為首尾兩點的高度差。具體公式如下，式中xi，yi分別表示第i個數(shù)據(jù)點fi的橫縱坐標值。

Xi=

i=2,…,n-1;

(1)

Yi=yi+1-yi-1,i=2,…,n-1。

(2)

這樣選取映射參數(shù)的優(yōu)勢如下：

(1)兩個參數(shù)之間不存在直接的相關(guān)性。

(2)兩個參數(shù)的取值范圍相同。

(3)兩個參數(shù)的變化可以完備的表現(xiàn)出3個數(shù)據(jù)點的所有變化形式。

(4)對于增大、保持、減小等不同模式的變化敏感，且映射后相反模式之間的距離最大。

1.1.2 密度峰值聚類(DPCA) 聚類采用Alex Rodriguez和Alessandro Laio提出的基于密度峰值的聚類方法[20]。其主要思想是尋找被低密度區(qū)域分離的高密度區(qū)域。對于一個數(shù)據(jù)集，聚類中心被一些低局部密度的數(shù)據(jù)點包圍，而且這些低局部密度的點距離其他有高局部密度的點的距離都比較大[21]。在這樣的模型中，DPCA主要有兩個需要計算的量：第一，局部密度ρi；第二，與高密度點之間的距離δi。局部密度ρi的定義為：

ρi=∑jχ(dij-dc)。

(3)

距離δi的定義為：

(4)

圖2 樣本聚類結(jié)果Fig.2 Sample of clustering results

1.1.3 復原 在第一步中，原始時間序列被轉(zhuǎn)化成了諸多長度為3的片段，并將它們映射到二維平面中，現(xiàn)在需要將其復原，也就是根據(jù)聚類的結(jié)果把它們從二維平面轉(zhuǎn)化為時間序列，復原公式如下：

(5)

式中Z表示復原后的序列，若第一個點被分為class1，那么Z1,2,3=f1。

圖3展示了不同算法對含噪聲時間序列的預處理結(jié)果。通過對比可以看到，本文提出的算法不僅保留了關(guān)鍵的峰值點，而且成功的消除了無關(guān)波動的影響，保留了數(shù)據(jù)關(guān)鍵的波動模式。這是因為雖然映射后的結(jié)果與基元表示法[22]都能完備的表現(xiàn)所有波動模式，但經(jīng)過聚類以后，一部分波動模式會根據(jù)其波動幅度大小以及形狀的不同，更加靠近數(shù)據(jù)主要的波動模式(如上升、下降、保持等)，從而達到消除冗雜波動的目的。而部分符號化的復原方式則很好的避免了在PAA算法中，處理后的數(shù)據(jù)無法保留原數(shù)據(jù)的關(guān)鍵峰值，仍有很多冗雜波動。另外，PAA算法預處理后的數(shù)據(jù)并不符合高斯分布，因此無法直接使用SAX算法進行后續(xù)處理。

圖3 不同預處理算法的結(jié)果Fig.3 The results of different pretreatment algorithm

1.2 復雜網(wǎng)絡(luò)構(gòu)建

(6)

然后定義滑動窗口n，滑窗內(nèi)n個符號構(gòu)成的短模式組作為網(wǎng)絡(luò)的節(jié)點，模式組之間的轉(zhuǎn)換構(gòu)成網(wǎng)絡(luò)的連邊，連邊的權(quán)重ω即為兩個模式組之間轉(zhuǎn)化的次數(shù)，構(gòu)建有向加權(quán)復雜網(wǎng)絡(luò)，表1顯示了復雜網(wǎng)絡(luò)的構(gòu)建流程，滑窗的選取有所重疊，一定程度上包含了前一片段的信息，使得片段之間具有記憶性和傳導性同時不缺乏多樣性[24]。

表1 情況下的復雜網(wǎng)絡(luò)構(gòu)建流程Table 1 The process of complex networks building when n=3

2 Lorenz系統(tǒng)的實驗與分析

2.1 Lorenz映射系統(tǒng)

Lorenz映射系統(tǒng)是最早發(fā)現(xiàn)的混沌運動的耗散系統(tǒng)，許多學者對它也進行了去噪方面的研究[25]。Lorenz系統(tǒng)方程表示見公式7，當ρ>24.74時，系統(tǒng)呈現(xiàn)非周期的混沌態(tài)。使用該系統(tǒng)產(chǎn)生的時間序列，驗證本文方法的有效性和先進性。

(7)

2.2 數(shù)據(jù)與網(wǎng)絡(luò)構(gòu)建

選取參數(shù)為δ=18，ρ=28，β=8/3的Lorenz系統(tǒng)中連續(xù)1000個x值進行分析。首先使用本文提出的預處理算法對數(shù)據(jù)進行預處理，結(jié)果見圖4，然后使用1.2節(jié)中的方法進行復雜網(wǎng)絡(luò)的構(gòu)建，結(jié)果見圖5。

對網(wǎng)絡(luò)結(jié)構(gòu)進行分析可知，原始數(shù)據(jù)的網(wǎng)絡(luò)(見圖5(a))主要集中在上升下降兩種模態(tài)，且節(jié)點轉(zhuǎn)換模式主要集中在這兩種模態(tài)的自轉(zhuǎn)換中。而使用本文方法構(gòu)建的網(wǎng)絡(luò)(見圖5(b))在結(jié)構(gòu)上與原始數(shù)據(jù)類似，但出現(xiàn)了一種新的模態(tài)“MMM”。從時間序列的角度來看，它表示的是原始數(shù)據(jù)中上升或下降幅度較小的部分,即系統(tǒng)的平衡點附近；從曲線運動的角度來看，這一模態(tài)反應的是曲線圍繞吸引子進行圓周運動時方向改變的部分以及即將轉(zhuǎn)入另一吸引子的部分(見圖6)[26]。因此本文算法構(gòu)建的網(wǎng)絡(luò)在準確描述數(shù)據(jù)的波動模式的基礎(chǔ)上，對于Lorenz映射系統(tǒng)的物理意義也有著更深入的挖掘。

圖4 Lorenz系統(tǒng)輸出x值及其預處理后的數(shù)據(jù)Fig.4 x value output by Lorenz system and data preprocessed by Pro-DCPA

(其中節(jié)點大小代表節(jié)點強度高低，連線深淺代表權(quán)值大小。The higher the strength and weight ,the larger the node and the deeper the line in the net.)圖5 復雜網(wǎng)絡(luò)拓撲圖Fig.5 The topological of complex network

圖6 lorenz系統(tǒng)XOY平面圖Fig.6 XOY plane of the Lorentz chaos system

2.3 噪聲對網(wǎng)絡(luò)的影響

對原始數(shù)據(jù)添加信噪比為8 dB的高斯白噪聲(見圖7)進行實驗，分析噪聲對不同網(wǎng)絡(luò)結(jié)構(gòu)產(chǎn)生的影響，從而驗證本文算法構(gòu)建的網(wǎng)絡(luò)對噪聲的魯棒性。

對于原始數(shù)據(jù)構(gòu)建的網(wǎng)絡(luò)，加入噪聲后，網(wǎng)絡(luò)的節(jié)點數(shù)明顯增多，而且多個連邊的權(quán)值較高(見圖8(a))，對網(wǎng)絡(luò)結(jié)構(gòu)造成了較大影響。而用本文方法構(gòu)建的網(wǎng)絡(luò)加入噪聲后(見圖8(b))，雖然模態(tài)有所增加，但其連邊權(quán)值明顯很低，網(wǎng)絡(luò)仍保持原有的結(jié)構(gòu)。

①并行服務(wù)快捷高效。應用服務(wù)子系統(tǒng)一方面充分利用目前服務(wù)器多核處理器的優(yōu)勢，在IIS服務(wù)器和空間數(shù)據(jù)服務(wù)器中都采用了并行服務(wù)方式，能夠同時啟動多個進程并行處理多個數(shù)據(jù)服務(wù)請求，提高了Web應用服務(wù)的效率；另一方面利用ArcGIS的SilverLight客戶端能夠?qū)崿F(xiàn)多線程并行計算的特點，對于用戶提交的復雜應用請求可以分解為多線程計算處理任務(wù)，分步提交顯示計算結(jié)果，滿足了用戶Web查詢快捷高效的應用要求。

使用節(jié)點自轉(zhuǎn)換概率為指標，定量說明網(wǎng)絡(luò)的穩(wěn)定性。節(jié)點的轉(zhuǎn)換概率，即節(jié)點i轉(zhuǎn)換為節(jié)點j的概率，定義為prij：

圖7 Lorenz系統(tǒng)輸出x加噪數(shù)據(jù)及其預處理后的數(shù)據(jù)Fig.7 x value output by Lorenz system with noise and data preprocessed by Pro-DCPA

圖8 復雜網(wǎng)絡(luò)Fig.8 Complex networks

(8)

式中：ωij為節(jié)點i轉(zhuǎn)換到j(luò)的權(quán)重；nsi表示節(jié)點i的強度。當i=j時，該式表達的則是節(jié)點的自轉(zhuǎn)換概率。表2為使用不同預處理算法構(gòu)建網(wǎng)絡(luò)中關(guān)鍵節(jié)點的自轉(zhuǎn)換概率?？梢钥闯?，對于未加噪數(shù)據(jù)來說，本文算法構(gòu)建的網(wǎng)絡(luò)與使用原始數(shù)據(jù)構(gòu)建的網(wǎng)絡(luò)差異不明顯，但加入噪聲以后，隨著滑窗長度n的增加，原始數(shù)據(jù)構(gòu)建的網(wǎng)絡(luò)中，關(guān)鍵節(jié)點自轉(zhuǎn)換概率明顯下降，但使用本文算法構(gòu)建的網(wǎng)絡(luò)仍舊保持較高的自轉(zhuǎn)換概率，特別是在n=10的時候，相較原始數(shù)據(jù)提升了15%以上。

表2 網(wǎng)絡(luò)關(guān)鍵節(jié)點自轉(zhuǎn)換概率Table 2 The transition probabilities of the major node in the different net

注：表中3H表示3個連續(xù)H構(gòu)成的節(jié)點。Here “3H” means the node “HHH”.

綜上，通過Lorenz系統(tǒng)的實驗與分析表明，使用本文算法構(gòu)建的網(wǎng)絡(luò)，不僅可以更加準確地表現(xiàn)系統(tǒng)的波動模式，還可以挖掘出更深的物理意義，而且對尺度和噪聲都具有更強的魯棒性。

3 AUV數(shù)據(jù)波動模式分析

3.1 數(shù)據(jù)獲取及網(wǎng)絡(luò)構(gòu)建

研究選取AUV下潛狀態(tài)下的實際速度數(shù)據(jù)進行實驗，采樣頻率為1 Hz，如圖9所示，實驗選取深度增加部分的速度數(shù)據(jù)。然后使用本文提出的方法構(gòu)建有向加權(quán)網(wǎng)絡(luò)(見圖10)。網(wǎng)絡(luò)中節(jié)點表示了數(shù)據(jù)的波動情況，從理論上來說，通過不同模式的組合，共有33=27種可能的模式組出現(xiàn)，然而實際只出現(xiàn)了其中20種，使用原始數(shù)據(jù)構(gòu)建的網(wǎng)絡(luò)會有更多的模式出現(xiàn)，但更多的模式往往是因為一些無關(guān)的波動所引起。

圖9 AUV實際深度與速度Fig.9 The actual speed and depth of AUV

圖10 AUV下潛階段實際速度數(shù)據(jù)構(gòu)建的網(wǎng)絡(luò)(滑窗長度n=3)Fig.10 The topological of complex network based on actual speed during AUV submergence stage(sliding window n=3)

3.2 網(wǎng)絡(luò)特性分析

基于構(gòu)建好的網(wǎng)絡(luò)，利用節(jié)點分布的統(tǒng)計學特性對網(wǎng)絡(luò)進行分析，挖掘數(shù)據(jù)中的關(guān)鍵信息，掌握AUV的運動狀態(tài)。首先觀察每個模式所占的比例，了解數(shù)據(jù)整體結(jié)構(gòu)。圖11顯示了網(wǎng)絡(luò)中各模式的占比，其中M模式所占比例最高，這表明多數(shù)數(shù)據(jù)處于“保持”狀態(tài)，AUV運行狀態(tài)相對平穩(wěn)。

圖11 各模式所占比例Fig.11 The proportion of fluctuating models

接下來以節(jié)點為單位，使用節(jié)點強度nsi對網(wǎng)絡(luò)進行分析，節(jié)點強度不僅考慮了節(jié)點的的相鄰個數(shù)，同時也考慮了與相鄰節(jié)點的權(quán)重

nsi=∑j∈Nωij。

(9)

式中：nsi為第i個節(jié)點的強度；ωij表示節(jié)點i和j之間的權(quán)重。

表3顯示了使用Pro-DPCA算法構(gòu)建的復雜網(wǎng)絡(luò)中，強度前5的節(jié)點占比情況。其中“MMM”占了最大的比例66.3%，遠高于其他節(jié)點。從圖12可以看出，節(jié)點強度服從冪律分布，并呈現(xiàn)長尾特性，這表明主要波動模式集中在少數(shù)模態(tài)中，其他模態(tài)只占了很小的比例，對于該時間序列則主要集中在“保持”這一模態(tài)。這反映出AUV在下潛過程中基本保持勻速行駛，與實際情況相符合。

表3 節(jié)點強度及所占比例(n=3,Pro-DPCA算法)Talbe 3 Proportion of node strength(n=3,Pro-DPCA)

接下來選擇不同的滑窗長度n，并使用不同的預處理算法進行網(wǎng)絡(luò)構(gòu)建，進行對比實驗，觀察Pro-DPCA算法與其他預處理算法以及原始數(shù)據(jù)所構(gòu)成網(wǎng)絡(luò)之間的差異，分析網(wǎng)絡(luò)在不同尺度下是否依然可以準確表現(xiàn)數(shù)據(jù)的波動模式。

表4列舉了3種算法在不同滑窗長度n下的關(guān)鍵節(jié)點占比情況，其中PAA算法選取的片段長度為3，因此節(jié)點數(shù)為原始的1/3。通過對比可以看出，使用本文算法構(gòu)建的復雜網(wǎng)絡(luò)可以準確的挖掘出數(shù)據(jù)中的大量“保持”狀態(tài)，原始數(shù)據(jù)由于波動較多，很難表示出數(shù)據(jù)應有的波動模式，而PAA算法雖然對波動進行了平滑，但效果仍不理想。其中在n=3的時候PAA算法構(gòu)成的網(wǎng)絡(luò)達到了27個節(jié)點，表現(xiàn)了所有可能的波動模式，這也就意味著滑窗長度n=3這個尺度對于片段長度為3的PAA算法處于一個飽和的狀態(tài)，進一步增加片段長度會使該算法在這里有更好的表現(xiàn)，但是這樣也平滑掉了更多的關(guān)鍵峰值，而且會將序列維度進一步降低，對于其他的分析起到不利的作用[27]。另一方面，隨著滑窗長度n的增加，PAA算法中的關(guān)鍵節(jié)點占比下降迅速，當n=10的時候不足1%，而使用本文算法構(gòu)建的網(wǎng)絡(luò)仍能有60%的比例，可以更好的適應滑窗尺度的變化。

(其中R為節(jié)點按強度大小排序的秩。HereRis the serial numbers that are sorted by their node strength value.)
圖12 節(jié)點強度分布
Fig.12 The distribution of node strength

表4 關(guān)鍵節(jié)點占比Table 4 The proportion of node strength in different net

注：表其中節(jié)點強度指的是網(wǎng)絡(luò)中關(guān)鍵節(jié)點(連續(xù)n個”M”)的強度。The node strength is refer to the strength of the major node(ntimes of “M” continuity).

綜上，通過實驗和分析表明，本文提出的方法可以在含有大量冗雜波動的AUV數(shù)據(jù)中挖掘出主要的波動模式，準確表現(xiàn)出AUV的航行狀態(tài)，并且對于尺度的變化有著更強的魯棒性。

實際上，不僅限于AUV，在現(xiàn)實生活中的很多領(lǐng)域，例如風場信號、通信信號、工業(yè)中的滾動軸承信號等等，這些信號的時間序列也都是可以使用復雜網(wǎng)絡(luò)方法分析系統(tǒng)狀態(tài)的。而這些實際數(shù)據(jù)中往往也伴都隨著噪聲的干擾，本文的算法同樣適用于這些領(lǐng)域。因此，本文提出的算法應用不僅局限于AUV數(shù)據(jù)的去噪，還有著其他更為廣泛的實際應用。

4 結(jié)語

本文使用復雜網(wǎng)絡(luò)的方法描述復雜非線性系統(tǒng)的狀態(tài)，并提出了一種映射與密度峰值聚類結(jié)合的改進預處理算法，以此優(yōu)化復雜網(wǎng)絡(luò)的構(gòu)建。在實驗中使用Lorenz系統(tǒng)和AUV實際航行數(shù)據(jù)對本文方法進行了驗證。結(jié)果表明，與傳統(tǒng)方法相比，本文的方法可以更為準確地描述時間序列的波動情況，深入挖掘數(shù)據(jù)中的關(guān)鍵信息，且對噪聲干擾和尺度變化具有更強的魯棒性。本文的研究為去除時間序列數(shù)據(jù)中的噪聲干擾這一問題提供了新的可行方案，在數(shù)據(jù)挖掘和復雜系統(tǒng)建模等領(lǐng)域有十分廣泛的應用前景。