一種共軸雙旋翼飛行器懸停控制聯合仿真

陳漢， 李科偉， 鄧宏彬， 危怡然， 趙瑾

(1.北京理工大學 機電學院, 北京 100081; 2.中國兵器工業集團有限公司 第203研究所, 陜西 西安 710065)

0 引言

無人機近年來受到越來越多科學家及工程師的關注，尤其是微小型無人機。微型旋翼機比全尺寸直升機更靈活，其尺寸緊湊的特點以及懸停、轉彎和向各方向移動的能力使小型旋翼航空器非常適合動態環境下的操作。與固定翼飛機相比，微型旋翼飛機具有顯著優勢，特別是在需要保留飛機靜止(懸停)時或在嚴格受限的環境中操縱。在戰場環境中，微型旋翼無人機可以執行多種任務，例如偵查環境、電子對抗、通訊中繼甚至直接攻擊目標等[1]。

由于可操作性、多功能性、穩定性和易于控制等優勢，再考慮到實際的工程項目需求，飛行器能夠攜帶質量較大的任務載荷，本文選擇微型共軸旋翼機而不是單旋翼的配置方案。與目前最流行的四旋翼飛行器對比，共軸雙旋翼飛行器的結構更加緊湊，盡管目前存在一些可折疊四旋翼飛行器，但四旋翼飛行器無法滿足折疊后的尺寸要求，而本文所涉及的共軸雙旋翼飛行器其翼片折疊后特征尺寸顯著減小，更加便攜，適合放置在狹長的火箭彈艙室內運輸并投放，且折疊旋翼比折疊機身結構設計更簡單可靠，因此本文采用共軸雙旋翼的配置方案[2]。

然而，四旋翼飛行器的相關控制算法相比之下更加簡單，研究已經非常成熟，而共軸雙旋翼飛行器的相關研究則較少，尤其是國內關于微小型共軸雙旋翼飛行器控制算法的文獻幾乎沒有。目前，國外相關的研究工作有：美國Drexel大學Husnic[1]對一種微型共軸雙旋翼直升機設計了控制器，并利用MATLAB/Flightgear進行了聯合仿真；Maryland大學Lee[3]對雙旋翼涵道飛行器的懸停和前進飛行性能進行了研究；Duke大學Giovanetti[4]對共軸雙旋翼的最小功率控制和轉子優化設計進行了研究；Purdue大學Arama[5]對一種小型共軸雙旋翼飛機設計了混合控制策略，在仿真中得到了較好的結果。另外Gecgel[6]、孫茂等[7]、Wicha[8]在共軸雙旋翼的空氣動力特性和優化功率損耗方面均做出了獨立研究。國內的相關研究工作有：南京航空航天大學寧承威[9]對共軸傾轉旋翼機的飛行動力學建模與控制進行了研究；北京理工大學黃祥斌[10]對一種采用空氣舵控制的球形共軸旋翼飛行器進行了總體設計，包括控制算法和飛行實驗；清華大學袁夏明等[11-12]對共軸式無人直升機建模與魯棒跟蹤控制進行了深入細致的研究，完善了相關理論；石征錦等[13]對Fuzzy-PID算法在共軸雙旋翼姿態控制的應用進行了研究；匡銀虎等[14]研究了模糊自適應算法在多旋翼飛行器控制中的應用；李悅[15]對共軸雙旋翼涵道飛行器的氣動特性和仿真方法進行了獨立研究。

從控制器設計和仿真技術來看，國內外相關研究普遍基于對共軸雙旋翼的簡化或近似的動力學模型對控制算法進行設計，并在MATLAB/Simulink仿真平臺上對控制算法進行性能測試。這種做法存在兩方面問題：首先，只能依靠經驗公式和近似理論求解上旋翼和下旋翼干擾的問題，在樣機測試之前無法準確判斷不同系統狀態下的升力大小和方向，模型仍存在系統誤差。如果控制器對模型的依賴較強，即使在仿真中可能取得較好的成果，在樣機實驗中也會因為未建模的誤差、裝配和加工精度、舵機精度等問題造成姿態發散。其次，MATLAB中難以對飛行器的復雜操縱機構進行準確建模，目前仿真方法使用的都是簡化數學模型，將其整體考慮為一個剛體加上外部控制力的方式，忽略了機體上各部件的相對運動。顯然這種做法是不準確的，沒有建立起舵機動作量和控制算法輸出控制量之間的聯系，控制器輸出量和實際情況差距很大。

