基于最大似然估計和混合梯度優化的射手模型辨識

吳駿雄， 林德福， 王輝， 袁亦方

(1.北京理工大學 無人機自主控制技術北京市重點實驗室, 北京 100081；2.北京特種機電研究所, 北京 100012)

0 引言

光纖圖像制導武器使用光纖圖像與指令雙向傳輸體制，可在復雜的交戰條件下攻擊坦克、艦船、地面或低速空中目標。彈載圖像導引頭能夠將圖像探測器敏感的戰場實時圖像傳輸到控制地面站，射手可根據顯示器中的圖像選擇攻擊目標以及攻擊部位，并通過圖像跟蹤器自動跟蹤，也可以通過控制手柄進行手動跟蹤直至命中目標，即所謂的“人在回路”制導方式[1]。因此，了解射手對導引頭跟蹤誤差的響應并對射手本身進行建模，不僅可以更深入地了解射手的感知過程，而且能夠對制導控制系統設計[2]以及培訓模擬設備評價起重要作用[3]。

對人進行建模的研究始于20世紀50年代[4]，其中最具代表性的是McRuer等提出的精確模型，該模型已在諸多領域應用[5]。隨著現代控制理論的發展，Kleinman等[6]建立了人的最優控制模型，Hess[7]提出了能充分體現人的信號處理過程的Hess結構模型。近來，將智能技術應用到人的建模問題中逐漸成為研究熱點，其中模糊控制和神經網絡技術已經成功應用到飛行員模型辨識中[8-9]。

人響應的具體特征是由模型參數決定的[10]，可以利用參數辨識技術從實驗數據中得到。目前人的模型參數辨識方法可分為兩類：一類是針對人對時不變控制對象的響應；另一類是針對人對環境或者控制對象動力學突然變化時的響應[11]。傅里葉系數[12]和最大似然估計[13]是辨識時不變控制對象人的行為典型方法，前者需先估計非參數化的系統頻率響應，再通過擬合非參數化模型得到模型參數；后者直接通過擬合時域數據來辨識模型參數。與傅里葉系數法相比，最大似然估計方法僅使用1個步驟來辨識模型，減小了參數估計的偏差和方差。針對時變控制對象人的響應辨識問題主要有最大似然、小波變換[11]和Kalman濾波[14]等方法。“人在回路”制導體制下射手的控制對象是導引頭內回路，在飛行過程中其動力學特性幾乎不變。因此，本文采用最大似然估計方法來辨識射手模型。

基于上述分析，本文擬針對射手控制行為進行建模辨識。目前國內外關于射手模型的文獻較少，未有詳細的關于射手建模的相關文獻。其中趙軍民等采用輔助變量法對模型參數進行辨識[2]，但其認為模型的神經延遲環節固定，實驗條件也僅考慮在跟蹤橫向運動的坦克目標，未能充分激勵射手的響應。因此，本文將設計具有較寬頻率范圍的激勵信號，而且為了提高辨識精度，采用輸出誤差結構的最大似然法。此外，針對參數估計過程梯度優化算法對初值敏感且易發散的特性，采用混合梯度算法進行優化，即利用遺傳算法逼近最優值，利用梯度算法的快速性加快算法運行，并使用單純形法保證算法的穩定性。

1 參數估計問題

1.1 控制任務

光纖圖像制導系統采用大容量光纖傳輸系統，能夠將電視導引頭的圖像實時傳輸到地面站，且能在武器地面站實現傳統的導引頭圖像跟蹤能力，使得跟蹤器可以重復使用。相比于彈載圖像跟蹤器，地面站圖像跟蹤器的功能可以設計得更加強大和全面，同時降低武器使用成本。人工控制模式下導引頭的控制框圖如圖1所示。

圖1中：qt為目標的慣性空間視線角；qs為導引頭光軸相對于慣性空間的角度；Δq為視場角誤差信號；Hp(s)為射手模型的傳遞函數，s為拉普拉斯算子；n為射手模型中的殘差信號；u為射手的輸出；K為信號放大器；G(s)和H(s)分別為穩定回路校正網絡和速率陀螺。當目標偏離光軸時，射手通過顯示屏觀測到視場角誤差Δq，并操控手柄使得光軸指向目標，實現對目標的跟蹤。射手在導引頭跟蹤回路中相當于1個反饋控制器，輸入信號為屏幕上顯示的視場角誤差。因此射手的目標跟蹤任務可認為是單回路補償任務，控制對象即為圖1中的虛線框內部分。不失一般性，可假設G(s)≈1，H(s)≈1，且導引頭穩定回路的帶寬約為20～30 Hz，穩定回路的傳遞函數可近似為100/(s+100)[15]，則控制對象的傳遞函數Hc(s)可表示如下：

