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

姜敏敏羅文茂趙 力

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

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

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

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

半隱式遞推方法在非線性問題中得到了很多研究[12－15]，但是半隱式遞推格式必須根據不同非線性問題的特性來優化設計，需在穩定域、精度、計算量等方面進行考量，而沒有標準算法。

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

1 半隱式遞推格式

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

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

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

式中:下標n表示第n步的數據，下標n＋1 表示第n＋1 步的數據，H為遞推步長。

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

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

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

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

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

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

式中:下標n＋0.5 表示第n＋0.5 步的數據，上標′表示第2 步遞推的結果。

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

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

式中:下標n＋1 表示第n＋1 步的數據。

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

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

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

2 仿真研究

2.1 耦合Lorenz 振子

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

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

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

圖1 龍格庫塔法的相圖

圖2 半隱式法的相圖

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

表1 計算結果對比

2.2 耦合Duffing 振子

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

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

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

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

圖3 龍格庫塔法求解結果

圖4 半隱式法求解結果

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

3 實驗研究

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

圖5 局放信號采集場景

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

圖6 實驗采集的局放信號

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

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

圖7 局放信號去噪求解結果

4 結論

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