基于修正輔助變量法的船舶操縱響應模型辨識

曾道輝， 蔡成濤,2,3

(1.哈爾濱工程大學 智能科學與工程學院，黑龍江 哈爾濱 150001； 2.船海裝備智能化技術與應用教育部重點實驗室，黑龍江 哈爾濱 150001； 3.電子政務建模仿真國家工程實驗室，黑龍江 哈爾濱 150001)

隨著船舶智能化和自主化水平的提升，使得船舶在操縱性能變化的動態環境和多任務場景約束下的安全高效航行更依賴于在線建立的船舶操縱運動模型。船舶操縱運動模型在用于辨識的結構方面，可具體分為分離型[1]、整體型[2]、響應型[3]和時間序列[4-5]多種模型結構形式。而船舶操縱響應模型作為船舶操縱運動模型中的一種，其中的操縱性指數物理意義明顯，能夠直觀反映船舶的操縱性能，是一種描述船舶舵角與航向角映射關系的動態響應模型。通過對船舶操縱響應模型的在線辨識建模，不僅能夠為船舶的自動舵和航向自適應控制器提供精確的數學描述，而且能預報船舶的未來航向以及通過辨識的操縱性指數實時反映船舶的操縱性能。因此，實現船舶操縱響應模型參數的在線辨識，對船舶的航向預報和航向控制有著重要意義。

為了辨識船舶操縱響應模型中的參數，目前已經有許多系統辨識方法用于船舶操縱響應模型的在線辨識中。文獻[6]以一階船舶操縱響應模型為辨識對象，設計了一種基于標準最小二乘法的在線參數辨識試驗平臺，結合實船試驗數據驗證了在線辨識試驗平臺設計的正確性和合理性。文獻[7]通過引入P型學習率，提出了一種快速收斂迭代學習最小二乘算法，用于辨識船舶的一階非線性響應模型和二階非線性響應模型，文中仿真和水池實驗均表明，該算法相比標準最小二乘法在收斂速度和辨識精度都有所提高。文獻[8]針對船舶二階非線性響應模型，利用多新息最小二乘法對其進行參數辨識，辨識結果表明該算法辨識的參數更接近真實值，能更準確地描述船舶運動的動態特性。文獻[9]針對裝載吊艙式推進器的無人水面艇，利用最小二乘法辨識了其一階船舶操縱響應模型，對辨識模型進行回轉試驗的結果同實船回轉試驗結果的比對證明了辨識得到模型的正確性。文獻[10]針對船舶實際航行中存在海洋環境擾動和數據欠激勵的情況下，提出了一種滿秩分解最小二乘算法對船舶一階線性響應模型的參數進行辨識，辨識結果表明了該算法能抑參數發散，并能夠獲得更高的參數辨識精度。文獻[11]對最小二乘支持向量機算法中的偏置項引入了懲戒因子進行了改進，大大改善了辨識模型的精度；文獻[12]利用滾動時間窗內的數據實時辨識參數的策略，增強了對船舶操縱響應模型的參數跟蹤能力；文獻[13]針對最小二乘支持向量機引入自適應加權的策略進行改進，并結合實船試驗驗證了該方法的有效性。文獻[14]將其應用于船舶的航向控制，增強了船舶航向保持精度、降低了航向調節時間和超調幅度；文獻[15]從解析模型的角度出發，結合格雷碼和精英選擇的策略對參數進行非線性尋優，使得該算法在小樣本數據情況下就可進行高效準確的辨識，極大地提升了辨識效率。文獻[16]針對多新息擴展卡爾曼濾波算法引入了遺忘因子進行改進，有效較低了歷史累積數據對參數辨識的影響，文中結合船模數據的辨識結果表明，該改進算法相比標準擴展卡爾曼濾波辨識方法在船舶操縱響應模型參數辨識方面更加精確。

