基于倒向隨機微分方程的光纖陀螺穩(wěn)定性評估方法

秦偉亮，魏法杰

（北京航空航天大學 經(jīng)濟管理學院，北京 100191）

光纖陀螺是一種基于Sagnac原理的純固態(tài)慣性儀表，具有環(huán)境適應性好、可靠性高、壽命長、綜合性能優(yōu)等特點[1]，較適合應用于空間領域。國外公開報道光纖陀螺產(chǎn)品壽命已超過15萬小時，國內(nèi)光纖陀螺在衛(wèi)星中連續(xù)工作時間也已超過數(shù)萬小時。目前光纖陀螺已在商業(yè)衛(wèi)星、軍用衛(wèi)星等領域取得了較廣泛應用。但由于空間環(huán)境的特殊性，光纖陀螺性能指標會隨著通電時間的積累而逐步退化[2]，進而影響到衛(wèi)星的使用壽命。光纖陀螺在空間環(huán)境條件下的參數(shù)穩(wěn)定性和衛(wèi)星服役的匹配性，已日益成為光纖陀螺競爭力的一個重要表現(xiàn)。如果陀螺壽命遠超過衛(wèi)星壽命，則陀螺成本相對會高，競爭力降低；如果陀螺壽命達不到衛(wèi)星壽命要求，則直接影響衛(wèi)星使用。因此，開展對光纖陀螺參數(shù)穩(wěn)定性的評估，以及對宇航光纖陀螺出廠指標參數(shù)進行控制的課題研究，對光纖陀螺的質量控制和推廣應用有較強的意義。

當前，對宇航光纖陀螺的評估主要集中在可靠性和壽命評估等方面，所采用的主要方法是地面等效壽命試驗或加速壽命試驗，并利用數(shù)理統(tǒng)計方法對產(chǎn)品剩余壽命進行預測，進而得出產(chǎn)品的指標首達設定閾值時的穩(wěn)定性時間分布[3]。

目前研究的特點是：根據(jù)光纖陀螺在空間環(huán)境條件下的退化機理，建立產(chǎn)品的退化模型，并通過試驗數(shù)據(jù)和仿真求解模型參數(shù)，利用參數(shù)和產(chǎn)品性能退化的首達時間獲得產(chǎn)品壽命估計。這是一種隨時間變化的正向評估方法，可用于光纖陀螺在出廠后預期壽命的后評估。但正向評估方法的周期較長、經(jīng)濟代價較大，不適合于對每只使用的光纖陀螺都進行評估。基于此，本文在評價方法上采取新的途徑：依據(jù)光纖陀螺在空間環(huán)境條件下的退化機理，通過引入倒向微分方程的理論，建立了相應模型，并通過模型求解，實現(xiàn)對產(chǎn)品的性能變化特征進行評估，實現(xiàn)對出廠參數(shù)進行控制，使其參數(shù)穩(wěn)定性同衛(wèi)星壽命相匹配。

1 倒向隨機微分方程基本理論

倒向隨機微分方程概念最早是由法國隨機控制和隨機分析專家Bismut[4]在1978年研究隨機最優(yōu)控制過程中提出來的。Bismut提出的概念是與正向系統(tǒng)方程相對偶的線性倒向隨機微分方程。1990年我國學者彭實戈院士（Peng）和法國學者Pardoux共同提出了非線性倒向隨機微分方程，并證明了其解的存在唯一性[5]。此后，倒向隨機微分方程的研究取得了長足的進步，在隨機控制、金融數(shù)學、遞歸效用和風險敏感效用等領域獲得了廣泛應用。

為引出倒向隨機微分方程，先對比分析在區(qū)間[0,T]上正向和倒向常微分方程形式：

式中，b(·)、g(·)是給定的函數(shù)；x0、yT是給定的參數(shù)。

式(1)的定解條件在初始時刻T=0時刻給出，稱它為正向常微分方程。式(2)的定解條件在終了時刻t=T時刻給出，稱它為倒向常微分方程。

