基于解析模型的小樣本操縱性KT指數辨識方法

張煒靈，蔡 烽，王 驍

(海軍大連艦艇學院， 遼寧 大連 116018)

自從野本謙作提出利用Z形試驗求取操縱響應KT指數的方法之后，KT方程在船舶操縱建模中得到廣泛應用。隨著計算機的發展，為提升系統辨識率，對KT指數的辨識從最原始的作圖法逐步演變。[1]鄭文龍等[2]通過自然回歸離散模型的最小二乘法建立一階和二階非線性KT連續模型；張彬等[3]利用分步辨識的方法對一階和二階非線性KT方程的參數進行辨識；SERGE等[4]采用經典的遺傳算法結合自航模試驗辨識船舶運動模型的水動力導數;YIN等[5]運用偏最小二乘法辨識整體型船舶操縱運動數學模型水動力導數；謝朔等[6-7]利用新息最小二乘法和卡爾曼濾波辨識KT方程參數。

就目前而言，幾乎所有的系統辨識方法的目標函數都為式(1)和式(2)，是基于微分形式的誤差平方和。該目標函數雖然表面上可使數據與模型相匹配，但本質上簡單地將角速度、角加速度和角加加速度作為獨立變量代入，忽略了相互之間的關聯，其處理顯然不夠完善。

(1)

(2)

本文在研究KT操縱響應模型的基礎上，求解在連續線性操舵情況下一階和二階操縱響應解析模型。直接利用航向角數據構建目標函數，采用基于格雷碼和精英選擇的遺傳算法進行非線性尋優。通過與最小二乘法相比較，證明本文提出的算法能在“理論上最優”的基礎上將航向角平均偏差降低約40%。為加快辨識效率，在辨識過程中調整航向角的時間步長，進行多次辨識驗證，說明該方法只需要小樣本數據就可進行高效準確的辨識，極大地提升辨識效率。

1 基于解析模型構建優化函數

在實際操縱過程中，假設舵角是線性變化的，即某次操縱過程中舵角變化率恒定。

操舵情況示意見圖1。假設操舵的時間向量為ti(i=1,2,…)，對應的舵角值為δi(i=1,2,…)，則對于t∈[ti,ti+1)，有

(3)

圖1 操舵情況示意

1.1 一階KT響應方程解析模型

一階KT操縱響應方程[8]為

(4)

將其適當變化為

(5)

φ(t)=F2(t)-F2(ti)+Ni-1+φ0，t∈[ti,ti+1)

(6)

(7)

(8)

(9)

(10)

式(6)～式(10)中：N0=0，M0=0。

1.2 二階KT響應方程解析模型

二階KT操縱響應方程為

(11)

將其適當變化為

(12)

(13)

(14)

對r進行分步積分后，即可得航向角為

φ(t)=F3(t)-F3(ti)+Li-1+φ0，t∈[ti,ti+1)

(15)

(16)

(17)

式(15)～式(17)中:L0=0。

1.3 構建優化函數

對于某次試驗而言，操作過程已確定，即ti(i=1,2,…)和對應的δi(i=1,2,…)都已確定。根據給出的操縱響應解析模型為

φ1(t)=F2(t)-F2(ti)+Ni-1+φ0

(18)

φ2(t)=F3(t)-F3(ti)+Li-1+φ0

(19)

式(18)和式(19)中：φ1(t)為一階解析模型；φ2(t)為二階解析模型。

在給定ti和δi之后，需再確定KT指數，航向角就能被確定。將解析模型改寫成關于KT指數的函數為

(20)

(21)

利用解析的航向數據和試驗數據的誤差平方和作為KT指數優化的評價依據，為將其與輪盤賭相結合，實現偏差越小，適應度越強，個體選擇的概率越大的目標，將偏差的倒數作為優化的目標函數，優化目標越大越好。

(22)

(23)

2 基于遺傳算法進行非線性尋優

2.1 格雷碼編碼