上述的文獻方法在模型參數辨識精度、收斂速度和參數實時跟蹤性能上均有貢獻，但都并未考慮受到有色噪聲干擾的實際量測數據會給參數辨識帶來有偏估計的影響。實際船舶航行過程中，由于海洋環境高頻擾動和傳感器測量噪聲等因素，會使得辨識的船舶操縱性指數存在精度上的誤差，這一缺陷給船舶的航向控制和預報都帶來了一定的影響，因此如何快速無偏地辨識船舶的操縱性指數是實現船舶操縱響應模型在線辨識中的關鍵問題。

考慮到實際過程中，噪聲對系統的影響是不可避免的。此時，要獲得系統模型參數的一致無偏估計, 遞推輔助變量最小二乘法(recursive instrumental variable least squares method, RIVLS)是一種有效的辨識方法。本文在RIVLS算法的基礎上，結合動態濾波的思想對其進行改進，提出了一種遞推修正輔助變量最小二乘法(recursive refined instrumental variable least squares method, RRIVLS)。該算法能夠對實時觀測得到的數據進行一次濾波預處理后，再結合輔助變量對模型參數進行辨識。辨識試驗以二階船舶操縱響應模型為辨識對象，在實船Z形試驗數據的基礎上，將遞推修正輔助變量最小二乘法、輔助變量最小二乘法、增廣最小二乘法(recursive extended least squares method, RELS)和標準最小二乘法(recursive least squares method, RLS)的辨識模型進行比較。

1 船舶操縱響應模型

船舶操縱響應模型是通過描述船舶航向角相對于舵角動態響應關系建立的模型。本文辨識的響應模型為二階線性響應模型，其方程為：

(1)

式中：r為船舶艏搖角速度；δr為外界環境干擾力矩造成的等效壓舵角；K為舵角增益系數；T1、T2和T3為時間常數。

考慮到艏搖角速度r不易通過傳感器直接獲取，因此對其利用航向角ψ前向差分得到：

(2)

式中h為采樣周期。

(3)

(4)

根據方程轉換法，令H=T1+T2，G=T1T2，同時結合式(2)～(4)，對式(1)進行相應的變換，可得到如下的受控自回歸模型(controlled auto-regressive, CAR)：

r(k)=a1r(k-1)+a2r(k-2)+

b1δ(k-1)+b2δ(k-2)+c

(5)

式中：a1、a2、b1、b2和c為待辨識參數，與K、T1、T2、T3和δr的關系為：

(6)

式(5)的模型可以描述船舶在風浪不太大的環境下的船舶操縱運動特性，對其加入噪聲干擾的影響，可得到如下的CAR模型：

r(k)=a1r(k-1)+a2r(k-2)+

b1δ(k-1)+b2δ(k-2)+c+e(k)

(7)

式中e(k)為船舶所受的有色噪聲，將該有色噪聲假設為自回歸過程(AR模型)，其模型的具體結構形式為：

(8)

式中：v(k)為均值為0的不相關隨機白噪聲；C(z-1)為單位后移算子z-1的多項式：

C(z-1)=1+c1z-1+c2z-2+…+cncz-nc

式中nc為噪聲向量的階數。

若當前采樣的時刻為k，式(7)可改寫為：

z(k)=φT(k)θ+e(k)

(9)

其中：

z(k)=r(k)

φ(k)=[r(k-1)r(k-2)δ(k-1)δ(k-2) 1]T

式(9)對應的系統矩陣形式為：

Zk=Φkθ+Ek

(10)

式中：

2 輔助變量最小二乘法的原理和改進

2.1 遞推輔助變量最小二乘算法的原理

最小二乘算法是一種根據最小二乘準則，在參數空間尋找一組最優模型參數估計值用于擬合系統模型的參數辨識方法，即令系統測量的輸出值與辨識模型輸出值之間的誤差平方和最小。式(10)的參數估計值的最小二乘解可寫為：

(11)

觀察式(11)最右側的第2項可知，只有當該項的值依概率收斂為0時，即：

(12)