(1)

(1)式中的增益K決定了跟蹤回路的帶寬和穩定裕度。導引頭跟蹤回路帶寬遠小于穩定回路帶寬，為保證導引頭系統回路具有足夠的穩定裕度，在這里選取增益K=2.

1.2 射手模型

本文采用精確模型[16]來描述射手行為。如圖1所示，射手模型的響應包含傳遞函數Hp(s)的響應和殘差信號n兩部分，其中殘差信號代表射手模型的非線性部分。射手對顯示運動提示作出響應的線性部分形式如下：

(2)

式中：Kp為視覺感知的增益；τL為視覺感知超前時間常數；τl為視覺感知滯后時間常數；ωnm為系統自然頻率；ξnm為系統阻尼；τd為視覺感知時間延時。 (2)式表明神經肌肉系統的動力學限制了射手響應，神經肌肉系統認為是1個2階質量彈簧阻尼系統。

McRuer的交叉定理[16]表明，人會調整自身的配平動力學來適應控制對象的動力學，使得系統開環傳遞函數在交叉頻率附近具有類似單積分器系統的頻率特性。對于(1)式所示控制對象的動力學，射手需要在更高的頻率處進行配平[17]。因此圖1中射手傳遞函數表示如下：

(3)

綜上所述，共有5個參數需要被估計，參數向量θ=[Kp,τL,τd,ξnm,ωnm]T。

2 射手模型辨識

下面基于最大似然估計方法來得到射手模型的參數，為簡化整個參數估計過程，采用輸出誤差的結構來構造似然函數。為了解決非線性優化問題的全局最優問題，并避免出現局部最優現象，即找到觀測值出現的最大概率，在使用經典的梯度優化算法之前，使用遺傳算法來決定梯度算法的初值。

2.1 最大似然估計方法

最大似然估計是一種時域參數估計方法，在眾多領域都有著應用[18]。最大似然估計是一類概率性貝葉斯估計方法，通過引入觀測量的條件概率密度，構造1個以觀測數據和未知參數為自變量的似然函數，以觀測值出現的概率最大作為估計準則，獲得系統模型的參數估計值。

在應用最大似然估計方法時，需要將射手傳遞函數模型表示為狀態空間形式。在轉化過程中，射手模型中的純延時環節為非線性環節，可以將延時環節近似為1個高階的傳遞函數，其中1種近似方法為Padé近似[19]，其具體形式如下：

(4)

式中：T為延遲時間；ka=0，…，λ，λ為Padé近似的階次。為保證延時模型在高頻處的精確描述，且不使近似射手模型階次過高，本文采用5階延時環節Padé近似[13]：

(5)

將(5)式代入(3)式中，可得近似的射手傳遞函數為

(6)

式中：b0，…，b6和a0，…，a6分別為對應的傳遞函數系數。

由于引入了Padé近似，射手傳遞函數的系數為未知辨識參數的非線性函數，例如(6)式中分母的第1項和最后1項系數為

(7)

將射手傳遞函數(6)式轉換為狀態空間的標準型，可得

(8)

式中：x為狀態向量；狀態矩陣A(θ)、控制矩陣B(θ)和輸出矩陣C(θ)的表達式分別為

(9)

(8)式給出的射手模型僅在輸出方程中包含1個額外的噪聲項，即假設沒有過程噪聲的影響，從而可以采用輸出誤差結構的最大似然估計，大大簡化了辨識步驟。在文獻[20]關于辨識人的模型研究中，一般假設人的響應殘差部分是1個零均值高斯噪聲，并通過實驗驗證該假設。本文同樣假設射手殘差具有相同的特性，即認為(8)式的模型中殘差是高斯量測噪聲，其均值和方差定義如下：

(10)

式中：n(k)為殘差信號n的離散采樣時間序列，k為序列號，k=1,…,m，m為采樣點個數。