基于以上現狀，本文主要在以下兩方面實現了創新：第一，采用適應性和魯棒性更好的滑模控制算法，應用于控制執行機構復雜的共軸反槳式無人機上，使姿態控制更加迅速穩定；第二，將Adams/MATLAB聯合仿真技術應用于飛行器控制仿真，解決具有復雜操縱機構的無人機難以進行數學建模和數值仿真的問題，利用聯合仿真技術，較好地對實際飛行器動力學系統包括舵機及其傳動機構進行建模，可信度更高[16]。從而使得本文設計的懸停控制算法具備可靠性和有效性，并具有可視化的優點，通過動畫保存可以直觀地找出設計中的錯誤，成果展示的效果也更好。

1 飛行器操縱原理

本文所涉及的共軸雙旋翼飛行器由任務載荷、電池倉、機身、下槳盤、上槳盤等5大部分組成，所有結構的三維建模均在Solidworks 2018軟件中完成，如圖1和圖2所示。

圖1 翼片折疊后Fig.1 Folded wings

圖2 共軸雙旋翼飛行器結構組成Fig.2 Structure of coaxial dual-rotor aircraft

由圖2可知，飛行器電池倉和任務載荷都可以看作是飛行器載荷，它們與機體固連；上槳盤和下槳盤的槳葉均可折疊，飛行時通過離心力展開并達到水平位置，槳葉與槳盤鉸接副的活動范圍為90°. 下面介紹機身上舵機控制航向的方式，操縱機構如圖3所示。

圖3 操縱機構Fig.3 Control mechanism

圖3中所有的連桿兩端均為球鉸約束，由步進電機帶動絞盤和連桿，與上方的操縱環連接，操縱環內是1個球面軸承，使操縱環中心點與飛行器中心軸線重合，但是可以傾斜及沿軸線上下運動。根據幾何知識，空間中的3個點可以確定1個平面，因此操縱環的位置和姿態可以通過3個舵機完全確定。當操縱環平面法向量不與機身中心軸線重合時，上槳葉迎角是周期性變化的，槳葉在轉動過程中在連桿帶動下其迎角與下方操縱環的姿態對應，在操縱環較高處槳葉的迎角也會變大，從而使產生的升力變大，在操縱環較低的位置則相反。最后的效果就是上槳盤所產生的所有槳葉升力合力方向與操縱環平面的法向量方向完全一致。

當飛行器向前飛行時，槳葉開始揮舞運動，這種運動較為復雜，會給控制器設計帶來挑戰。本文所涉及的飛行器在結構設計上限制了揮舞運動的幅度、減輕了其影響，并在動力學建模中考慮了揮舞運動帶來的周期性變距問題，采用近似方法對其進行建模，同時采用滑模控制解決建模不精確問題[17]。

2 飛行器動力學模型

2.1 無上旋翼和下旋翼干擾及揮舞運動的模型

為建模方便，現引入2個坐標系，即大地坐標系OXYZ和機體坐標系oxyz. 其中OXYZ坐標系可視為慣性坐標系，oxyz坐標系為機體坐標系，且oxyz的原點位于飛行器質心位置。飛行器坐標系以及簡化模型如圖4所示。

圖4 飛行器坐標系以及簡化模型Fig.4 Aircraft coordinate system and simplified model

圖4中，φ、θ、ψ表示歐拉角，即分別為滾轉角、俯仰角和偏航角，F1、F2分別為下旋翼和上旋翼的升力，R為槳盤半徑，L為上旋翼到下旋翼安裝中心的距離。

為了便于分析，在確保研究結果有效的前提下，本文提出以下3點假設：

1)將整個共軸雙旋翼飛行器視作1個剛體；

