用于傳感信號檢測的混沌振子數(shù)值并行求解方法*

姜敏敏羅文茂趙 力

(1.南京信息職業(yè)技術(shù)學(xué)院，江蘇 南京 210023；2.東南大學(xué)信息科學(xué)與工程學(xué)院，江蘇 南京 210096)

近年來，將混沌理論引入信號處理領(lǐng)域有了較大發(fā)展，在保密通信[1]、傳感信號檢測[2]、信號去噪[3]、圖像處理[4]等方面涌現(xiàn)出了眾多應(yīng)用成果。在應(yīng)用混沌振子時，一般采用單個振子、振子陣列和耦合振子的形式。 其中，耦合振子在傳感信號檢測領(lǐng)域中有眾多的研究成果，如文獻[5－10]所示。 在這些文獻中，普遍采用非剛性常微分方程形式的耦合混沌振子，如Lorenz 振子、Duffing 振子、Duffing-Van der Pol 振子等。 將耦合振子應(yīng)用于傳感信號檢測的過程中，一般要進行微分方程的數(shù)值求解。隨著耦合振子的增多，微分方程數(shù)量將變得龐大，所需求解運算量急劇增加。 傳感信號檢測的一個重要問題是求解的實時性，所以求解速度是制約耦合混沌振子應(yīng)用的一個因素。

常微分方程的數(shù)值求解一般分為顯式和隱式遞推方法[11]。 對于剛性常微分方程一般采用隱式方法數(shù)值求解。 對于非剛性常微分方程，一般采用顯式的定步長4 階龍格庫塔法進行數(shù)值求解。 采用定步長是因為變步長龍格庫塔法的計算量大。 而采用4 階精度則是因為低階的求解精度不夠，而高階方法超出了應(yīng)用的實際需求，會浪費算力。 所以采用非剛性常微分方程形式的耦合混沌振子進行傳感信號檢測的文獻中普遍采用定步長4 階龍格庫塔法。

但是，定步長4 階龍格庫塔法的遞推格式中，一個方程的求解需順序計算4 次，運算速度較慢。 所以，提高數(shù)值計算速度是一個需要研究的問題。

半隱式遞推方法在非線性問題中得到了很多研究[12－15]，但是半隱式遞推格式必須根據(jù)不同非線性問題的特性來優(yōu)化設(shè)計，需在穩(wěn)定域、精度、計算量等方面進行考量，而沒有標(biāo)準(zhǔn)算法。

針對傳感信號檢測時耦合混沌振子的數(shù)值求解問題，本文給出了一種半隱式的并行計算方法，其計算速度是同階次定步長龍格庫塔法的一倍。 仿真和實驗證明，該方法對于常見的非剛性耦合Lorenz、Duffing 振子的求解結(jié)果與四階龍格庫塔法類似，且在實際傳感信號檢測任務(wù)中采用二階方法就能得到有效結(jié)果。

1 半隱式遞推格式

為了更具廣泛性，考慮兩個帶有時間變量的2階耦合振子，其表達式為:

式中:x、y是其中一個振子的狀態(tài)變量，u、v是另一個振子的狀態(tài)變量，f1、f2、g1、g2是函數(shù)，t是時間變量。

在2 階精度時，其半隱式遞推格式分為兩次遞推，分別為并行計算的式(2)和式(3):

式中:下標(biāo)n表示第n步的數(shù)據(jù)，下標(biāo)n＋1 表示第n＋1 步的數(shù)據(jù)，H為遞推步長。

最終的遞推結(jié)果為式(2)和式(3)的相應(yīng)結(jié)果的平均值。 如變量x的第n＋1 步的遞推結(jié)果xn＋1為式(2)和式(3)得到的xn＋1相加后除2，其余變量以此類推。 為了后續(xù)推導(dǎo)過程描述方便，將以上2 階遞推結(jié)果標(biāo)記為T1(xn＋1，yn＋1，un＋1，vn＋1)。

該半隱式方法具有2 階精度。 可以看出，式(2)中的第一、三個公式是顯式遞推格式，第二、四個公式是隱式遞推格式。