最大似然估計方法中，需要尋找模型參數向量θ的估計值，該估計值使得似然函數L(θ)具有最大值。L(θ)定義為預測誤差的條件概率密度函數：

L(θ)=f(v(1)，v(2)，…，v(k)，…，v(m)|θ),

(11)

式中：v(k)為預測誤差，定義為量測的離散射手輸出信號u(k)和離散射手模型輸出信號(k)之差，

v(k)=u(k)-(k).

(12)

由(10)式定義的殘差性質，可得似然函數為

(13)

圖2所示為輸出誤差辨識策略的基本框架，該策略需要最小化1個二次型罰函數，該函數與輸出誤差相關，其中輸出誤差為射手實際輸出和在相同輸入情況下辨識模型的仿真輸出之差。因此，罰函數定義為實際數據和仿真數據的匹配度，即參數估計值和實際估計值之間的差別。輸出誤差辨識策略采用迭代搜索算法在參數空間中尋找最優參數值，使得罰函數值最小，則最優參數值即為最終的參數估計值。

在應用最大似然估計方法時，若估計的參數使得似然函數最大，則為模型參數向量θ的估計值；而在輸出誤差策略下使用最大似然估計方法，需要用求取似然函數對數形式的最小值代替求取似然函數的最大值，從而簡化為一個更加直接的優化問題。當對數似然取得全局最小時，得到的參數向量即為最大似然估計值，表示為ML. 因此，輸出誤差策略下的最大似然估計表示為

(14)

(14)式概括了本文射手模型的參數估計問題，它定義了1個強非線性優化問題。下面采用混合梯度優化算法來得到該優化問題的全局最優解。

2.2 混合梯度優化算法

傳統求解最大似然極小值的方法有Levenberg-Marquardt、高斯- 牛頓和Newton-Raphson等方法。本文提出一種混合梯度優化算法，采用遺傳算法和高斯- 牛頓梯度優化算法相結合的方式來估計射手傳遞函數的參數，具體策略如圖3所示。

在沒有先驗信息情況下，該混合優化算法首先應用遺傳算法進行初始迭代，然后用高斯- 牛頓和單純形法算法來優化(14)式。由于較大的初始誤差可能會導致參數估計發散，在計算初始時刻使用遺傳算法以減小算法對初始值的敏感。此外，高斯- 牛頓算法具有較快的收斂速度，但當離真實解較遠時會產生不準確的梯度信息、導致發散。因此，為增強算法的魯棒性，當高斯- 牛頓算法發散時，切換為Nelder-Mead單純形法進行迭代，直至高斯- 牛頓算法能夠收斂。

2.2.1 遺傳算法

遺傳算法固有的隨機性以及在初始參數向量集內的最優解開始搜索，導致該算法有很高的概率找到優化問題的全局最優值[21]。遺傳算法首先需要創建1個初始種群，即初始參數向量集，該種群在參數的上下界隨機選取，在此用矩陣表示為

(15)

通過目標函數(14)式評估種群中每個個體的適應度。圖4所示為對當前種群每個個體運用3個基因操作從而產生新種群的過程，具體步驟如下：

1) 復制(elite)：從當前種群中選擇適應度強的個體，并復制給下一代；

2) 交叉(crossover)：從當前種群中任選2個個體，隨機選擇基因中交叉點的位置，模擬進化過程中的繁殖現象，得到全新基因的個體；

3) 變異(mutation)：當前個體基因小概率隨機選擇變異點，并發生變化從而產生新的個體。

不斷重復上述步驟，直至終止條件滿足。

由于遺傳算法使用概率搜索技術，導致1個搜索點到另1個搜索點的轉移方式和關系具有不確定性，也導致每次迭代出的結果不一致。因此，從測試數據中來估計射手模型的參數，并不適合僅僅使用遺傳算法。

2.2.2 無約束高斯- 牛頓優化算法

遺傳算法迭代得到的模型參數估計值，在很大程度上接近于優化問題的全局最優解。為進一步得到精確的參數估計值，將遺傳算法的解作為高斯- 牛頓優化算法的初始參數估計值，這種基于梯度的優化算法是解決最大似然優化問題的一種經典求解方法[22]。高斯- 牛頓優化算法的參數迭代更新方程由下式給出：

