高速近距無人機協(xié)同編隊控制器設計研究

張佳龍, 閆建國, 張普, 齊曉巧, 呂茂隆

1.西北工業(yè)大學 自動化學院, 陜西 西安 710129; 2.西安電子科技大學 機電工程學院, 陜西 西安 710071; 3.空軍工程大學 裝備管理與安全工程學院, 陜西 西安 710051

自從無人機(unmanned aerial vehicle，UAV)在1991年海灣戰(zhàn)爭中得到成功運用以來，經(jīng)過幾十年的發(fā)展，無人機技術(shù)已相對成熟，并在各個領(lǐng)域中發(fā)揮了其獨特的作用[1]。單架無人機在作戰(zhàn)范圍、殺傷半徑、摧毀能力以及攻擊精度等方面受到了限制，影響整個作戰(zhàn)任務的成功率；另外，一旦單架無人機中途出現(xiàn)故障，必須立即中斷任務返回，但在戰(zhàn)爭中有可能貽誤無人機而破壞整個作戰(zhàn)計劃。因此，多無人機協(xié)同編隊能最大限度發(fā)揮單機優(yōu)勢，實現(xiàn)多無人機協(xié)同編隊飛行(multi-UAV coordinated formation flight)控制、決策和管理，從而提高無人機完成任務的效率，拓寬無人機使用范圍，達到安全、高可靠性地執(zhí)行各種任務的目的[2]。

多無人機協(xié)同編隊在理論研究和工程應用上具有重要的意義，并已取得了不少研究成果。目前，無人機協(xié)同編隊問題的研究已引起國內(nèi)外學者的廣泛關(guān)注[3-6]。無人機可看作是一種空中機器人，不考慮其實際的飛行動力學特性，可將其抽象為一種具有自主能力的智能體[7]。此時無人機協(xié)同穩(wěn)定性控制就轉(zhuǎn)化為多智能體系統(tǒng)的穩(wěn)定控制。目前，有關(guān)多智能體系統(tǒng)的研究主要集中在一致性問題[8-12]、群集問題[13-14]和編隊問題[15-16]。對于多智能體系統(tǒng)的編隊控制來說，其穩(wěn)定性是指多智能體系統(tǒng)的運動同步[17]。文獻[18]研究了具有非線性和隨機切換拓撲智能體網(wǎng)絡的分布式同步問題，多智能單元的多無人機協(xié)同作戰(zhàn)集中式結(jié)構(gòu)。協(xié)同編隊控制中，隊形保持是一個主要問題。為實現(xiàn)無人機群以預期的姿態(tài)，恒定的間距飛行，需要控制每架無人機速度，航向角等飛行信息，使其達到一致性目的。2012年，祁圣君等人提出一種無人機近距編隊飛行控制方法，并以雙機近距離編隊為基礎(chǔ)對控制率設計和算法實現(xiàn)研究[19]。2015年，美國"全球鷹"實現(xiàn)了兩無人機近距編隊飛行，“空中軟式加油”對接，是一個典型的近距協(xié)同編隊飛行過程。大中型固定翼無人機協(xié)同編隊高速飛行過程中，集群中單機出現(xiàn)故障，需要重新近距或超近距編隊飛行，已經(jīng)成為研究的熱點。故本文提出解決在多無人機高速編隊飛行中近距協(xié)同控制的問題。

本文針對無人機編隊在執(zhí)行作戰(zhàn)任務時，機群中單機出現(xiàn)故障，需退出編隊飛行(如無人機目標被擊中、單機故障不適宜繼續(xù)在編隊中，隊形被破壞需重新組隊編隊)，即系統(tǒng)協(xié)同穩(wěn)定性被打破，不能完成作戰(zhàn)任務，提出了一種無人機編隊協(xié)同飛行控制器的應急故障容錯設計方法，提高飛行協(xié)同編隊的穩(wěn)定性，能夠及時調(diào)整編隊中無人機的重新編排和信息分配，達到多無人機協(xié)同編隊飛行的目的。