對于式(1)所示的微分方程組，其有4個方程，如果采用定步長2 階龍格庫塔法求解，每個時間步解算4個方程，共需順序執(zhí)行2個時間步才能得到結(jié)果。 所以每次遞推運算的時長是2個時間步。 式(1)如果采用以上半隱式遞推算法，式(2)和式(3)可以并行計算，每次遞推運算的時長是1個時間步。所以2 階半隱式并行算法的運算速度是2 階龍格庫塔法的一倍。

對于精度要求高的情況下，需提高遞推階數(shù)。以4 階半隱式遞推格式為例，其遞推分兩步分別執(zhí)行。

第一步的遞推方法如上述式(2)、式(3)描述的2 階方法，其遞推結(jié)果為T1(xn＋1，yn＋1，un＋1，vn＋1)。

第二步的遞推分為兩個小步驟，第1個小步驟的遞推格式如式(4)和式(5)所示:

式中:下標(biāo)n＋0.5 表示第n＋0.5 步的數(shù)據(jù)，上標(biāo)′表示第2 步遞推的結(jié)果。

第1個小步驟的最終遞推結(jié)果為式(4)和式(5)的相應(yīng)結(jié)果的平均值。 如變量x的第n＋1 步的遞推結(jié)果xn＋1為式(4)和式(5)得到的x′n＋0.5相加后除2，其余變量以此類推。 將以上第1個小步驟的遞推結(jié)果標(biāo)記為T21(xn＋0.5，yn＋0.5，un＋0.5，vn＋0.5)。

第二步的第2個小步驟的遞推格式如式(6)和式(7)所示:

式中:下標(biāo)n＋1 表示第n＋1 步的數(shù)據(jù)。

第2個小步驟的最終遞推結(jié)果為式(6)和式(7)的相應(yīng)結(jié)果的平均值。 如變量x的第n＋1 步的遞推結(jié)果xn＋1為式(6)和式(7)得到的x′n＋1相加后除2，其余變量以此類推。 將以上第2個小步驟的遞推結(jié)果標(biāo)記為T22(xn＋1，yn＋1，un＋1，vn＋1)。

在以上的遞推過程中，T1和T22的計算過程是完全獨立的，可以并行計算。T1的計算是1個時間步，T22的計算需2個時間步，所以其計算時長是2個時間步。 作為對比，定步長4 階龍格庫塔法的計算需要順序執(zhí)行4個時間步。 所以，在4 階精度下，本文方法的計算速度是定步長4 階龍格庫塔法的一倍。

以上是2 階和4 階的半隱式遞推方法，對于更高階的遞推方法也可以類似構(gòu)造。 另外，如果微分方程組的方程個數(shù)大于4，其遞推格式也可以類似構(gòu)造，一般要求隱式和顯式格式數(shù)量相等。

2 仿真研究

2.1 耦合Lorenz 振子

Lorenz 振子作為一種常見的混沌振子，其應(yīng)用非常廣泛。 文獻[16]提出一種非剛性耦合Lorenz振子系統(tǒng)，其表達式如下所示:

式中:x1、y1、z1是其中一個振子的狀態(tài)變量，x2、y2、z2是另一個振子的狀態(tài)變量，α＝10，β＝24.76，γ＝8/3，ε＝0.05。 在仿真中設(shè)定遞推步長H＝0.01，遞推10 000個點。

通過定步長4 階龍格庫塔法和4 階半隱式方法計算，得到的相圖分別如圖1 和圖2 所示。 通過對比可以看出，兩種方法得到的相圖是一致的，說明本文方法對于非剛性耦合Lorenz 振子的計算是精確的。

圖1 龍格庫塔法的相圖

圖2 半隱式法的相圖

為了進一步對比計算精度，將兩種方法對于變量x1的計算結(jié)果在不同的遞推點處進行逐點對比，其結(jié)果如表1 所示。 通過對比可以看出兩種方法在小數(shù)第8 位上有細(xì)小的差異，說明本文的半隱式方法在應(yīng)對一般問題時具有足夠的精度。

