秦偉亮，魏法杰

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

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

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

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

1 倒向隨機(jī)微分方程基本理論

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

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

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

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

引入隨機(jī)因素，即建立正、倒向隨機(jī)微分方程。兩者對(duì)比，其結(jié)構(gòu)形式為：

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

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

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

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

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

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

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

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

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

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

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

2.1 基本模型

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

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

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

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

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

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

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

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

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

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

2.2 模型求解

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

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

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

2.3 模型參數(shù)估計(jì)

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

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

式中，N（·）是個(gè)標(biāo)準(zhǔn)的正態(tài)分布。

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

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

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

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

3 數(shù)值計(jì)算和應(yīng)用

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

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

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

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

3.2 參數(shù)計(jì)算

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

對(duì)式(20)展開可得：

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

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

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

圖1 模型估計(jì)值和同實(shí)測(cè)值對(duì)比圖Fig.1 Comparison of model estimates and measured value

3.3 應(yīng)用示例

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

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

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

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

蒙特卡羅方法對(duì)隨機(jī)過程的布朗運(yùn)動(dòng)抽樣生成的一般方法為：

● 采樣點(diǎn) 0 =t0＜t1＜t2＜…＜tn；

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

k= 1 ,2,… ,n。

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

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

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

4 結(jié)論和展望

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

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