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

秦偉亮，魏法杰

（北京航空航天大學 經濟管理學院，北京 100191）

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

與正向隨機微分方程相比，倒向隨機微分方程在結構上有較大不同，它們的主要區別和聯系表現在：

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

2）正向隨機微分過程關注如何認識一個客觀存在的隨機變化過程，而倒向微分方程關注在隨機干擾條件下如何使一個產品或系統達到預期的要求，也就是通過將來時刻給定的一個目標值來獲得當前（產品出廠）時刻的值。正向隨機微分方程的值代表今天確定的值變為明天一般不確定的狀態以研究其統計規律，而倒向隨機微分方程的值則表示將明天的目標變成今天的解以制定今天的決策[6]。

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

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

2.1 基本模型

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

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

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

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

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

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

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

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

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

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

2.2 模型求解

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

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

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

2.3 模型參數估計

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

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

式中，N（·）是個標準的正態分布。

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

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

其中，轉移密度函數為：

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

3 數值計算和應用

3.1 光纖陀螺采樣數據獲取

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

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

表1 陀螺零位測試數據表Tab.1 FOG bias data

3.2 參數計算

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

對式(20)展開可得：

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

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

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

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

3.3 應用示例

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

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

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

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

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

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

● 標準正態分布 (0,1)N抽樣產生

k= 1 ,2,… ,n。

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

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

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

4 結論和展望

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

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