表1 計算結(jié)果對比

2.2 耦合Duffing 振子

Duffing 振子目前在傳感信號檢測領(lǐng)域應(yīng)用非常廣泛。 為了檢測脈沖信號，文獻[6]提出一種非剛性耦合Duffing 振子系統(tǒng)，其表達式如下所示:

式中:x1是其中一個振子的狀態(tài)變量，x2是另一個振子的狀態(tài)變量，s(t)為待檢測信號，n(t)為噪聲，F(xiàn)＝0.2，ω＝5×106，k＝10，ξ＝0.7，計算步長1 ns。

將以上表達式改寫為1 階微分方程組形式:

為了驗證該耦合振子對脈沖信號的檢測能力，仿真了一些疊加高斯白噪聲的方波信號，其波形如圖3(a)所示。 將該信號加入式(10)所示的耦合Duffing振子中，通過定步長4 階龍格庫塔法求解得到的結(jié)果(?x1－?x2)如圖3(b)所示。 同樣利用4 階半隱式方法求解，得到的結(jié)果如圖4(a)所示。 通過對比可以看出，兩種方法得到的結(jié)果是一致的，說明本文方法對于非剛性耦合Duffing 振子的計算是精確的。

圖3 龍格庫塔法求解結(jié)果

圖4 半隱式法求解結(jié)果

更進一步，本文探討了2 階精度的半隱式法對非剛性耦合混沌振子的應(yīng)用可能性。 應(yīng)用2 階半隱式法對耦合Duffing 振子的求解結(jié)果如圖4(b)所示。 從該圖可以看出，雖然2 階半隱式方法的求解精度不高，相比4 階精度方法有一定的數(shù)值差異，但是落實到具體的傳感信號檢測應(yīng)用中，可以發(fā)現(xiàn)2階半隱式方法同樣可以得到與4 階方法相似的信號檢測結(jié)果。 通過該算例可以得出2 階半隱式方法可以應(yīng)用于傳感信號檢測、可以進一步提升運算的速度的結(jié)論。

3 實驗研究

為了驗證本文方法在實際傳感信號檢測中的正確性，在高壓開關(guān)柜上搭建了局部放電信號檢測的實驗場景。 通過帶寬為300 MHz～3 GHz 的微帶天線采集高頻電磁信號，示波器采用LeCroy HDO9204，其帶寬為2 GHz。 實驗場景如圖5 所示。

圖5 局放信號采集場景

實驗采集到的局放信號波形如圖6 所示。

圖6 實驗采集的局放信號

將采集的局放信號輸入式(9)的Duffing 振子系統(tǒng)，通過4 階龍格庫塔法求解，得到去噪后的輸出信號如圖7 所示。

從圖7 可以看出，去噪后局放信號被凸顯出來，說明利用Duffing 振子對局放信號能有效去噪。 同時，對比4 階龍格庫塔法和本文提出的4 階半隱式方法對局放信號的去噪結(jié)果，可以看出兩種數(shù)值求解方法得到的結(jié)果是相似的。 這說明本文半隱式方法對實際傳感信號的求解也是成功的。

圖7 局放信號去噪求解結(jié)果

4 結(jié)論

本文針對傳感信號檢測所用非剛性耦合混沌振子給出了一種半隱式并行數(shù)值算法，該算法的計算速度是相同階數(shù)龍格庫塔法的一倍。 通過耦合Lorenz 振子和耦合Duffing 振子的算例，證明了本文方法和相同階數(shù)龍格庫塔法的求解精度是相似的。并且通過算例發(fā)現(xiàn)，在利用耦合混沌振子進行傳感信號檢測的應(yīng)用中，2 階半隱式方法也能得到很好的求解結(jié)果。 因此，將耦合混沌振子應(yīng)用于傳感信號檢測時可以嘗試應(yīng)用2 階半隱式方法求解，這樣可以實現(xiàn)更小的計算代價和更高的求解速度。 此外，局放信號的檢測實驗結(jié)果也表明了本文方法可以用于實際的傳感信號檢測環(huán)境。