引入隨機因素，即建立正、倒向隨機微分方程。兩者對比，其結構形式為：

式(3)是正向隨機微分方程，式(4)是倒向隨機微分方程基本形式。兩式聯(lián)立，構成了一組正倒向隨機微分方程組，其中，b(·)、g(·)、σ(·)、z(·)是給定的函數(shù)。

在數(shù)學上，分析求解隨機微分方程一般需要測度論知識。(Ω,Ft,Pt)是一個含有信息流Ft的完備概率空間，F(xiàn)t是由(Bs,s≤t)生成的σ代數(shù)，即Ft=σ(Bs,s≤t)，Bt是Ft上的標準的 Brown運動，Pt是概率空間上的概率測度。在式(4)中，ω是一個確定的可測得值，X(t)是狀態(tài)變量，Y(t)和z(t)是倒向隨機過程中同時間t相關的、需要同時解決的兩個隨機過程。g（X(t) ,Y(t) ,z(t) ,t）是倒向隨機微分方程的生成元，是關于t、X(t)、Y(t)、z(t)的函數(shù)。

在工程實際上，一般過程都滿足理論要求條件：給定時間區(qū)間[0,T]上，隨機過程連續(xù)平方可積且有界，也即Y(t)、z(t)都是連續(xù)平方可積且有界的。在符合該條件時，正、倒向隨機微分方程都滿足解的存在唯一性[5]。

倒向隨機微分方程解得存在唯一性和生成元表示定理奠定了倒向隨機微分方程的理論基礎。Pardoux-Peng[6]建立的非線性Feynman-Kac公式，將求解生成元中的Y(t)和z(t)隨機過程與一個擬線性偏微分方程組的解建立了關系，從而獲得倒向隨機微分方程的唯一解表示。在工程應用上，Y(t)是可觀測的隨機過程，z(t)可用數(shù)學期望形式表示為：

式中，Δ表示為采樣時間間隔，為條件數(shù)學期望，O(Δ)為Δ的高價小量。為便于表達，式中及下文將Y(t)簡化表示為Yt，X(t)表示為Xt，z(t)表示為zt。

從式(5)中可以看出，z(t)的含義是Y(t)波動特性的表示。

與正向隨機微分方程相比，倒向隨機微分方程在結構上有較大不同，它們的主要區(qū)別和聯(lián)系表現(xiàn)在：

1）正向隨機微分方程有兩個決定函數(shù)b(·)和σ(·)，其解是個隨機過程X(t)，而倒向隨機微分過程是由生成元g(X(t),Y(t),z(t),t)決定，它的解是一組隨機過程Y(t) ,z(t)；

2）正向隨機微分過程關注如何認識一個客觀存在的隨機變化過程，而倒向微分方程關注在隨機干擾條件下如何使一個產(chǎn)品或系統(tǒng)達到預期的要求，也就是通過將來時刻給定的一個目標值來獲得當前（產(chǎn)品出廠）時刻的值。正向隨機微分方程的值代表今天確定的值變?yōu)槊魈煲话悴淮_定的狀態(tài)以研究其統(tǒng)計規(guī)律，而倒向隨機微分方程的值則表示將明天的目標變成今天的解以制定今天的決策[6]。

3）為了獲得倒向隨機微分方程的解，可建立一組正倒向隨機微分方程組。根據(jù)隨機最優(yōu)控制理論和最大值原理，正向隨機微分方程描述了系統(tǒng)的狀態(tài)方程，引入一類倒向隨機微分方程可作為狀態(tài)方程的伴隨方程，共同形成一組正倒向隨機微分方程組[7]。通過對方程組的求解獲得倒向隨機微分方程的值。

2 光纖陀螺倒向微分方程模型

2.1 基本模型