1 任務描述

假設在執(zhí)行一次作戰(zhàn)任務中，由n(n>1)架無人組成的編隊分別從不同的機場起飛，前往指定區(qū)域執(zhí)行作戰(zhàn)任務。在任務執(zhí)行之前，需要考慮敵方的綜合實力，如防空系統(tǒng)，對地攻擊系統(tǒng)，攔截戰(zhàn)術(shù)彈道和巡航導彈的能力，禁飛區(qū)以及無人機自身約束性能等，這些是我方作為出動兵力進行精確打擊的參考。為了提升打擊行動的成功率，我方無人機編隊勢必要求相互配合，相互協(xié)同，而且集群間通過數(shù)據(jù)鏈共享信息，進行統(tǒng)一決策，協(xié)同分工，保持“長機-僚機”之間的距離，盡可能按比例縮減無人機之間的相對位置和速度。這對無人機在復雜戰(zhàn)場環(huán)境中安全高效、快速準確執(zhí)行作戰(zhàn)任務是至關(guān)重要的。

對于垂直打擊任務中，協(xié)同編隊控制是必然且極其重要的一個環(huán)節(jié)。進行無人機協(xié)同編隊控制器設計，首先攜帶雷達的預警機(即“長機”)對無人機群進行信息感知，并對信息進行融合，其次在執(zhí)行任務過程中，其中一架無人機故障，無人機編隊重新編排，預警機將已感知信息進行融合處理，(重新)分配，快速決策，使每架無人機生成期望的航跡進行重新規(guī)劃，然后利用所研究的協(xié)同控制理論和隊形設計方法實現(xiàn)協(xié)同飛行。本文設計過程中，認為“長機-僚機”之間的通訊順暢，以及僚機之間信息傳遞可達到實時性。最后通過搭建無人機協(xié)同虛擬仿真平臺，驗證所提出的控制方法在穩(wěn)定飛行中的可行性。圖1為無人機協(xié)同編隊飛行控制的流程圖。

圖1 多無人機協(xié)同編隊飛行結(jié)構(gòu)圖

2 模型建立

該無人機長約5 m，翼展寬約2 m，起飛重量約30 kg，包括航電系統(tǒng)設備，PC-104飛行計算機，傳感器，定位接收器，數(shù)據(jù)交換器等。該無人機使用渦輪發(fā)動機，提供最大起飛重量250 N，巡航速度大約為42 m/s，大迎角失速臨界值為30 m/s，最大燃油量為10 L，最大巡航時間為20 min。按“長機-僚機”編隊飛行[20]。無人機協(xié)同編隊幾何模型如圖2所示：

圖2 無人機協(xié)同編隊幾何模型

2.1 無人機縱向氣動力模型

具有推力矢量的連續(xù)性數(shù)學模型參數(shù)，通過系統(tǒng)辨識技術(shù)進行參數(shù)估計。在氣動模型參數(shù)識別研究中，該無人機縱向連續(xù)數(shù)學模型如(1)式和(2)式所示：

0-4.1170.7780

0-33.884-3.5730

0010·

V

α

q

θ+20.168

0.554

-39.085

0iH

(1)

-67.334-4.117-7.9490

20.533-0.655-1.9960

0100·

β

p

r

φ+0.272-0.774

-101.84533.474

-6.261-24.363

00δA

δB

(2)

式中，iH表示失速角;V表示無人機相對于空氣的速度;α,q和θ分別表示無人機的迎角,俯仰角速度以及俯仰角;δA和δB為舵偏角。

該數(shù)學模型描述無人機編隊縱向機動以恒定速度Vo=42 m/s飛行,飛行高度為500 m,αo≈3°。這種模型為單機模型,是建立在無人機非線性模型簡化后(小擾動假設后的線性化后)的模型。

2.2 無人機非線性模型