格雷碼的一位等位基因變化只會引起十進制的微小變化，使遺傳算法的交叉變異操作的尋優效率大大提升。為既方便染色體的交叉變異操作，又能提高遺傳算法的尋優效率，采用格雷碼的編碼方式。[9]

假設有一個二進制編碼為

B=bmbm-1…b2b1

(24)

式(24)對應的格雷碼為

G=gmgm-1…g2g1

(25)

格雷碼和二進制碼轉換為

(26)

式(26)中：“⊕”為異或運算。

在辨識過程中，每個參數的編碼長度為l，每組參數的基因組合拼接成1條染色體，其編碼長度為

L=n·l

(27)

式(27)中：n為需要優化的參數數量。

由于各參數的取值范圍各不相同，為能準確辨識，根據艦船的基本情況，設定合理的取值范圍[a，b]。在解碼時，先將格雷碼解碼成二進制編碼，在二進制編碼轉換成十進制編碼時，進行伸縮變換，將值域控制在各參數的取值范圍內，即

(28)

(29)

式(28)和式(29)中：x為從二進制轉換成的十進制形式參數；x′為經過伸縮變換后的十進制參數；A和B為x的下限和上限；a和b為x′的下限和上限。

2.2 精英選擇

用比例算子對種群進行選擇并對當代最優個體進行精英保留。首先，將上一代最優個體計入當代個體中，計算每一個體的適應度，再標記此時的最優個體，同時剔除當代適應度最低的個體，以維持種群數量。按照個體占種群總體適應度的比例進行歸一化，將其作為個體選擇復制到下一代的概率。

3 辨識結果比較

在進行交叉和變異時，對指定參數的染色體區域進行操作，不同的參數之間不進行交叉和變異操作，其中：交叉概率為pc;變異概率為pm。

在進行辨識時，種群數量M=50，交叉概率pc=0.4，變異概率pm=0.01，編碼長度為l=15。從二進制轉換為十進制的范圍為0～(215-1)。

1) 一階辨識時回轉性指數K設為0～0.1;應舵性指數T設為0～20;壓舵舵角δ0設為-2～2。

2) 二階辨識時K為0～0.1;T1設為0～20;T2設為0～10;T3設為0～10;壓舵舵角δ0設為-2～2;時間步長設為1 s;遺傳代數設置為200代。

在長山群島西南海域對某艦進行多次Z形試驗,以其中一次10°/10°Z形試驗為數據基礎，進行操縱性KT指數辨識。試驗的操舵t和對應δ為

(30)

隨機生成數量為50的初始種群，種群對應的參數代入到解析模型中進行仿真，對應50條隨機的航向角曲線，見圖2和圖3。

經過遺傳算法的優化，種群進化200代之后，曲線幾乎集中在同一區域，表明種群的適應性已到達最優。航向角曲線見圖4和圖5。

圖4和圖5在一定程度上說明了辨識結果的收斂性，200代種群的適應性已普遍很高。辨識過程的收斂情況見圖6和圖7。

一階辨識得到的最優參數為：K=0.075 9;T=17.487 8;δ0=-0.772 3;航向角平均偏差為0.491 1°，最大偏差為1.289 5°。二階辨識得到的最優參數為：K=0.072 2;T1=14.398 2;T2=4.067 1;T3=2.115 5;δ0=-0.784 9;航向角平均偏差為0.424 8°;最大偏差為1.180 9°。最優航向角曲線見圖8和圖9。

圖6 一階辨識優化趨勢 圖7 二階辨識優化趨勢

圖8 一階最優航向角曲線 圖9 二階最優航向角曲線

利用最小二乘法對式(1)和式(2)進行參數辨識。將最小二乘法的辨識結果與本文的辨識結果相比較，結果見表1。本文的辨識方法能有效降低航向角的最大偏差，并將平均偏差縮小約40%，有效提高辨識的準確性。同時，由于在辨識之前已設置參數的值域范圍，在實數域中進行尋優，能有效避免最小二乘估計時出現參數為復數和負數的情況。

表1 辨識結果和偏差比較表

4 小樣本辨識分析