繼續分析式(12)，可變形為：

(13)

(14)

(15)

圖1 輔助變量構建原理Fig.1 Diagram of instrumental variable construction

φ(k)=[r(k-1)r(k-2)δ(k-1)δ(k-2)1]T

(16)

φ*(k)=[r*(k-1)r*(k-2)δ(k-1)δ(k-2)1]T

(17)

K(k)=P(k-1)φ*(k)[1+φT(k)P(k-1)φ*(k)]-1

(18)

(19)

P(k)=[I-K(k)φT(k)]P(k-1)

(20)

(21)

(22)

2.2 遞推輔助變量最小二乘算法的改進

A(z-1)r(k)=B(z-1)δ(k)+c+e(k)

(23)

其中：

(24)

如對式(23)兩邊同乘一個預濾波器F(z-1)，可得：

A(z-1)rf(k)=B(z-1)δf(k)+cf+ef(k)

(25)

其中：

(26)

顯然，若正確選擇預濾波器F(z-1)可實現對式(8)中多項式C(z-1)的對消，從而改善有色噪聲e(k)的統計特性。預濾波器F(z-1)的數學模型表示為：

F(z-1)e(k)=v(k)

(27)

式中：

F(z-1)=1+f1z-1+f2z-2+…+fnfz-nf

(28)

該濾波器是未知的、有限階的線性濾波器。自回歸有色噪聲e(k)在預濾波器F(z-1)的作用下將會逐步白化，這在一定程度上能夠降低系統的噪信比，提高自適應濾波器的遞推精度。相比常規的輔助變量法，這種處理問題的方法稱作修正的輔助變量法[17]。綜上所述，從式(16)～(28)，可以將遞推修正輔助變量最小二乘法(RRIVLS)歸納成：

φ(k)=[r(k-1)r(k-2)δ(k-1)δ(k-2) 1]T

(29)

φf(k)=[rf(k-1)rf(k-2)δf(k-1)δf(k-2) 1]T

(30)

(31)

(32)

(33)

(34)

(35)

(36)

(37)

(38)

(39)

(40)

(41)

(42)

(43)

(44)

(45)

(46)

式(29)～(46)的遞推過程表明，修正的輔助變量法的基本思想是先對觀測得到的數據進行一次濾波預處理后，再結合輔助變量法對模型參數進行辨識。其中所用的濾波模型實際上是一種動態模型，在整個迭代過程中不斷靠偏差信息來調整濾波器參數，使其逐步逼近于一個較好的濾波模型，從而改善噪聲的統計特性。理論上如果濾波模型選擇得合適，再經過迭代后便可對數據進行較好的白化處理，但是實際過程中的噪聲過程復雜且充滿不確定性，這使得數據白化處理的可靠性下降。因此RRIVLS算法在數據預濾波處理的基礎上將其與原有的輔助變量法相結合，進一步提升參數辨識精度，實現對模型參數的無偏估計。需要注意的是，當濾波器階數設置為nf=0，則F(z-1)=1，此時RRIVLS算法將退化為RIVLS算法。

3 試驗驗證

為驗證RRIVLS對船舶操縱響應模型的參數辨識效果，在相對平穩的海況下，以本研究團隊自主建造的科研試驗船在海洋環境下進行10°/10°和20°/20°Z形試驗。科研試驗船的船型和主尺度參數分別如圖2和表1所示。

圖2 科研試驗船Fig.2 Research vessel

表1 科研試驗船的主尺度參數Table 1 Principal dimension parameters of research vessel