無人機非線性動態(tài)系統(tǒng)模型的控制理論設計分析是比較困難的,大多數(shù)非線性系統(tǒng)模型參數(shù)采用最小二程法估計。無人機非線性模型如下[21-23]:

=f(x,δ,G,FA(x,δ),MA(x,δ))

(3)

式中，x為無人機狀態(tài)矢量,y和δ均為輸出矢量,G為慣性坐標系下的幾何參數(shù)矢量,f和y為剛體坐標系下的函數(shù)模型。氣動力FA和氣動力矩MA通常用氣動參數(shù)表示:CD為阻力系數(shù)、CY為側(cè)向力系數(shù)、CL為升力系數(shù)、Cl為滾轉(zhuǎn)力矩系數(shù)、Cm為俯仰力矩系數(shù)、Cn為偏航力矩系數(shù)。

FA=SCD(x,δ)

CY(x,δ)

CL(x,δ)MA=SCl(x,δ)

Cm(x,δ)

Cn(x,δ)

3 控制器設計

編隊飛行控制包括兩方面的問題:縱向和橫向控制器設計,本文針對縱向控制器設計。在慣性坐標系下“長機-僚機”的相對位置如圖3所示:

圖3 “長機-僚機”編隊飛行位置幾何模型

3.1 內(nèi)環(huán)控制設計

內(nèi)環(huán)控制設計是為了盡可能的增大反饋增益實現(xiàn)預期的航跡,同時保持較合理的航姿,它基于(1)式和(2)式?？v向內(nèi)環(huán)控制器是對期望的滾轉(zhuǎn)角來設計的,它被擴展到外環(huán)控制器,有以下關(guān)系:

iH=Kqq+Kθq(θ-θd)

(4)

這些參數(shù)通過根軌跡的方法辨識。Kq=0.12,Kθ=0.50,允許閉環(huán)系統(tǒng)的阻尼比ξ=0.54,閉環(huán)系統(tǒng)增益裕度GM=15.98 dB,相裕度PM=870。

由此產(chǎn)生的橫側(cè)向內(nèi)環(huán)控制器是一種線性控制器用來追蹤所設定的迎角,通過外環(huán)控制器來控制,增加無人機橫側(cè)向穩(wěn)定,使用以下的關(guān)系:

δA=Kpp+Kφ(φ-φd)

(5)

δR(s)=Krss+ωor(s)

(6)

根據(jù)根軌跡方法,滾轉(zhuǎn)角反饋增益Kp=0.04,Kφ=0.35,使得閉環(huán)增益ξ=0.35,閉環(huán)系統(tǒng)增益裕度GM=13.4 dB,相裕度PM=78°,偏航角反饋增益Kr=0.16。機翼負扭矩常數(shù)ωo=1.8,閉環(huán)滾轉(zhuǎn)阻尼ξ=0.70,增益裕度GM=20.30 dB,相裕度PM=95.1°。

3.2 外環(huán)控制器設計

外環(huán)控制器包括兩部分,垂直和水平兩部分。垂直方向是一個線性控制器,通過控制高度誤差和導數(shù),調(diào)節(jié)飛機的操穩(wěn)性。控制器提供預期的螺旋角θd輸入內(nèi)環(huán)控制器中:

θd=Khh+Khs

(7)

水平是一種非線性動態(tài)反演控制器,基于控制器使用前向和后向誤差作為輸入,該控制器要求滾轉(zhuǎn)角和迎角指令作為內(nèi)環(huán)控制器指標:

δT

θd=f(x-xL,f

L)

(8)

非線性動態(tài)反演控制方法,在平面動力學范疇取消非線性,通過使用反饋線性化模型。在一些特定的假設條件下,這種控制方法可以使用一種簡單的線性技術(shù)來實現(xiàn)。對于特定問題,非線性系統(tǒng)的特征是以前向和后向距離誤差f和l作為輸出,設定的迎角θd和滾轉(zhuǎn)角δT作為輸入。根據(jù)反饋線性化理論,直到輸入有確定的信號后,才有輸出。輸入和輸出導數(shù)是可逆的,針對此問題,一階和二階導數(shù)輸出如下:

VLxy-Vxy(cos(x-xL))+ΩLf

-l

(9)

Vxysin(x-xL)-VxyVω1cos(x-xL)×

qVtan(φd)

Tb+KTδT+VxyVω2-sin(x-xL)

cos(x-xL)-

ΩLVxycos(x-xL)

sin(x-xL)+Lf

l+ΩL

(10)

式中

ω1=1mcosαcosβ

ω2=Sm(CDcosβ-CYsinβ)+gsinγ

編隊在飛行過程中,協(xié)同轉(zhuǎn)彎的條件為:

Ω=?=gVtanφ

(11)

由此2×2矩陣相關(guān)的輸入和二階導數(shù)的輸出是可逆的,由此產(chǎn)生的逆關(guān)系:

qVtan(φd)

T+KTδT=

1Vxycos(x-xL)sin(x-xL)

Vω1sin(x-xL)-Vω1sin(x-xL)d

ΩL

ω2ω1+lsin(x-xL)-cos(x-xL)

-Vω1lcos(x-xL)-Vω1lsin(x-xL)×

ΩLVxy+lsin(x-xL)-fcos(x-xL)

-Vω1lcos(x-xL)-Vω1lsin(x-xL)×LVxy

(12)

φd=arctan

1gcosγ[dcos(x-xL)+dsin(x-xL)]+VgΩL

(13)

非線性動力學前饋通道描述為:

δT=mKTcosγ[dsin(x-xL)-dcos(x-xL)]+

1KT12ρoV2S(CDo+CDααo)+msinγ-Tb-

mKTcosγΩL[cos(x-xL)+sin(x-xL)]

(14)

輸入信號為δTφT,忽略非線性特性,得出以下關(guān)系:

(15)

上述系統(tǒng)本質(zhì)是一個線性系統(tǒng)組成的2個通道,每個通道都有集成串聯(lián)的通道,這種系統(tǒng)可使用線性控制補償器,對于這種特殊的控制器,設計如下控制律:

d=-ls-ll

(16)

選擇以上的增益以便與之前設計的控制器等效,由此產(chǎn)生3個通道的增益值為:

Kh=3.225,Khs=1.759

(17)

4 仿真分析

4.1 無人機飛行動態(tài)特性仿真

無人機編隊飛行過程中,仿真初始條件設置為:H=2 000 m;V=100 m/s;α=8°;γ=0°;δe=4°;δz=0°;T=50 kN[24];Iy=205 130 kgm2;m=15 000 kg;S=37.14 m2;=3.35 m;xT=7.2 m;-25°≤δe≤25°;-20°≤δz≤20°。仿真曲線如圖4所示:

由圖4a)～f)可知,無人機以100 km/h的速度進入機動,隨著迎角的增加,無人機高度增加,速度下降很快。約5.3 s時,無人機迎角達到55°,遠遠超過常規(guī)無人機機動中最大限制迎角30°,無人機速度降到最小,約為42 km/h,高度增加到2.07 km。此后操控無人機減小迎角,減小速度,無人機高度下降,速度增加,約在8.3 s時,無人機迎角改變44°,速度恢復到58 km/h,隨后無人機逐漸增大速度,恢復到正常穩(wěn)定飛行狀態(tài),與實際機動基本接近。