目前，在衛(wèi)星領域應用的光纖陀螺一般采用干涉式閉環(huán)信號處理方案，它在衛(wèi)星環(huán)境條件下的性能退化特性在統(tǒng)計上可用隨機維納過程描述。光纖陀螺光源作為有源器件，其輸出功率穩(wěn)定性和噪聲特性直接影響到陀螺的參數(shù)穩(wěn)定性[2]，在陀螺的性能退化和壽命評估中具有重要作用。對光源輸出功率和噪聲的動力學特性分析表明，在穩(wěn)態(tài)條件下的均勻介質中，光量子的運動軌跡是個馬爾科夫過程，其輸出功率和噪聲特性符合Langevin模型。

Langevin模型是一個正向隨機微分方程，其一般的方程表達式為：其中，X(t)是狀態(tài)變量，μ和σ是狀態(tài)變量的漂移系數(shù)和擴散系數(shù)，B(t)是標準的Brown運動。在描述光源的輸出特性時，X(t)就表示光源的光功率，μX(t)表示光功率的隨機漂移特性，σdB(t)表示光功率的噪聲特性。

由式(6)可以看出，Langevin模型中變量X(t)的輸出特性同X(t)隨時間的變化率相關（在工程上也可理解為同變量的前一時刻X(t- 1 )的特性相關），這是同隨機維納過程模型的最大差別。

同光源的輸出特性模型相比，干涉式閉環(huán)光纖陀螺的輸出特性較復雜，影響因素較多，尚沒有統(tǒng)一的動力學模型。但從光纖陀螺數(shù)理統(tǒng)計分析和可靠性研究角度，參照文獻[3]所提到的研究方法，可以對陀螺的輸出狀態(tài)進行解析，選擇陀螺靜態(tài)條件下的輸出特性進行研究，以反映陀螺自身特性變化引起的輸出性能變化。本文選擇在常溫環(huán)境下陀螺的零位輸出狀態(tài)特性進行研究。

干涉式閉環(huán)光纖陀螺在工程實現(xiàn)上，是利用當前時刻和前一時刻的光源干涉強度差來實現(xiàn)閉環(huán)和信號處理輸出的。借鑒光源輸出的動力性模型，我們將陀螺的零位輸出特性歸納為Langevin方程模型，是符合邏輯和工程特性的，但模型的具體符合度還需要利用數(shù)據(jù)進行檢驗和驗證。因此本文以Langevin方程為基礎，以光纖陀螺理論零位輸出為狀態(tài)變量，構建其狀態(tài)方程。以狀態(tài)方程為基礎，下一步的重點是找到其伴隨方程，完成正倒向隨機微分方程組的建立。

構建倒向隨機微分方程，尚沒有工程模型可以直接借鑒引用。按照倒向隨機微分方程理論，生成元函數(shù)g(X(t),Y(t),z(t),t)的構成形式?jīng)Q定了倒向微分方程的模型。但在理論上，非線性方程的生成元結構非常復雜，沒有通用的模型可直接引用，因此若要建立光纖陀螺的倒向微分方程，只能從方程的機理和工程實際過程入手，構建相應模型并對模型進行檢驗驗證。

由生成元g(X(t),Y(t),z(t),t)特性可知，Y(t)和z(t)是同倒向微分方程的解相關聯(lián)的兩個隨機過程，z(t)是Y(t)波動特性的反應，X(t)是狀態(tài)變量，而Y(t)是可觀測的。另外我們建立光纖陀螺倒向微分方程的目的，就是在輸出波動的情況下（波動是由于陀螺應力變化、環(huán)境影響及測試干擾等造成的），利用未來的陀螺零位值確定當前的零位值。因此將陀螺實際零位Y（t）作為方程的輸出參量，將生成元g(X(t),Y(t),z(t),t)線性化，設為陀螺狀態(tài)變量X(t)的線性函數(shù)，即Y(t) =atX(t) +bt，（其中at為權重系數(shù)，bt為同步偏差），就可完成方程的設立。