試驗的采樣序列長度為1 000，采樣周期為0.2 s。Z形試驗的數據曲線如圖3所示，其中的艏搖角速度通過式(2)差分得到。注意到圖3中Z形試驗數據中的艏搖角速度在某些時刻出現了抖動，這是由于試驗的航向角數據是在真實海洋環境下進行采集，所采集的航向角數據不可避免的會受到海洋環境的高頻擾動和傳感器測量噪聲的影響，因此艏搖角速度的差分過程會將該類噪聲進一步放大，從而致使艏搖角速度出現抖動，這也是式(7)中所提及的有色噪聲，受到有色噪聲干擾的實際量測數據會給參數辨識帶來有偏估計的影響。因此，本文所提出的RRIVLS算法針對受到有色噪聲干擾的艏搖角速度數據進行辨識建模，力圖獲得參數的無偏估計。

圖3 Z形試驗數據曲線Fig.3 Zigzag test data curve

3.1 船舶操縱響應模型參數辨識結果及分析

利用RRIVLS(nf=1)的式(29)～(46)結合上述的Z形試驗數據對式(7)中的模型參數進行辨識，并將其與RIVLS、RELS和RLS算法進行比對，得到的參數辨識曲線圖分別如圖4、圖5所示。各算法的初始設置情況如下所示：

圖4 基于10°/10°Z形試驗數據的參數辨識曲線Fig.4 Parameter identification curve based on 10°/10° zigzag test

圖5 基于20°/20°Z形試驗數據的參數辨識曲線Fig.5 Parameter identification curve based on 20°/20° zigzag test

由圖4和圖5可知，在10°/10°Z形試驗中，四種算法在遞推至400次左右時，所辨識的參數均已呈現收斂趨勢；在20°/20°Z形試驗中，RRIVLS和RIVLS辨識的參數在遞推至200次左右均呈現收斂趨勢，而RELS算法辨識的參數a1以及RLS算法辨識的參數b1和b2并未呈現收斂趨勢，這是因為RLS算法由于自身的固有缺陷，容易受有色噪聲的影響，致使對參數的辨識過程逐步偏離真值，而RELS算法由于需要考慮對噪聲進行建模，若噪聲模型的階次選擇不當，也會致使參數無法收斂。RRIVLS和RIVLS采用輔助變量的方式，一定程度上抑制了有色噪聲干擾對參數辨識的影響，保證了參數的無偏估計性，且RRIVLS算法能夠對觀測得到的數據進行一次濾波預處理，使得參數的收斂過程相比RIVLS算法更為平滑。

綜合來看，各算法的收斂結果均有所差異，因此需要后續試驗進行辨識結果的驗證比對。各算法在10°/10°Z形試驗數據和20°/20°Z形試驗數據上的最終參數辨識結果分別如表2和表3所示。

表2 式(7)的參數辨識結果(10°/10°Z形試驗)Table 2 Parameter identification results of formula (7) (10°/10° zigzag test)

表3 式(7)的參數辨識結果(20°/20°Z形試驗)Table 3 Parameter identification results of formula (7) (20°/20° zigzag test)

對表2和表3的參數辨識結果根據式(6)計算，對應的K、T1、T2、T3和δr的計算結果如表4和表5所示。

表4 式(1)的參數計算結果(10°/10°Z形試驗)Table 4 Calculation results of parameters in formula (1) (10°/10° zigzag test)

表5 式(1)的參數計算結果(20°/20°Z形試驗)Table 5 Calculation results of parameters in formula (1) (20°/20° zigzag test)

考慮到系統辨識方法在辨識參數較多時，可能存在“參數相消”效應，即所辨識的參數雖然偏離各自的真值，但同樣可以使誤差函數極小化，從而使得模型的輸入輸出特性與試驗數據擬合得很好。因此為了驗證表4和表5的參數辨識結果的準確性，這里選取野本標準法[18]計算的K、T參數為標準進行比對，野本標準法是一種通過Z形試驗數據曲線在圖解分析的基礎上，直接求取船舶操縱響應模型K、T參數的計算方法，盡管該方法計算的參數精度也存在一定的誤差，但在一定程度上也能為表4和表5的參數辨識結果提供一定的參數真值理論參考標準。為了將表4和表5的辨識結果與野本標準法進行K、T參數的比對，同時考慮到船舶在進行低頻運動時，可將式(1)的二階線性響應模型近似降階為一階線性響應模型[19]：

