姜敏敏羅文茂趙 力

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

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

常微分方程的數(shù)值求解一般分為顯式和隱式遞推方法[11]。 對(duì)于剛性常微分方程一般采用隱式方法數(shù)值求解。 對(duì)于非剛性常微分方程，一般采用顯式的定步長(zhǎng)4 階龍格庫(kù)塔法進(jìn)行數(shù)值求解。 采用定步長(zhǎng)是因?yàn)樽儾介L(zhǎng)龍格庫(kù)塔法的計(jì)算量大。 而采用4 階精度則是因?yàn)榈碗A的求解精度不夠，而高階方法超出了應(yīng)用的實(shí)際需求，會(huì)浪費(fèi)算力。 所以采用非剛性常微分方程形式的耦合混沌振子進(jìn)行傳感信號(hào)檢測(cè)的文獻(xiàn)中普遍采用定步長(zhǎng)4 階龍格庫(kù)塔法。

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

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

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

1 半隱式遞推格式

為了更具廣泛性，考慮兩個(gè)帶有時(shí)間變量的2階耦合振子，其表達(dá)式為:

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2 仿真研究

2.1 耦合Lorenz 振子

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

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

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

圖1 龍格庫(kù)塔法的相圖

圖2 半隱式法的相圖

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

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

2.2 耦合Duffing 振子

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

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

將以上表達(dá)式改寫(xiě)為1 階微分方程組形式:

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

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

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

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

3 實(shí)驗(yàn)研究

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

圖5 局放信號(hào)采集場(chǎng)景

實(shí)驗(yàn)采集到的局放信號(hào)波形如圖6 所示。

圖6 實(shí)驗(yàn)采集的局放信號(hào)

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

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

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

4 結(jié)論

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