2)將共軸雙旋翼飛行器質心與中心對稱軸視作重合；

3)由于旋翼體積和質量很小，可忽略旋翼陀螺效應對系統動力學特性的影響[18]。

現定義V=(p,q,r,u,v,w)為機體速度向量，即共軸雙旋翼在機體坐標系下的速度向量，其中p、q、r、u、v、w分別為機體相對于x軸、y軸、z軸的角速度以及質心相對于x軸、y軸、z軸的速度。

基于歐拉- 龐卡萊方程建立共軸雙旋翼無人機動力學模型，在機體坐標系下其方程的一般形式為

(1)

式中：Q為機體質心相對于大地坐標系下的廣義坐標向量；U為控制力；M(Q)為慣性矩陣；C(Q,V)為根據歐拉動力學方程得出的變換矩陣；F(V,Q,U)為機體坐標系下包含空氣動力、重力以及控制輸入的總和。

(1)式中部分矩陣的具體表達式如下：

(2)

(3)

式中：m為飛行器質量；Ix、Iy和Iz分別為飛行器相對于x軸、y軸、z軸的轉動慣量。

(4)

式中：F2x和F2y分別為F2在機體坐標系內x軸和y軸上的投影。將z1作為調整上旋翼系統總距離和升力大小的基準，讓其隨著控制器輸出的F2期望值變化，則可以求解出z2和z3.

圖5 上旋翼拉力分解Fig.5 Decomposition of upper rotor pulling force

飛行器質心處的受力情況由(5)式給出：

(5)

式中：

(6)

U1為垂直高度控制量，U2為滾轉輸入控制量，U3為俯仰控制輸入量，U4為偏航控制量，η為旋轉力矩系數(與旋翼受到的阻力和電機軸摩擦系數有關)；V(Q)為飛行器機體坐標系向大地坐標系轉換的矩陣，為簡化書寫，記s()=sin (),c()=cos()，

(7)

旋翼拉力可由(8)式計算得到：

(8)

式中：S為槳葉迎風面積；ρ為空氣密度；Ct為旋翼拉力系數；ω為旋翼旋轉角速度。根據(5)式和(6)式聯立后將等式兩側同時從機體坐標系轉換到大地坐標系，可以得到大地坐標系下的廣義坐標向量對時間的導數：

(9)

共軸雙旋翼飛行器在無風或風速很小狀態下低速飛行時，空氣阻力對系統的影響較小，可以忽略不計。同時，假設共軸雙旋翼在飛行過程中的滾轉角和俯仰角很小，且其變化率也充分小，則系統的數學模型可以簡化為

(10)

由(10)式可知，飛行器質心位置的3個坐標x、y、z這3個通道之間的耦合非常嚴重，在考慮旋翼間氣動干擾對入流動力學的影響后，其力與力矩的動態變化過程更復雜，在飛行器本體動力學的非線性之外進一步疊加了控制力產生過程的嚴重非線性，同時也加重了耦合問題。

2.2 上旋翼和下旋翼干擾及揮舞運動的模型

為了反映上旋翼和下旋翼之間的氣動干擾，根據葉素理論和Pitt-Peters動態入流模型建立上旋翼和下旋翼入流干擾及揮舞運動模型。該方法相比渦流理論或計算流體力學(CFD)方法的形式簡潔、計算方便，與實驗數據也吻合較好。需要指出的是，該模型建立在以下假設基礎上：槳葉為剛性，忽略槳葉根部截斷效應、槳尖損失與揮舞鉸外伸量，且不考慮非定常效應，認為入流速度在槳盤平面內均勻分布[12]。Pitt-Peters動態入流模型的核心思想是將誘導速度的動態變化和空氣阻力變化通過(11)式聯系起來：

(11)

式中：λi=[λ0,λs,λc]T，λ0、λs和λc分別表示誘導入流比的時均、1階橫向分量和1階縱向分量；C=[Ct,Cr,Cp]T，Cr和Cp分別為旋翼滾轉力矩和俯仰力矩系數；Mi、Vm和L分別為入流動力學的慣性矩陣、質量流量參數矩陣和靜態傳遞矩陣。誘導速度的相互影響用(12)式表示為