式中T=T1+T2-T3。

對表4和表5的參數結果按照上述進行換算，其換算結果和野本標準法的計算結果分別如表6和表7所示。

表6 K、T值的計算結果表(10°/10°Z形試驗)Table 6 Calculation result table of K and T values (10°/10° zigzag test)

表7 K、T值的計算結果表(20°/20°Z形試驗)Table 7 Calculation result table of K and T values (20°/20° zigzag test)

由表6和表7的計算結果可知，野本標準法分別在10°/10°Z形試驗和20°/20°Z形試驗所計算的K、T和δr參數值具有一定的誤差，這是因為計算所采用的Z形試驗數據是在真實海洋環境條件下采集的，而海洋環境的實時變化都會對船舶操縱性參數產生一定的影響。因此，野本標準法通過不同的Z形試驗數據所計算得到的K、T和δr參數值具有一定的誤差是正常的。同時注意到表6和表7中RRIVLS、RIVLS和RELS所辨識參數換算的K、T值與野本標準法計算的K、T值都較為接近，只有RLS算法所辨識參數換算的K、T值與野本標準法的誤差較為顯著；而對于參數δr，RRIVLS、RIVLS和RELS辨識的參數δr值與野本標準法計算的δr值也都較為接近，而RLS算法所辨識的δr值與野本標準法計算的δr值在10°/10°Z形試驗中的誤差較為顯著。綜合來看，RRIVLS、RIVLS和RELS在不同試驗中辨識的參數換算結果K、T和δr參數值具有一定的合理性，但并不能證明T1、T2和T3的合理性。因此需要在后續開展模型驗證試驗，從而以試驗的角度證明T1、T2和T3的辨識準確性。

3.2 模型驗證

為了進一步對比表4和表5中各算法辨識參數結果的精度,從擬合和泛化2個角度出發對辨識模型的輸出特性進行驗證：

1)分別將表4和表5中各參數辨識結果代入式(1)中，再分別以10°/10°Z形試驗和20°/20°Z形試驗舵角數據為輸入，對式(1)的微分方程采用四階龍格庫塔法進行求解，從而可得到對應的模型輸出曲線。該角度可評價解算得到的航向角與實際測量航向角的擬合精度。各Z形擬合試驗的航向角比對曲線和擬合誤差曲線分別如圖6和圖7所示。

圖6 10°/10°Z形擬合試驗的航向角比對曲線和擬合誤差曲線Fig.6 Heading angle comparison curve and fitting error curve of 10°/10° zigzag fitting test

圖7 20°/20°Z形擬合試驗的航向角比對曲線和擬合誤差曲線Fig.7 Heading angle comparison curve and fitting error curve of 20°/20° zigzag fitting test

2)利用10°/10°Z形試驗得到的參數辨識結果即表4的參數辨識結果代入式(1)中，對其進行100 s的20°/20°Z形仿真試驗，仿真過程利用四階龍格庫塔法進行求解，從而可得到對應的20°/20°Z形航向角泛化曲線。該角度可評價表4中各算法辨識模型的泛化精度。20°/20°Z形泛化試驗的航向角比對曲線和泛化誤差曲線如圖8所示。

圖8 20°/20°Z形泛化試驗的航向角比對曲線和泛化誤差曲線Fig.8 Heading angle comparison curve and generalization error curve of 20°/20° zigzag generalization test

上述所提及誤差計算為：