在大迎角機動過程中,平尾偏角、推力矢量舵偏角均未超過23°。剛開始,無人機速度逐漸增大,平尾偏角和推力矢量舵偏角緩慢增大,3 s時,推力矢量舵偏角達到最大30°,此時無人機速度也達到最大,處于臨界失速狀態(tài)。7 s后,無人機做大機動,平尾偏角和推力矢量舵偏角處于動態(tài)變化,最后逐漸趨于0°,說明無人機在機動過程中保持良好的動態(tài)穩(wěn)定性。俯仰角速度q達到最大至30°/s,在剛開始的1 s,無人機處于水平狀態(tài),俯仰角速度為0,隨著速度增大,3.2 s時達到最大;之后逐漸減小,7 s時速度最小,無人機處于失速狀態(tài),7 s之后,速度逐漸恢復到穩(wěn)定飛行狀態(tài),俯仰角速度也逐漸趨向于零。此過程說明無人機在機動過程中能夠保持航向穩(wěn)定。

圖4 無人機編隊飛行過程姿態(tài)變化動態(tài)特性

仿真結(jié)果顯示,無人機在做大迎角縱向機動過程中,隨著俯仰角、迎角的增大,飛行速度降低,當達到最大迎角時,速度從100 m/s降到了42 km/s,在7 s時俯仰角和迎角恢復,速度也開始增加?？梢?現(xiàn)代無人機在縱向大迎角機動過程中能量損失非常大,損失率約75.93%,此時迎角必須及時恢復到穩(wěn)定飛行狀態(tài)時的大小,因此現(xiàn)代無人機對發(fā)動機性能提出了更高要求,要求發(fā)動機具有較大的推重比,同時協(xié)同飛行控制系統(tǒng)的設計也同樣重要。

4.2 編隊飛行動態(tài)特性分析

無人機編隊穩(wěn)定性控制器在設計階段是非常關(guān)鍵的,但仿真研究表明閉環(huán)系統(tǒng)穩(wěn)定性和追蹤效果,是由測量和模型誤差所引起的。下面對編隊協(xié)同飛行進行穩(wěn)定性進行分析:

其中編隊:lc=-20 m,fc=20 m,hc=20 m

仿真初始數(shù)據(jù)[20]:

CDo=0.009,CDα=0.508,CDq=0

CDiH=-0.034,CLo=-0.049,CLα=3.258

CLq=-0.001,CLiH=0.190,Cmo=0.027

CDα=-0.474,Cmq=-3.449,CmiH=-0.364

(18)

本文所設計的控制對于穩(wěn)態(tài)誤差具有很好的魯棒性。對于閉環(huán)系統(tǒng)穩(wěn)定評估,在允差±5%范圍內(nèi)選取6個參數(shù),分別為τTd,CLo,CLα,Cmα,Clα和Cnβ。在以下仿真圖中,仿真結(jié)果和實際飛行數(shù)據(jù)進行比較,尤其“長機”的飛行參數(shù)對于“僚機”作為參考,僚機仿真的初始條件設置相同,此時編隊控制幾乎接近真實僚機。仿真圖如圖5所示:

圖5 無人機編隊飛行誤差仿真圖

在飛行過程,長機和僚機之間的垂直高度誤差,呈現(xiàn)一種顫振趨勢,之后趨于穩(wěn)態(tài)。在剛開始,長機和僚機之間存在信息傳遞堵塞,曲線的振蕩幅度較大,而后趨于固定值。在無人機編隊執(zhí)行任務過程,為防止無人機之間因距離過小而發(fā)生碰撞,所設計的控制器使無人機之間的橫側(cè)向距離誤差處于安全范圍內(nèi)某一固定值。由圖5b)可知,在剛開始,長機和僚機之間的橫側(cè)向距離不一致,呈現(xiàn)一種偏差幅度較大的趨勢,通過反饋調(diào)節(jié)偏差逐漸減小,而后緩慢趨于穩(wěn)定。在無人機編隊協(xié)同飛行過程中,前向反饋和后向反饋調(diào)節(jié)對無人機的穩(wěn)定性起著重要的作用,它可以保持無人機編隊幾何隊形不變,同時每架無人機具有相同位置和姿態(tài)。如圖5c)可知,在剛開始,后向反饋偏振幅值較小,同時后向反饋比前向反饋能夠在較短時間達到姿態(tài)穩(wěn)定。在無人機近距編隊控制過程中,所設計的閉環(huán)系統(tǒng)控制器,能夠使無人機之間的相對誤差呈現(xiàn)一種先增大后逐漸減小,最終緩慢趨于穩(wěn)定的趨勢。由圖5d)可知,編隊飛行前向反饋和后向反饋,在剛開始出現(xiàn)小幅度的振蕩,這是由于攜有雷達的“長機”從接受“僚機”距離誤差信號,到進行信號處理然后分配僚機,在此期間,會出現(xiàn)信息傳遞堵塞現(xiàn)象,則“長機-僚機”和“僚機-僚機”出現(xiàn)協(xié)同飛行中距離不一致。表1為文獻[20]和本文提出高速近距控制器仿真數(shù)據(jù)對比。