λi=λ0,i+Kjiλ0,j+λfs，

(12)

式中：Kji表示旋翼間誘導速度干擾系數，與兩組旋翼的直徑和間距等設計參數有關；i、j∈{u,l},i≠j，下標l表示下旋翼，u表示上旋翼；λi為入流比；λ0,i、λ0，j為靜止狀態下無窮遠處到某一葉片的入流比；λfs為飛行器運動帶來的入流比。根據上述入流干擾模型，可以將旋翼系統的拉力和扭矩表示如下：

(13)

(14)

式中：T和Q分別為整個旋翼系統產生的拉力和扭矩；δ為葉素槳距角；Cd為槳葉阻力系數；v0,i為無窮遠處到某一葉片的誘導速度；σ=Nbc/(πR)表示槳盤實度，c為槳葉弦長，Nb為槳盤上的槳葉根數；a0為翼片升力線斜率；A為槳盤面積。

另外，考慮到槳葉的剛度，揮舞運動對另外俯仰和滾轉通道的正交力矩為

(15)

式中：k9=NbSβγIβΩ2/8，Sβ為葉片剛度系數，γ為槳葉洛克數，Iβ為揮舞轉動慣量；αa為縱向揮舞角；αb為橫向揮舞角。

綜上所述，考慮揮舞運動和旋翼氣動干擾后，共軸雙旋翼飛行器的動力學模型修改為(其他項不變)

(16)

3 控制器設計

3.1 總體控制結構

圖6所示為共軸雙旋翼飛行器系統的控制結構示意圖，整個控制系統可分為姿態控制和位置控制。

圖6 控制結構圖Fig.6 Control structure

由于共軸雙旋翼的三軸姿態控制之間不存在耦合關系，且姿態控制響應最快，本文優先對姿態進行控制，保證其精確控制。考慮到共軸雙旋翼的欠驅動問題，即4個輸入控制6個自由度，直接通過PID算法建立姿態角與水平位置之間的數學關系，通過水平位置的信息解算出所需姿態角的信息，再根據姿態角和水平控制量計算高度控制量，進而通過控制姿態角和位置高度實現軌跡跟蹤的目的[19]。滾轉角與俯仰角的期望值分別表示如下：

(17)

式中：kp、ki、kd分別為比例、積分、微分的調節系數；φd、θd、ψd分別為期望的俯仰角、滾轉角、偏航角；Xd、Yd為期望值；X、Y可以由機載傳感器通過測量得到。

3.2 滑模PID控制算法

因為通常采用線性滑模面的滑模控制器會有明顯的抖振，且無法通過調整參數來消除抖振，所以本文在滑動模態控制算法的基礎上進行改進，將滑模面函數設計成誤差項的比例、微分、積分3項相加的形式，這種算法即為滑模PID控制算法[19]。這種做法的目的在于改善滑模面的性能，相當于在滑模面附近建立一個“緩沖帶”，以有效地抑制抖振。滑模控制器的一般設計步驟為：1)選取適當的PID滑模面函數s；2)設計合適的控制律，使該系統能夠到達且保持在期望的滑模面s=0上。

以共軸雙旋翼俯仰通道為例，定義跟蹤誤差為

eφ=φd-φ，

(18)

然后選取滑模面函數為

(19)

式中：λ1、λ2為控制器參數。

再根據趨近律推導控制律方程，選取指數趨近律，即

(20)

式中：ε、k為控制器參數，ε=0.2，k=1.

選取Lyapunov函數為

(21)

則有

(22)

由此可知，四旋翼系統漸進穩定，跟蹤誤差最終收斂到0，且滿足滑動模態的可達條件。

因此，由

(23)

可得

(24)

式中：c1、c2為控制器參數。

同理可得滾轉角、偏航角和高度的控制律分別為

(25)

式中：d1、d2、h1、h2為控制器參數。

4 Adams/MATLAB聯合仿真

4.1 仿真模型建立