由圖6和圖7可知，RRIVLS、RIVLS和RELS所辨識的模型相比RLS辨識的模型擁有更高的擬合精度，這一點可以從Z形試驗的航向角比對曲線和航向角擬合誤差曲線直觀看出。而由圖8的航向角泛化誤差曲線可以得出，各辨識模型的誤差曲線隨著時間推移在不斷的增大，這是因為所比對的20°/20°Z形試驗數據是在實際海洋環境中采集的，其中的操縱性參數必然會受到海洋環境的影響而產生變化，因此利用表4的參數進行仿真試驗所產生的誤差放大現象是正常的。同時，注意到在仿真曲線的前60 s，利用表4中RRIVLS、RIVLS和RELS算法辨識的參數結果所模擬的試驗數據曲線仍然能夠反映出實際中20°/20°Z形試驗的航向角變化特性，且RRIVLS相比其余3種辨識模型的泛化精度更高。這也側面說明了RRIVLS算法辨識得到參數的準確性和合理性。而RLS算法受到有色噪聲的影響，導致辨識的參數偏離真值，因此在圖6～8中的誤差表現顯著。

為了更具體地評價辨識模型的擬合精度和泛化精度，對圖6～8中的航向角誤差曲線利用均方根誤差(root mean squared error, RMSE)和最大絕對誤差(maximum absolute error, MAE)進行評價：

MAE=max{|ei|}i=1,2，…,N

式中N為樣本數。各Z形試驗的RMSE值和MAE值分別如表8～10所示。

表8 10°/10°Z形擬合試驗的RMSE值和MAE值評價表Table 8 RMSE value and MAE value evaluation table of 10°/10° zigzag fitting test

表9 20°/20°Z形擬合試驗的RMSE值和MAE值評價表Table 9 RMSE value and MAE value evaluation table of 20°/20° zigzag fitting test

表10 20°/20°Z形泛化試驗的RMSE值和MAE值評價表Table 10 RMSE value and MAE value evaluation table of 20°/20° zigzag generalization test

表8～10反映了圖6、圖7和圖8中各辨識模型輸出數據與實際航向角數據在不同誤差指標下的評價結果。評價結果表明，RRIVLS算法所辨識得到的模型在擬合角度和泛化角度的RMSE指標和MAE指標相比其余算法更小，符合圖6、圖7和圖8中的局部展示結果。且RRIVLS算法相比RIVLS算法，在RMSE指標和MAE指標上均有提升，從而驗證了針對RIVLS算法改進的有效性。

3.3 算法實時性分析

本文辨識船舶操縱響應模型的過程中，使用了實船10°/10°Z形試驗和20°/20°Z形試驗各1 000組數據進行辨識試驗，在Windows10系統下的Octave平臺中記錄了RRIVLS、RIVLS、RELS和 RLS 算法的計算耗時，具體如表11所示。

表11 各算法的計算耗時Table 11 The running time of each algorithm

從算法耗時統計結果可知，RRIVLS算法的平均單次耗時為0.02 ms，耗時為另外3種算法的2～3倍。這是由于RRIVLS算法額外考慮了對觀測數據的動態濾波，在辨識過程需要額外辨識濾波器的參數，導致了每次辨識循環需要進行2次參數辨識和一次數據濾波，這個過程會額外增加算法的計算耗時，犧牲了RRIVLS算法的實時性，用以提升算法的參數辨識精度，綜合來看RRIVLS算法的實時性在可接受范圍以內。

4 結論

1)遞推修正輔助變量最小二乘法所辨識的模型相比其余三者的辨識模型與船舶實際航向的擬合誤差和泛化誤差更小，其對應的均方根擬合誤差和均方根泛化誤差分別可達到1°和2°以下，最大絕對誤差分別在3°和4°以下。這說明了修正的遞推輔助變量最小二乘法在實船航向過程中能夠有效抑制有色噪聲的影響，提升了模型辨識的精度，從而為實現精確穩定的航向自適應控制或構建航向在線預報模型提供了基礎。

2) 本文的辨識工作僅考慮了二階線性船舶操縱響應模型，對于船舶的線性運動和中小幅運動比較適應，若船舶處于大舵角或強非線性運動過程中，利用該模型所描述的動態響應過程與實際過程可能會出現較大的偏差。因此，未來可以針對二階非線性船舶操縱響應模型展開辨識研究。