通過上述設定，其正倒向隨機方程組的數(shù)學模型為：

其中，X(t)為光纖陀螺的理論零位，是方程組的狀態(tài)變量；Y(t)是光纖陀螺的實際零位，是輸出變量；μ和σ是狀態(tài)變量的漂移系數(shù)和擴散系數(shù)，μX(t)是陀螺零位漂移的趨勢項，σdB(t)是陀螺零位的擴散項，二者都反映了陀螺零位變化的固有特征；z(t)是陀螺的波動特性反映，由陀螺的應力變化、環(huán)境影響及測試干擾等導致。

方程組在工程上的含義為在隨機波動情況下陀螺實際零位同陀螺理論零位之間的函數(shù)關系特性。

2.2 模型求解

根據(jù)模型(7)，光纖陀螺的零位輸出變量Y(t)和狀態(tài)變量X(t)滿足線性關系Y(t) =atX(t) +bt（其中at、bt為待求系數(shù)，設初始值a1= 1 ，b1= 0 ），且終端條件Y(t)=ω已知。為便于方程推導，不失一般性，對時間尺度歸一化，設T= 1 。

當t∈[0 ,T]，利用It?公式，可得出Y(t)關于X(t)的表達式[8]。

利用式(8)(9)，將對方程組的求解轉化為對下列隨機微分方程的參數(shù)μ和σ進行估計。

2.3 模型參數(shù)估計

對模型參數(shù)估計，需要獲得產(chǎn)品在一定時間內(nèi)的測試數(shù)據(jù)。對陀螺采樣測試時間進行歸一化，得到在時刻陀螺的測試采樣數(shù)據(jù)，設0 ≤s≤t≤1，

利用It?公式，可推導出給定tY和sY的條件分布為[8]：

式中，N（·）是個標準的正態(tài)分布。

對所有1 ≤k≤n，將陀螺輸出采樣數(shù)據(jù)簡化表示為用下列簡化符號表示：

本文采用最大似然估計法對參數(shù)進行估計。根據(jù)式(15)，關于觀測值的似然函數(shù)表達式為：

其中，轉移密度函數(shù)為：

文獻[8]證實了利用最大似然估計法所獲得的正倒向隨機微分方程估計參數(shù)具有強收斂性。在工程上，利用采樣數(shù)據(jù)，就可根據(jù)最大似然估計法獲得參數(shù)最優(yōu)估計，并可利用估計參數(shù)計算模型的結果，同實際數(shù)據(jù)相比對，用來驗證模型的正確性。

3 數(shù)值計算和應用

3.1 光纖陀螺采樣數(shù)據(jù)獲取

為了驗證光纖陀螺在衛(wèi)星上長期通電的參數(shù)穩(wěn)定性，我們在長壽命試驗室抽取同一批次的若干只陀螺進行了長期穩(wěn)定性測試試驗，并按照測試規(guī)范，逐月對陀螺的零位進行測試，表1是其中一只陀螺近9年的零位穩(wěn)定性測試數(shù)據(jù)匯總。

為驗證光纖陀螺正倒向隨機微分方程模型的正確性。在數(shù)據(jù)處理上，先利用 2009～2014年的數(shù)據(jù)對模型參數(shù)進行估計，然后再用2015～2019年數(shù)據(jù)對模型參數(shù)進行交叉驗證。

表1 陀螺零位測試數(shù)據(jù)表Tab.1 FOG bias data

3.2 參數(shù)計算

通過3.1節(jié)分析，將正倒向隨機微分方程組的參數(shù)估計問題轉換為一個隨機微分方程的參數(shù)估計問題，進而根據(jù)最大似然估計方法，得出參數(shù)的一致性估計。參數(shù)估計問題轉化為對如式(20)所示方程求根：

對式(20)展開可得：