首先將Solidworks2018軟件中建立的裝配體進行必要的簡化，按照相對運動關系，將原本的202個零部件中在飛行時不存在相對運動的一組生成單個零件，則簡化后整個飛行器由18個零件組成。重新將18個零件進行裝配后，再依次將每個零件另存為.x_t格式，并導入Adams View中添加約束關系、力和力矩、輸入及輸出變量，最終輸出動力學模型文件。由于Adams不具備空氣動力學仿真的能力，公式中的誘導速度、槳葉阻力系數和揮舞角均取自該模型的空氣動力學仿真和實驗，然后將兩副旋翼的拉力和扭矩作用點放置在槳盤中心，力的數值參照(11)式以函數的形式定義。

Adams中的動力學模型如圖7所示。輸出的文件是Simulink可讀取的.slx格式，文件名為adams_sys，打開后如圖8所示。動力學模型以adams_sub模塊的形式與Simulink中搭建的控制模型連接，其輸出就是飛行器當前的3個姿態角和質心空間坐標。

圖7 Adams中的動力學模型Fig.7 Dynamics model in Adams

圖8 Simulink中的動力學子模塊Fig.8 Dynamics submodule in Simulink

在Simulink中將動力學模塊的輸入量和輸出量通過控制算法建立起關系，在MATLAB/Simulink軟件中將模塊組裝成控制器，并添加適當的自定義函數。由于控制律中存在關于飛行器質量和轉動慣量等信息，需要在程序運行之前對這些信息進行初始化，本文采用運行m文件的方法，將所有初始化的常量寫入同一文件中，運行后上傳到Workspace中，運行Simulink的模型后便可以被調用，從而方便了統一修改和管理，比其他方法更加簡便[20]。

整個仿真模型搭建完成后如圖9所示，橙色模塊即為Adams中的動力學模型。

圖9 Simulink中的仿真模型Fig.9 Simulation model in Simulink

由于直接對飛行器上的3個舵機進行控制，這里對飛行器動力學模型輸入的量為3個舵機偏轉角度τ1、τ2、τ3以及上、下2個電動馬達的轉速ω1和ω2. 3個舵機選用相同型號規格的步進電機，2個電動馬達也完全相同，暫不考慮它們內部的響應問題。上文中已經建立了4個控制量U1、U2、U3、U4的關系式，下面對其進行轉化，設計一定的規則，將U1、U2、U3、U4轉化為τ1、τ2、τ3、ω1和ω2共5個輸入量。

(26)

式中：Δz1、Δz2、Δz3分別為3個舵機對上旋翼升力面相應支撐點造成的z軸方向位移范圍;Δτ為相應的步進電機轉動范圍，聯立(4)式即可求出5個輸入量。

4.2 懸停仿真結果

系統仿真的初始值為

X0=[0 m，0 m，0 m，0 rad，0 rad，0.1 rad，0 m/s，
0 m/s，0 m/s，0 rad/s，0 rad/s，0 rad/s]，

式中：從左至右依次為系統位置x、y、z以及姿態角φ、θ、ψ、速度和角速度。

控制器中的參數按照固定步長在一定范圍內循環搜索測試，逐步縮小范圍和步長，最終確定了某種意義上的最優參數組合，如表1所示。

表1 控制器參數表

飛行器結構尺寸參數和氣動參數如表2所示。

表2 飛行器動力學參數

誘導速度場實際上并不是空間均勻分布的，并且與旋翼轉速有關，這里只取1個平均值。當設置目標位置坐標與初始位置相同時，仿真得到的三維運動軌跡以及位置和姿態隨時間變化曲線如圖10和圖11所示。

圖10 原地懸停軌跡Fig.10 In-place hover trajectory

圖11 懸停仿真結果Fig.11 Hover simulated results

三維軌跡圖中紅色圓圈代表懸停的目標位置(0 m，0 m，0 m)，從位移隨時間變化曲線中可以看出：飛行器的懸停位置基本穩定，并隨時間變化逐漸趨于穩定，后半段的較穩定狀態時飛行器的位置三軸坐標均在(-0.05 m，0.05 m)范圍內；姿態角一直存在小幅度震蕩，這是因為用PID算法計算姿態角動態誤差較大以及對系統建模仍然存在一定誤差。從位置誤差曲線可以看出，在后半段較穩定狀態時飛行器的懸停位置誤差在(-0.05 m，0.05 m)范圍內，姿態角與姿態角期望值之間的誤差幾乎為0°，表明系統姿態跟蹤性能較好。

