于芳星， 姬 波， CHENG Quanrun， 盧紅星， 柳宏川

(1.鄭州大學 信息工程學院，河南 鄭州 450001； 2.伍倫貢大學 電子計算機與通信工程學院，澳大利亞 新南威爾士州 伍倫貢市 2522)

0 引言

由光反饋效應逐漸形成的自混合(SMI)干涉理論被廣泛應用于位移計算、震動檢測、形貌重構、粒子運動測量等領域[1]。Fischer等[2]首次提出光反饋干涉儀中的雙腔結構(雙腔OFI系統(tǒng))由激光二極管(LD)和兩個外部高反射金鏡組成。該模型在研究具有許多自由度的非線性系統(tǒng)中的動力學時，具有十分理想的效果。Yu等[3]和Fan等[4]先后研究了多重光反饋以及多模激光器的自混合干涉理論，提出了含預反饋的激光自混合干涉型位移測量結構，可以在弱反饋條件下得到類鋸齒波以進行判向。近年來，研究人員已在許多應用中使用了雙腔OFI結構：Jiang等[5]和Geng等[6]提出了一種基于雙腔OFI的多普勒速度測量方法； Zhang等[7]和Chen等[8]使用雙腔OFI實現(xiàn)了二維振動測量和多個目標的運動檢測；Mezzapesa等[9]提出了一種使用雙腔OFI來實現(xiàn)納米級位移傳感的方法。

Lang等[10]首次提出一種描述單模半導體激光器的速率方程：Lang-Kobasyashi(L-K)方程，它可以完整地描述具有外部光反饋(EOF)的單模激光的動態(tài)行為，成為研究光反饋自混合干涉效應動態(tài)問題的最經典方法之一。L-K方程本質是一個常微分方程組，具有變系數(shù)、非線性的特征，因此無法獲得解析解(符號解)，需要對其進行數(shù)值求解。目前，已有的研究對雙腔OFI系統(tǒng)中的L-K方程主要通過歐拉法[11]、梯形公式法[12]、改進歐拉法[13]、四階龍格-庫塔法[14]進行求解。但這些方法存在著求解精度不夠高、計算速度不夠快的問題，導致生成的光電信號不能夠十分準確地描述激光的動態(tài)行為，降低了系統(tǒng)的測量精度。因此，為了進一步提高光電信號雙腔OFI系統(tǒng)的動態(tài)行為求解精度，本文提出一種求解L-K方程的六階龍格-庫塔算法，通過選取更多的區(qū)間點計算積分曲線的斜率平均值，使其更接近于真實值以提高精度，并將其應用于雙腔OFI系統(tǒng)的移動物體運動檢測仿真軟件中進行仿真實驗。

1 雙腔OFI系統(tǒng)的L-K方程

如圖1所示，雙腔OFI系統(tǒng)由激光二極管(LD)、光電二極管(PD)、聚焦透鏡(Lens)、目標1(Target-1)、目標2(Target-2)組成。其中Taerget-1為被測量目標，其表面的反射率非常低，LD與Target-1之間的腔稱為被測量腔。Target-2表面具有較高的反射率可以用于提供預反饋，LD與Target-2之間的腔稱為控制腔。

圖1 雙腔OFI系統(tǒng)Figure 1 Dual-cavity OFI system

L-K方程包括3個聯(lián)立的延遲微分方程(DDE)，通過修改L-K微分方程，可以擴展得到雙腔L-K方程：

E(t-τ2)·cos(ω0τ2+φ(t)-φ(t-τ2))；

(1)

(2)

(3)

雙腔L-K方程中參數(shù)的物理意義和數(shù)值見表1。

表1 L-K方程中參數(shù)的物理意義Table 1 Physical meaning of the parameters in the L-K equation

2 雙腔OFI系統(tǒng)L-K方程的龍格-庫塔算法

2.1 龍格-庫塔算法

在自然科學和社會科學的許多領域，經常會遇到一階微分方程的初值問題：

(4)

在一系列離散的點x1,x2,x3，…，xn上求出未知函數(shù)y(x)之值y(x1),y(x2),y(x3),…,y(xn)的近似值y1,y2,y3，…，yn，即為所求的微分方程初值問題的數(shù)值解[15]。龍格-庫塔算法是目前最為常見的一種微分方程初值問題的數(shù)值解法，相比于歐拉法、梯形公式法、改進歐拉法具有較高的精度和穩(wěn)定性，而且容易編程，所以被廣泛應用于各個領域。

龍格-庫塔算法的基本思想是在(xn,xn+1)之間取一些積分曲線的若干點的切線斜率，再進行一次(或多次)算數(shù)(或加權)平均后產生新的斜率，隨后按這個斜率從(xn,yn)以直線代曲線向前推進一步。根據(jù)差商代替導數(shù)的思想，由微分中值定理可知：

(5)

注意到一階微分方程y′(x)=f(x,y) 可得

y(xn+1)=y(xn)+hf(xn+θh,y(xn+θh))。

(6)

2.2 四階龍格-庫塔算法

當在區(qū)間內取4個點計算平均斜率時，即為四階龍格-庫塔算法。四階龍格-庫塔公式并不唯一，其中最經典的為[16]

(7)

2.3 六階龍格-庫塔算法

六階龍格-庫塔算法需要計算7個k值，根據(jù)上述推導過程，可以求得一種六階龍格-庫塔公式[17]：

(8)

2.4 L-K方程的龍格-庫塔算法

根據(jù)六階龍格-庫塔公式，由式(1)可以得到L-K方程的龍格-庫塔公式：

(9)

根據(jù)式(9)，可以得到L-K方程六階龍格-庫塔算法的偽代碼如下所示。

算法L-K方程的六階龍格-庫塔算法。

輸入：E、φ、N關于t的微分方程和初始值；