在提出以上基于解析模型的辨識方法之后，只要在辨識過程中仿真的航向角數據與試驗的航向角數據對應的時刻是相同的，目標函數內所含的航向角數據在理論上可取任意值，即在辨識過程中完全可幾秒取1個數據，甚至幾十秒取1個數據進行辨識。步長為30 s的初始曲線和優化后的曲線見圖10和圖11。

基于以上試驗數據，取不同的時間步長進行多次辨識，以最終辨識得到的KT指數對應的航向角平均偏差作為辨識好壞的衡量標準，見表2。航向角均差隨時間步長的變化趨勢見圖12。

表2 不同樣本數量辨識的偏差情況

圖12 航向角偏差趨勢

由表2和圖12可知：辨識的時間步長減小到約50 s后，即當辨識數據樣本大于等于6個之后，辨識得到的KT指數對應的航向角平均偏差逐漸開始收斂。當時間步長縮小到30 s，即數據樣本增加到10個之后，樣本數量增加并不能大幅度提高辨識精度。在利用該方法進行辨識時，時間步長只需取30 s，即只需要10個航向角數據點就能得到理想的效果，且其結果明顯優于最小二乘法。

5 不確定度分析

考慮到在試驗過程中各種數據都存在一定的誤差，這些誤差最終都會通過辨識過程傳遞到KT指數。為衡量最終辨識得到的KT指數的可靠程度，引入不確定度的概念定量分析KT指數的數據質量。

5.1 不確定度元

5.1.1航向角不確定度

參照艦船陀螺羅經的設備指標說明書，在動基座羅經狀態下，勻速直航航向偏差εφ1≤±0.8°×secφ。艦船所處維度為38°50′，主羅經的偏差εφ1≤±1.027 0°。主羅經傳遞到分羅經的傳遞偏差εφ2≤0.2°。假設這2種偏差在偏差極限內服從正態分布，根據國標中不確定度計算原理[10]:主羅經的標準不確定度μφ1=0.342 3°，自由度νφ1=50；傳遞偏差引起的不確定度μφ2=0.066 7°，自由度νφ2=50。航向角的合成不確定度為

(31)

根據韋爾奇-薩特思韋特公式計算得到航向角的有效不確定度的自由度為

(32)

航向角數據的有效不確定度對應的有效自由度取νφ=50。

5.1.2 舵角不確定度

自動操舵儀的操舵盤最小讀數為1.0°，讀數偏差εδ1≤±0.5°。自動舵的靈敏度為0.5°，由其引起的操舵舵角的偏差εδ2≤±0.5°。假設讀數偏差服從均勻分布，靈敏度偏差服從正態分布，則根據不確定度計算原理，讀數偏差引起的標準不確定度為μδ1=0.288 7°，自由度為νδ1=12；靈敏度偏差引起的標準不確定為μδ2=0.166 7°，自由度為νδ2=50。δ的合成不確定度為

(33)

由式(33)計算得到舵角的有效不確定度的自由度為

(34)

δ的有效不確定度對應的有效自由度取νδ=20。

5.2 不確定度合成

將基于解析解辨識一階操縱性KT指數的過程也看成函數，即

[K,T]=f1(φ1,…，φm,δ1,…，δn)

(35)

[K,T1,T2,T3]=f2(φ1,…,φm,δ1,…,δn)

(36)

根據國標中的不確定度傳遞原理[10]，將航向角數據的不確定度和舵角數據的不確定度代入基于解析模型辨識一階和二階KT指數的過程中，得到KT指數的標準不確定度見表3。

表3 KT指數不確定度結果

由于本文提出的辨識方法只用到航向角數據，少量的航向角數據的變化雖然會引起目標函數值的變化，但對其最優參數的影響很小。相比所有誤差都會被計算在內的最小二乘法，利用提出的辨識方法能有效降低辨識得到的KT指數的不確定度。由于該辨識方法得到的T2和T3的數值相對較大，其相對不確定度明顯下降。在樣本存在一定不確定度的前提下，有效證明了該辨識方法的可靠性。

6 結束語