表1 仿真數(shù)據(jù)對比分析

從文獻仿真數(shù)據(jù)和本文仿真數(shù)據(jù)對比可知：在無人機高速近距重編隊過程中，平均前向距離誤差，平均后向距離誤差以及平均垂直距離誤差均比文獻仿真誤差小，而且前向距離誤差減小了36.217%，后向距離誤差減小了5.2%，以及平均垂直距離誤差減小了47.155%，說明本文所設計的控制器可使得無人機編隊協(xié)同飛行過程中，能夠快速近距編隊重組，完成作戰(zhàn)任務。

5 結(jié) 論

本文所設計的控制器是基于非線性動態(tài)反演對非線性數(shù)學模型進行控制，通過內(nèi)外環(huán)協(xié)同控制，使得無人機編隊從集結(jié)隊形，到松散編隊再到編隊集結(jié)后收斂于穩(wěn)定飛行狀態(tài)，而且具有良好的魯棒性。通過選取合適的反饋增益，無人機能夠按照預期的軌跡和幾何編隊飛行，但對近距離或超近距離協(xié)同飛行而言，內(nèi)環(huán)控制器在靈敏度稍微有些遲緩，下一步需要繼續(xù)進行研究。

參考文獻：

[1] 樊瓊劍,楊忠,方挺,等. 多無人機協(xié)同編隊飛行控制的研究現(xiàn)狀[J]. 航空學報, 2009,4(4): 0683-0691

Fan Qiongjian, Yang Zhong,Fang Ting, et al. Research Status of Coordinated Formation Flight Control for Multi-UAVs[J]. Acta Aeronautica et Astronautica Sinica, 2009,4(4): 0683-0691 (in Chinese)

[2] 汪朝群. 紅外誘餌對紅外制導導彈的干擾特性及仿真[J]. 紅外與激光工程, 2001,30(4): 163-167

Wang Chaoqun. Some Characteristics of Infrared Jam and Its Simulation Technique on Infrared Guided Missile[J]. Infrared and Laser Engineering, 2001,30(4): 163-167 (in Chinese)

[3] 樊瓊劍, 楊忠, 馮茂巖. 一種基于雙環(huán)控制的多無人機緊密編隊飛行[C]∥2008年航空宇航科學與技術(shù)全國博士生學術(shù)論壇, 2008

Fan Qiongjian, Yang Zhong, Feng Maoyan, et al. A Multi-UAV Dased on Dual-Loop Control[C]∥Nation Ph D Academic Forum, 2008 (in Chinese)

[4] Corey J S, Rajeeva K. Adaptive Control of UAVs in Close-Coupled Formation Flight[C]∥Proceedings of the American Control Conference, 2000: 849-853

[5] Hummel D. The Use of Aircraft Wakes to Achieve Power Reduction in Formation Flight[C]∥Proceedings of the Fluid Dynamics Panel Symposium, 1996: 1777-1794

[6] Chichka D F, Speyer J, Park C G. Peak-Seeking Control with Plication to Formation Flight[C]∥Proceedings of the 38th IEEE Conference on Decision and Control, 1999: 2463-2470