利用 Matlab的求根數(shù)值解算方法，將表1中的2009～2014年陀螺逐月零位測試數(shù)據(jù)代入方程中，獲得最優(yōu)的參數(shù)估計值為= 0.6532 ，=0.9755。

將參數(shù)值代入方程式(10)中，就得到了光纖陀螺零位輸出變量tY的線性隨機微分方程。為了驗證模型的適宜性和參數(shù)的穩(wěn)定性，利用已知方程可求得模型的計算值，即光纖陀螺零位的估計值。將估計值與實際測量值相對比，可對模型的適宜性進行分析判斷。

模型和參數(shù)值驗證估計效果如圖1所示。圖中藍色點表示陀螺零位的實際測量值，紅色點表示利用方程推導出的陀螺零位估計值。經(jīng)計算擬合誤差的均方根值為0.0481 (°)/h。從驗證的效果看，其擬合優(yōu)度可決系數(shù)R2=0.7033，擬合度較好，在一定程度上可以用式(10)的方程模型來刻畫陀螺零位的變化。

圖1 模型估計值和同實測值對比圖Fig.1 Comparison of model estimates and measured value

3.3 應用示例

利用式(10)得到的光纖陀螺正倒向隨機微分方程模型，可用于求解在給定將來時刻目標值的情況下，如何得到當前（產(chǎn)品出廠）時刻的允許值。

以實際應用為例，用戶要求光纖陀螺在衛(wèi)星上連續(xù)工作 10 年后的零偏不大于 4.0 (°)/h（即ω=4.0 (°)/h），為了確保陀螺滿足衛(wèi)星壽命要求，需確定在出廠時刻（即T=0）允許的零位數(shù)值最大值。

求解光纖陀螺倒向隨機微分方程，可獲得所需結果。由于倒向隨機微分方程含有隨機過程，相關解法不如常規(guī)微分方程的解法成熟，為此相關學者對倒向微分方程的數(shù)值解法進行了大量研究[9]，主要從純數(shù)學角度研究如何提高算法精度，其中涉及到大量的數(shù)學條件期望求解算法等。

本文重點關注模型的建立和工程應用，數(shù)值精度不是重點，因此利用最基本的歐拉數(shù)值解法，通過對布朗運動進行蒙特卡羅方法模擬，將參數(shù)（,）代入進行求解。

蒙特卡羅方法對隨機過程的布朗運動抽樣生成的一般方法為：

● 采樣點 0 =t0＜t1＜t2＜…＜tn；

● 標準正態(tài)分布 (0,1)N抽樣產(chǎn)生

k= 1 ,2,… ,n。

經(jīng)過對微分方程的歐拉數(shù)值解法計算，在陀螺10年后的零位不大于4.0 (°)/h（即ω=4.0 (°)/h）情況下，陀螺出廠零位應不大于1.5 (°)/h。計算獲得的陀螺連續(xù)工作10年（120個月）的零位變化如圖2所示。

圖2 估計的陀螺零位變化圖Fig2 Variation curve of FOG estimate bias

為了清楚得到陀螺的零位變化趨勢，在數(shù)據(jù)處理上進行了平滑處理。結果可以看出，在終端時刻的要求不同，初始時刻的取值不同。

4 結論和展望

本文從光纖陀螺零位隨機漂移的規(guī)律和模型出發(fā)，建立光纖陀螺的倒向隨機微分方程模型，并利用最大似然估計法求解了模型參數(shù)，并對模型和參數(shù)進行了分析和驗證，最后以某光纖陀螺產(chǎn)品實際要求為例，獲得了滿足用戶零位穩(wěn)定性要求的出廠控制指標。

由于倒向隨機微分方程在工程領域的實際應用尚處于探索階段，本文也是首次將倒向隨機微分方程應用于光纖陀螺的質量控制領域，研究子樣還不夠多，后續(xù)還需要對模型的適應性進一步進行驗證，但該方法對拓展光纖陀螺乃至慣性產(chǎn)品質量管控的研究視野必將起到積極的作用。