程序運行完成后，會生成.res文件記錄仿真結果。在Adams中導入后可以觀察此次運行的動畫，通過動畫直觀感受飛行器的運動狀態，找出需要改進的位置。圖12采用圖像疊加的方法將動畫中從某一固定觀察角度且相同間隔時刻的畫面以不同透明度顯示到同一張圖片中，更加直觀地展示了飛行器運動過程。

圖12 懸停動畫截圖Fig.12 Hover animation screenshots

當設置目標位置坐標與初始位置不同時，設定目標位置坐標(2 m，2 m，2 m)，仿真得到的三維運動軌跡以及位置、姿態隨時間變化曲線和動畫截圖如圖13、圖14和圖15所示。

圖13 移動懸停軌跡Fig.13 Move-hover trajectory

圖14 移動懸停仿真結果Fig.14 Move-hover simulated results

圖15 移動懸停動畫截圖Fig.15 Move-hover animation screenshots

由圖13～圖15可以看到，在經過初始短暫的超調和震蕩后，依舊可以穩定在目標位置附近，說明本文所設計的控制器性能可以滿足要求。

4.3 抗干擾仿真結果

在Adams View中，在飛行器質心處添加3個正交力和3個正交力矩，數值大小利用隨機函數生成，用來模擬施加干擾力和力矩，測試控制器的抗干擾性能。干擾力和力矩的最大值分別為2 N和1 N·m，干擾施加的頻率設置為10 Hz，干擾作用的時間范圍是仿真開始后的20～25 s. 加入干擾后的仿真結果如圖16所示。

圖16 加入干擾后懸停仿真結果Fig.16 Hover simulated results under interference

由圖16可以看到，在持續作用的隨機干擾結束后，飛行器姿態可以在5 s后恢復到穩定狀態，質心位置可以在8 s后恢復到懸停的目標位置。在干擾作用5 s內，姿態角仍然可以穩定在較小的范圍內，質心位置有發散趨勢，這是因為干擾力變化頻率太快，超過了飛行器機動調整能力，使得飛行器姿態誤差一直保持在較大的范圍內，這種誤差經過積分放大后就導致了質心位置的偏移。如果干擾作用時間較短，類似于脈沖形式，則飛行器位置質心受到的影響就大大減小。加入脈沖干擾后懸停軌跡坐標如圖17所示。

圖17 加入脈沖干擾后懸停軌跡坐標Fig.17 Hover simulated results under pulse interference

5 結論

本文首先介紹了實驗室在研的一種共軸雙旋翼飛行器結構和工作原理；然后結合Pitt-Peters動態入流模型，建立了共軸雙旋翼飛行器的動力學模型，并基于滑模PID算法設計了對應的控制器；最后利用Adams/MATLAB聯合仿真技術進行了飛行器飛行控制仿真，分別對懸停，移動懸停，懸停中持續干擾和脈沖干擾4種情況的仿真結果進行了分析。所得主要結論如下：

1)本文所設計的控制器直接對飛行器舵機和電動馬達進行控制，使MATLAB中控制模型與實際的飛行器飛控芯片內代碼非常接近，從而使得仿真結果更具有可信性。

2)從仿真結果中可以看出，飛行器具備了較好的懸停穩定性、抗干擾性能和一定的機動能力。飛行器在控制器作用下的位置和姿態均得到了有效控制，懸停性能完全滿足要求。

3)在仿真結果的穩定階段(t>20 s)，姿態角仍然存在震蕩，意味著舵機和電機工作狀態也在高頻率的震蕩，會對舵機和電機的壽命造成一定的損害。主要原因為通過PID算法直接建立期望姿態和水平位置的映射關系時并沒有考慮到飛行器動力學屬性，下面的工作中重點考慮改進這一部分算法，使飛行器姿態得到更好的收斂。