[7] Lee G, Chong N Y. Decentralized Formation Control for Small-Scale Robot Teams with Anonymity[J]. Mechatronics,2009, 19(1): 85-105

[8] Hu A H, Cao J D, Hu M F, et al. Consensus of a Leader-Following Multi-Agent System with Negative Weights and Noises[J]. IET Control Theory & Applications, 2014, 8(2): 114-119

[9] Ablikim M, Achasov M N, Albayrak O, et al. Observation of a Charged Charmoniumlike Structure in e+e-→(D*D*)±π at

s=4.26 GeV[J]. Physical Review Letters, 2014, 112(13):132001

[10] Xie Y Y, Wang Y, Ma Z J. Delay Consensus of Leader-Following Multi-Agent Systems[J]. Acta Physicalsinica, 2014, 63(4): 40202-040202

[11] Ji L H, Liao X F. Consensus Analysis of Multi-Agent System with Multiple Time Delays[J]. Acta Physical Sinica, 2012, 61(15): 415-418

[12] Nan Jing. Solution of a Class of Rotational Relativistic Rotation Dynamic Equation Using the Generalized Variational Iteration Theory[J]. Acta Physica Sinica, 2013, 62(4): 221-229

[13] An B R, Liu G P. Using Time Delay Compensation Scheme in Dynamic Output Feedback Controller for Networked Control Systems[J]. Acta Physica Sinica, 2014, 63(9): 972

[14] Olfati-Saber R. Flocking for Multi-Agent Dynamic Systems: Algorithms and Theory[J]. IEEE Trans on Automatic Control, 2006, 51(3): 401-420

[15] Su Z, Bai Z S, Xu J, et al. Synthesis, Structure and Property of Cobalt (Ⅱ) Complexes with 3, 5-di Benzoic Acid[J]. Cryst Eng Comm, 2009, 11(5): 873-880

[16] Oh K K, Ahn H S. Formation Control of Mobile Agents Based on Inter-Agent Distance Dynamics[J]. Automatica,2011, 47(10): 2306-2312

[17] Afkhami A, Soltani-Felehgari F, Madrakian T, et al. Fabrication and Application of a New Modified Electrochemical Sensor Using Nano-Silica and a Newly Synthesized Schiff Base for Simultaneous Determination of Cd2+, Cu2+and Hg2+Ions in Water and Some Foodstuff Samples[J]. Analytica Chimica Acta, 2013, 771: 21-30

[18] Ablikim M, Achasov M N, Ai X C, et al. Observation of a Charged Charmoniumlike Structure in e+e-→(D*D*)±π at s=4.26 GeV[J]. Physical Review Letters, 2014, 112(13): 132001

[19] 祁圣君,張喆. 無人機近距編隊飛行建模與仿真[J].航空工程進展,2012,3(3):0366-0367

Qi Shengjun, Zhang Zhe. Modeling and Simulation of UAV Close Formation Flight Control[J]. Advances in Aeronautical Science and Engineering, 2012,3(3): 0366-0367 (in Chinese)

[20] Gu Y, Seanor B, Campa G, et al. Design and Flight Testing Evaluation of Formation Control Laws[J]. IEEE Trans on Control Systems Technology, 2006, 14(6): 1105-1112

[21] Roskam J. Airplane Flight Dynamics and Automatic Flight Controls[J]. Design Analysis & Research Corporation, 1995, 32(3): 0792-0799

[22] Stevens B L, Lewis F L, Johnson E N. Aircraft Control and Simulation: Dynamics, Controls Design, and Autonomous Systems[M]. England, Johnson Wiley, 2015

[23] Rauw M O. FDC 1.2-A Simulink Toolbox for Flight Dynamics and Control Analysis[J]. Gemeente Zeist, 2001, 1(99): 7

[24] Goman M, Khrabrov A. State-Space Representation of Aerodynamic Characteristics of an Aircraft at High Angles of Attack[J]. Journal of Aircraft, 1994, 31(5): 1109-1115