輸出：E、φ、N值。

① 設置E、φ、N的初始值E0、φ0、N0；

② 設置步長h；

③ for (k=1∶SamplingNumber) do

④ 計算微分方程中目標的相位變化值；

⑤ 根據(jù)六階龍格-庫塔公式計算k1、k2、k3、k4、k5、k6、k7；

⑥ then

⑦ 依次計算E、φ、N；

⑧ 更新E0、φ0、N0；

⑨ end for

用n表示采樣數(shù)的值，那么六階龍格-庫塔算法的時間復雜度為O(n)。

3 仿真軟件設計和實現(xiàn)

雙腔OFI系統(tǒng)仿真軟件通過MATLAB進行開發(fā)，用于測量移動物體在運動中的動態(tài)行為，觀察仿真結果是否符合預期并預判測量精度，最終指導硬件平臺上的實際物理測試實驗過程。

雙腔OFI系統(tǒng)仿真軟件進行模塊化后可以分為7個主要模塊，其模塊分解圖見圖2。

圖2 模塊分解圖Figure 2 Module diagram

雙腔OFI系統(tǒng)仿真軟件中的接口包括常量、變量、計算反饋強度、計算相位變化量、定義微分方程、龍格-庫塔算法、繪制圖像7個接口，接口定義示例見表2。

表2 常量接口Table 2 Constant interface

雙腔OFI系統(tǒng)仿真軟件的總體仿真流程為：通過六階龍格-庫塔算法求解L-K方程，然后根據(jù)所得的E(腔內電場強度)值繪制OFI信號關于時間t的圖像。

4 仿真實驗與結果分析

本實驗的主要工作是根據(jù)電場強度E的值體現(xiàn)L-K方程的求解精度，不同反饋強度下的實驗為對照組。因此，實驗的目標是觀察反饋強度對干涉效應的影響及不同反饋強度下的測量精度。k1、k2的值由實驗輸入?yún)?shù)反饋系數(shù)c1、c2求解得出。

4.1 單腔實驗結果

在單腔實驗中，設置c1=0，分別設置c2=0.3、c2=2、c2=3進行3組實驗。作為對照，實驗分別采用歐拉法、四階龍格-庫塔算法(RK4)、六階龍格-庫塔算法(RK6)進行求解。實驗對比結果見圖3、圖4，圖中橫坐標t代表時間，縱坐標OFI signal代表系統(tǒng)生成的光反饋信號強度。

圖3 單腔實驗中歐拉法與RK6的對比圖Figure 3 Comparison of Euler method and RK6 in single cavity experiment

圖4 單腔實驗中RK4與RK6的對比圖Figure 4 Comparison of RK4 and RK6 in single cavity experiment

表3和表4展示了單腔實驗中不同算法求解的部分E值和精度改進率。通過計算可得，RK6相比于歐拉法的平均精度改進率約為22.95%，RK6相比于RK4的平均精度改進率約為6.01%。

對比實驗結果分析表明，在歐拉法與RK6的對比實驗中，當c2較小時，其信號幅度也比較小，歐拉法和RK6的插補效果幾乎相同；當c2分別增大到2和3時，歐拉法會出現(xiàn)較明顯的偏差，RK6的效果明顯優(yōu)于歐拉法。在RK4與RK6的對比實驗中，當c2較小時，RK4和RK6的插補效果幾乎相同；當c2增大到3時，RK4會出現(xiàn)較明顯的偏差，RK6的效果明顯優(yōu)于RK4。

表3 單腔實驗中歐拉法與RK6的部分E值Table 3 Some E values of Euler method and RK6 in single cavity experiment

表4 單腔實驗中RK4與RK6的部分E值Table 4 Some E values of RK4 and RK6 in single cavity experiment

4.2 雙腔實驗結果

在雙腔實驗中，設置3組c1、c2的值分別為：c1=0.1，c2=0.3；c1=1.5，c2=0.3；c1=4，c2=0.3。實驗對比結果見圖5、圖6。

圖5 雙腔實驗中歐拉法與RK6的對比圖Figure 5 Comparison of Euler method and RK6 in dual-cavity experiment

圖6 雙腔實驗中RK4與RK6的對比圖Figure 6 Comparison of RK4 and RK6 in dual-cavity experiment

表5和表6展示了在雙腔實驗中不同算法求解的部分E值和算法精度改進率。通過計算可得，RK6相比于歐拉法的平均精度改進率約為22.94%，RK6相比于RK4的平均精度改進率約為6.01%。

表5 雙腔實驗中歐拉法與RK6的部分E值Table 5 Some E values of Euler method and RK6 in dual-cavity experiment

表6 雙腔實驗中RK4與RK6的部分E值Table 6 Some E values of RK4 and RK6 in dual-cavity experiment

對比實驗結果分析表明：在雙腔實驗中，其總體趨勢與單腔實驗相似，即RK6的效果明顯優(yōu)于歐拉法和RK4；但是當c2=0.3，同時c1逐漸增大時，信號波形浮動范圍進一步增大，會導致歐拉法、RK4和RK6的插補效果都變得惡劣。

5 結論

本文通過分析微分方程數(shù)值解法的原理，選取更多區(qū)間點計算積分曲線的斜率平均值，提出一種求解雙腔OFI系統(tǒng)中L-K方程的六階龍格-庫塔算法，并在移動物體運動檢測仿真軟件中進行了對比實驗。結果表明：六階龍格-庫塔算法與歐拉法相比，求解精度平均提高了約22%，與四階龍格-庫塔算法相比，求解精度平均提高了約6%。因此，通過六階龍格-庫塔算法求解L-K方程可以提高其求解精度，從而提高光電信號雙腔OFI系統(tǒng)對于移動物體的高靈敏度感測精度。