(i+1)=(i)+Δ,Δ

(16)

(17)

Mθ=θ(i)為Fisher信息矩陣，

(18)

(16)式中需要求解Fisher信息矩陣的逆矩陣，則要求Fisher信息矩陣可逆。通過求解Fisher信息矩陣可以得到Cramér-Rao的下界(CRLB)[23]，即參數估計可實現的最小方差，表示為Jc. 當未知模型參數過多時或數據中信息量不足時，信息矩陣M容易變成病態矩陣甚至奇異矩陣，導致迭代步長在某些方向上過大或者信息矩陣不可逆。因此，將信息矩陣M進行奇異值分解并求逆，有

(19)

式中：tks、uks分別為矩陣T和U的第ks列；δ1，…，δr為信息矩陣的r個奇異值。將奇異值從大到小排列，有

δ1≥δ2≥…≥δks≥δr.

(20)

一般而言，信息矩陣中最小的幾個奇異值所包含的矩陣信息較少，且當奇異值過小時容易導致信息矩陣病態甚至不可逆，故可將最小的幾個奇異值舍去。假設有η個奇異值滿足舍去的標準：

(21)

式中：ε為計算機的計算精度；δmax為最大奇異值。經過上述處理后，信息矩陣的逆矩陣階次從r階降低到r-η階，故該方法又稱為降階矩陣求逆法[24]。

此外，為防止高斯- 牛頓優化算法發散而終止迭代，將該算法切換成Nelder-Mead單純形法。Nelder-Mead單純形法是一種多維直接搜索的局部優化方法，在尋優過程中不必計算目標函數的梯度，只是針對一定圖形的頂點，按照一定規則進行搜索，該方法操作簡單且一直是收斂的，但是計算速度較慢，具體過程詳見文獻[25]。

(22)

3 仿真實驗及結果分析

3.1 激勵信號

下面介紹射手目標跟蹤補償實驗，目標視場角信號驅動目標在視場內運動，射手根據視場誤差角控制手柄，使光軸對準目標。為了防止射手預測信號的趨勢，激勵信號應表現一定的隨機性，McRuer等[16]在實驗中，將激勵信號設計為由若干正弦信號疊加而成，其形式為

(23)

式中：Akf為正弦信號幅值；ωkf為正弦信號頻率；φkf為正弦信號相位；kf=1,…,N，N為正弦信號個數。

為了得到較高精度的射手動力學模型，激勵信號必須具有充分寬的頻帶和高的信噪比，激勵信號的每個部分即Akf、ωkf和φkf都需要精確設計。因此將激勵信號的頻譜設計成1階低通濾波器形狀[12]，每個頻率點的幅值為

(24)

表1 激勵信號參數

3.2 實驗裝置和實驗過程

實驗在實驗室開發的仿真系統上進行，該系統主要包括顯示子系統、動力學仿真子系統和控制手柄，仿真系統原理圖如圖5所示。其中動力學仿真子系統在嵌入式系統RTX基礎上開發，用以采集手柄的輸出并實時仿真導引頭系統動力學，計算頻率為1 000 Hz；顯示子系統采用vista2D+微軟基礎類庫(MFC)方式開發，用以顯示視場角誤差，其顯示的信息通過用戶數據報協議(UDP)和動力學仿真子系統通信得到，然后射手根據顯示信息操縱手柄，將導引頭光軸對準目標以減小跟蹤誤差。

實驗中共有5個參試人員，在采集數據前，所有參試人員需進行相當次數的訓練，直到他們的跟蹤表現穩定在一定水平。參試人員的跟蹤水平可以通過誤差信號均方根來評價，每次訓練后參試人員都會被告知他們的跟蹤評分，以便在初始熟悉期間激勵改善他們的跟蹤表現，且在達到一定熟練程度后用以保持穩定的表現。為了提高辨識結果精度，每個參試人員重復5次實驗并采集相應的測量數據。

3.3 實驗結果

實驗中圖1所示的所有信號都被記錄下來，為了減少信號中噪聲部分的影響并改善射手模型參數估計的精度，將每個參試人員的5次重復實驗信號求平均值。圖6所示為平均視場角誤差信號Δq、射手輸出信號u、目標視場角(激勵信號)qt、光軸相對于慣性系夾角qs. 由圖6可知參試人員成功將光軸對準了目標。

3.3.1 算法性能分析

表2 遺傳算法參數域的上界和下界

遺傳算法是一種基于概率啟發式的搜索算法，在一定迭代步長限制下，每次迭代結果不完全相同。對單個參試人員的實驗結果重復若干次遺傳算法辨識運算，將迭代結果最小的10個解代入高斯- 牛頓優化算法進一步優化。圖7(a)給出了這10個解的全部迭代過程，由圖7(a)可知，重復10次遺傳算法最終給出的值都不同，但通過高斯- 牛頓算法優化后最終都得到了相同的優化結果，即得到全局最小對數似然函數值，其對應的模型參數即為辨識結果。因為兩個優化算法中殘差方差的初始估計值都設置為1，所以得到的目標函數值較大，經過前述交替迭代過程最終得到各自的估計值，具體結果如圖7(b)所示。

為了對比說明混合優化算法的優勢，選擇高斯- 牛頓優化算法結合單純形法的復合優化算法[24]來估計模型參數。在表2參數域內隨機選擇100個初始參數向量，并分別使用高斯- 牛頓優化算法結合單純形法的復合優化算法進行尋優，圖8以升序方式給出了相應的結果。從圖8可知，收斂至全局最小的參數向量小于一半，剩下的參數向量將導致收斂發散或者收斂至局部最小值。由此可見，采用混合優化算法能大大地減小初值選取工作，且能夠精確地收斂至全局最小值。

3.3.2 結果分析

表3 模型參數辨識結果(參試人員1)

表4給出了本次實驗所有參試人員的辨識結果，圖10為相應的頻率響應曲線。從圖10可知，所有參試人員有著一致的辨識結果。其中不同參試人員的傳遞函數略有不同，這是因為每個參試人員有著不同的控制策略且個體之間的差異導致的，這種現象在該類型實驗中普遍存在。

表4 所有參試人員的辨識結果

相關性系數ρu一般用于評價辨識模型的準確性，能夠用來表示辨識模型的輸出再現量測輸出的能力。如果相關系數接近于1，則表示辨識模型的輸出能很好地擬合量測輸出，辨識模型可以較真實地反映實際對象的動態特性；反之，如果相關系數接近于0，則表示辨識結果較差，未能捕捉到實際對象的動態特性。相關系數ρu的定義[26]如下：

(25)

表4中參試人員1的相關系數為86%，可知模型輸出能夠較好地擬合量測輸出，辨識模型能夠反映射手動態特性。圖11給出的參試人員1量測輸出和模型輸出的對比圖也證明了該結論。表5給出了所有參試人員的相關系數。由表5可見，本實驗得到的每個參試人員辨識模型相關性系數均大于70%，可以認為辨識模型能夠較好地擬合射手的動力學特性。

相關系數參試人員12345ρu^u/%86.1872.6980.9585.7492.52

3.3.3 殘差分析

射手輸出除了線性輸出部分，還包含模型為高斯噪聲的殘差，殘差可以由辨識線性射手模型輸出和實際量測射手的輸出之差得到。圖12(a)給出了參試人員1的殘差時域圖，其相應的功率譜密度(PSD)在圖12(b)中給出。由圖12可以看到，殘差信號并不為白噪聲，其特性近似為1個3階低通濾波器：

(26)

式中：Kn為殘差強度，Kn=0.26；ωn為殘差濾波器的轉折頻率，ωn=0.8 rad/s；ξn為殘差濾波器的阻尼系數，ξn=0.06.

射手模型殘差的概率密度函數如圖13所示。由圖13可見：有殘差信號的分布接近于零均值的高斯分布，同時也驗證了之前的假設；此外殘差信號譜特性并沒有影響到信號的正態性。

4 結論

本文針對光纖制導武器中射手模型估計問題，采用基于輸出誤差框架下最大似然估計方法對射手模型參數和殘差特性進行辨識。為提高優化算法的性能和魯棒性，引入混合梯度優化算法，并成功應用于實驗得到的數據中，經過多輪次迭代計算都能夠找到全局最優解。辨識結果表明：算法能夠從時域數據中精確地辨識得到射手模型參數，射手辨識模型能夠真實地反映射手的動態過程；射手的殘差信號基本符合正態分布，且具有色噪聲的頻率特性。