朱科風(fēng) 周慶華 王廣鵬 王樂(lè) 林思宇

文章編號(hào)： 2095-2163（2018）03-0127-05中圖分類號(hào)： 文獻(xiàn)標(biāo)志碼： A

摘要： 關(guān)鍵詞： （School of Physics & Electronic Science， Changsha University of Science & Technology， Changsha 410114， China）

Abstract： The traditional complementary filter algorithm which realizes the data fusion of multi sensors and compensates the integration error of the gyroscope has the advantage of easier calculation. But it also has some defects， such as large data fluctuation and calculating error. In this paper， the differentiation element is introduced to solve the problem of data fluctuation. An improved algorithm is proposed with PID controller to weaken data fluctuation. The experimental results show that the algorithm in this paper has the advantage of less computation and small data fluctuation. It can achieve stable state very quickly and track quadcopter effectively.

Key words：

基金項(xiàng)目：

作者簡(jiǎn)介：

通訊作者：

收稿日期： 引言

科技的不斷發(fā)展使得中小型四軸飛行器（以下簡(jiǎn)稱四軸）的續(xù)航能力越來(lái)越強(qiáng)，四軸的垂直起落、空中懸停、靈活轉(zhuǎn)向等特點(diǎn)也使得其在民用和軍事領(lǐng)域都得到了廣泛的應(yīng)用[1]。而四軸的姿態(tài)獲取是能穩(wěn)定控制四軸飛行的關(guān)鍵。

四軸的姿態(tài)獲取指的是獲取四軸相對(duì)于面的角位置，獲取方式主要依賴于慣性測(cè)量傳感器，其中包括陀螺儀、加速度計(jì)以及磁力計(jì)。但是能搭載在四軸上的傳感器往往精度不高，主要表現(xiàn)在陀螺儀的動(dòng)態(tài)特性會(huì)因積分運(yùn)算而產(chǎn)生累積誤差，且誤差隨時(shí)間增大；加速度計(jì)與磁力計(jì)的靜態(tài)特性容易受到噪聲影響[2]。在硬件條件受限的情況下，對(duì)姿態(tài)解算算法的研究成為了一個(gè)熱點(diǎn)，例如：由Madgwick提出的梯度下降法解決了數(shù)據(jù)融合，避免了地磁補(bǔ)償前對(duì)當(dāng)?shù)氐卮艌?chǎng)角度的測(cè)量，但是在高速運(yùn)動(dòng)情況下偏航角響應(yīng)存在延遲，會(huì)導(dǎo)致姿態(tài)解算出錯(cuò)[3]。Mahony互補(bǔ)濾波算法通過(guò)加速度計(jì)與磁力計(jì)測(cè)量的姿態(tài)數(shù)據(jù)來(lái)補(bǔ)償陀螺儀的積分誤差[4]，運(yùn)算量小，但精度不高。本文采用經(jīng)PID控制器改進(jìn)的互補(bǔ)濾波算法，使姿態(tài)解算響應(yīng)更快，穩(wěn)定性更強(qiáng)，更有利于四軸的控制。

1姿態(tài)描述

為了對(duì)四軸的姿態(tài)進(jìn)行描述，需要建立相應(yīng)的2個(gè)不同的坐標(biāo)系。飛行器的姿態(tài)是參考地面來(lái)描述的，需要以當(dāng)?shù)氐乩砦恢脼閰⒄?，以東、北、天方位為X、Y、Z軸建立導(dǎo)航坐標(biāo)系N[5]。陀螺儀、加速度計(jì)等傳感器是搭載機(jī)體上的，所測(cè)得的數(shù)據(jù)是以機(jī)體為參考，以前、右、上方位為x、y、z軸的載體坐標(biāo)系b來(lái)描述。2個(gè)坐標(biāo)系間的變換關(guān)系則如圖1所示。

圖1導(dǎo)航坐標(biāo)系N與載體坐標(biāo)系b的變換關(guān)系

Fig. 1Transformation relation between coordinate system N and coordinate system b

姿態(tài)的數(shù)學(xué)表示方式通常采用的是歐拉角和四元數(shù)[6]。載體相對(duì)于導(dǎo)航坐標(biāo)系的角位置關(guān)系可以通過(guò)導(dǎo)航坐標(biāo)系繞坐標(biāo)軸連續(xù)旋轉(zhuǎn)3次來(lái)確定：繞導(dǎo)航坐標(biāo)系Z軸旋轉(zhuǎn)φ角，再繞新坐標(biāo)系的Y軸旋轉(zhuǎn)θ角，再繞新坐標(biāo)系的X軸旋轉(zhuǎn)角，實(shí)現(xiàn)變換。φ、θ、分別為載體的偏航角、俯仰角、翻滾角，用歐拉角表示的由N系變換到b系的變換矩陣為CbN，數(shù)學(xué)表述將分別如式（1）、式（2）所示。

CbN=cos φcos θsin φcos θ-sin θ

sin sin θcos φ-sin φcos sin sin θsin φ+cos φcos sin cos θ

sin sin φ+sin θcos φcos sin θsin φcos -sin cos φcos θcos （1）X

Y

Z=CbNx

y

z（2）以歐拉角來(lái)表示變換矩陣直觀易懂，但是歐拉角在求取姿態(tài)時(shí)存在奇點(diǎn)，不能用于全姿態(tài)解算[7]。為了解決求取姿態(tài)時(shí)的奇點(diǎn)問(wèn)題，本文采用四元數(shù)法計(jì)算求取姿態(tài)。四元數(shù)法具有計(jì)算量小，解算無(wú)奇點(diǎn)等優(yōu)點(diǎn)，能夠滿足四軸全姿態(tài)解算的需求。

四元數(shù)由4個(gè)元構(gòu)成，定義q為一個(gè)常規(guī)四元數(shù)向量，數(shù)學(xué)公式具體如下：q=q0+q1i+q2j+q3k（3）其中，q0、q1、q2、q3為實(shí)數(shù)，i、j、k為兩兩正交的向量。當(dāng)僅是關(guān)注N系與b系的角位置關(guān)系時(shí)，兩坐標(biāo)系之間的變換可以等效為無(wú)中間過(guò)程的一次等效旋轉(zhuǎn)，這個(gè)旋轉(zhuǎn)可以由四元數(shù)向量給出表示[8]。矢量v=q1i+q2j+q3k包含了旋轉(zhuǎn)軸和旋轉(zhuǎn)方向的信息，標(biāo)量u=q0包含了旋轉(zhuǎn)角度的信息。根據(jù)四元數(shù)的約束條件，只有單位四元數(shù)才能表示旋轉(zhuǎn)，將四元數(shù)進(jìn)行單位化，研究推得變換公式如下：q=qq20+q21+q22+q23（4）由N系變換到的b系的變換矩陣用四元數(shù)表示為RbN，其計(jì)算公式可見(jiàn)如下：

RbN=

1-2q22+q232q1q2+q0q32q1q3-q0q2

2q1q2-q0q31-2q21+q232q2q3+q0q1

2q1q3+q0q22q2q3-q0q11-2q21+q22（5）

由RbN=CbN得出歐拉角與四元數(shù)關(guān)系，即有如下數(shù)學(xué)公式：θ=sin-1q0q2-q1q3

=tan-12q2q3+q0q11-2q21+q22

φ=tan-12q1q2+q0q31-2q22+q23（6）2改進(jìn)型互補(bǔ)濾波算法

2.1傳統(tǒng)的互補(bǔ)濾波算法原理

四軸通過(guò)陀螺儀、加速度計(jì)以及磁力計(jì)獲得原始的姿態(tài)數(shù)據(jù)。陀螺儀的動(dòng)態(tài)特性好，短時(shí)間內(nèi)測(cè)量的數(shù)據(jù)精度較高，但是由于零漂、溫漂的影響隨著時(shí)間的積累會(huì)形成較大的誤差；而加速度計(jì)與磁力計(jì)靜態(tài)特性好，但動(dòng)態(tài)特性差，容易受高頻噪聲影響。即在頻域上陀螺儀高頻特性好，低頻特性差；加速度計(jì)與磁力計(jì)低頻特性好，高頻特性差。通過(guò)互補(bǔ)的高通、低通濾波器融合3個(gè)傳感器的數(shù)據(jù)，可以實(shí)現(xiàn)較為準(zhǔn)確的姿態(tài)角解算。

根據(jù)傳感器的頻域特性，將陀螺儀的監(jiān)測(cè)數(shù)據(jù)通過(guò)高通濾波器的設(shè)計(jì)處理，同時(shí)將加速度計(jì)與磁力計(jì)的數(shù)據(jù)通過(guò)低通濾波器的綜合處理，2部分?jǐn)?shù)據(jù)相加即得到最終姿態(tài)。互補(bǔ)濾波原理可如圖2所示。

假設(shè)陀螺儀測(cè)量運(yùn)算后的數(shù)據(jù)為θg=θ+μL，加速度計(jì)與磁力計(jì)的測(cè)量數(shù)據(jù)為θma=θ+μH。其中，θ為無(wú)誤差姿態(tài)數(shù)據(jù)，μL為低頻噪聲，μH為高頻噪聲。選取一組傳遞函數(shù)H1s=ss+Cs、H2s=Css+Cs滿足H1s+H2s=1。并且，H1s為高通濾波器傳遞函數(shù)，H2s為低通濾波器傳遞函數(shù)，Cs的選取滿足θs對(duì)互補(bǔ)濾波器全通[9]。濾波器的估計(jì)量頻域響應(yīng)為：

θ^˙s=θgsH1s+θmasH2s=

θs+H1sμLs+H2sμHs（7）

在實(shí)現(xiàn)數(shù)據(jù)互補(bǔ)的同時(shí)為了消除靜態(tài)誤差，選取Cs=Kp+Kis，算法通過(guò)PI控制來(lái)重點(diǎn)調(diào)整濾波器的截止頻率與消除靜態(tài)誤差的時(shí)間。通過(guò)將加速度計(jì)與磁力計(jì)測(cè)得的姿態(tài)角θma與濾波器預(yù)估的姿態(tài)角求差，將差值經(jīng)過(guò)PI調(diào)節(jié)后補(bǔ)償陀螺儀的數(shù)據(jù)得到最終的姿態(tài)角，構(gòu)成一個(gè)如圖2所示的動(dòng)態(tài)的閉環(huán)系統(tǒng)，互補(bǔ)濾波器的估計(jì)量即如式（8）所示：

θ^˙=θg+Kpθma-θ^+Ki∫θma-θ^dt（8）

其中，θ^˙為互補(bǔ)濾波器估計(jì)的輸出量；θ^為反饋量；dt為采樣時(shí)間。

2.2PID改進(jìn)型互補(bǔ)濾波算法

互補(bǔ)濾波算法可以達(dá)到傳感器數(shù)據(jù)的充分融合，通過(guò)運(yùn)用PI控制器，有效地解決了互補(bǔ)濾波器的低通阻帶衰減較慢，受噪聲影響大的問(wèn)題[10- 11]。PI控制器可通過(guò)選取合適的Kp、Ki來(lái)調(diào)整濾波器的截止頻率和消除靜態(tài)誤差的時(shí)間[12]，但是周期性的誤差調(diào)整使得姿態(tài)數(shù)據(jù)不斷地振動(dòng)，導(dǎo)致解算的姿態(tài)角并不穩(wěn)定。四軸應(yīng)對(duì)突變的恢復(fù)能力較弱，不能迅捷恢復(fù)平衡。

為此，針對(duì)姿態(tài)數(shù)據(jù)的振動(dòng)問(wèn)題，本文設(shè)計(jì)引入微分環(huán)節(jié)D，構(gòu)成PID控制器。PID控制器可以通過(guò)微分環(huán)節(jié)反映誤差信號(hào)的變化趨勢(shì)，再據(jù)此展開(kāi)超前調(diào)節(jié)，抑制數(shù)據(jù)的振動(dòng)，并且增強(qiáng)系統(tǒng)的快速性。因此使用PID控制器可以抵消周期調(diào)整誤差所產(chǎn)生的振動(dòng)，還可以加快互補(bǔ)濾波算法的響應(yīng)速度。改進(jìn)算法的流程如圖3所示。

由圖3可知，ab、mb分別為加速度計(jì)與磁力計(jì)的測(cè)量值，ab=axayazT、mb=mxmymzT，歸一化后得ab^、mb^。vb^表示b系重力加速度，wb^表示b系的地磁強(qiáng)度，分別由N系的重力加速度001T與地磁強(qiáng)度bx0bzT經(jīng)過(guò)旋轉(zhuǎn)矩陣變換得到，旋轉(zhuǎn)變換矩陣則為式（5）的RbN。改進(jìn)算法通過(guò)計(jì)算加速度計(jì)、磁力計(jì)的測(cè)量數(shù)據(jù)與重力加速度、地磁強(qiáng)度的誤差e來(lái)補(bǔ)償陀螺儀的數(shù)據(jù)，補(bǔ)償后的陀螺儀數(shù)據(jù)通過(guò)一階龍哥庫(kù)塔方程對(duì)四元數(shù)進(jìn)行更新。各個(gè)參數(shù)的關(guān)系推導(dǎo)可表示如下：e=ab^×vb^+mb^×wb^（9）

經(jīng)過(guò)更新的四元數(shù)通過(guò)反饋環(huán)節(jié)實(shí)現(xiàn)對(duì)重力加速度與地磁強(qiáng)度的變換，再與傳感器數(shù)據(jù)經(jīng)過(guò)叉乘運(yùn)算求得誤差e，因此誤差e具有上一周期的陀螺儀數(shù)據(jù)信息。實(shí)驗(yàn)中，四軸姿態(tài)信號(hào)采樣周期為0.01 s，陀螺儀數(shù)據(jù)在短時(shí)間內(nèi)變化較小，上一周期的陀螺儀數(shù)據(jù)對(duì)下一個(gè)周期的陀螺儀數(shù)據(jù)具有參考意義。因此經(jīng)過(guò)上一周期四元數(shù)變換后求得的誤差e可以作為下一周期陀螺儀數(shù)據(jù)的預(yù)調(diào)量，構(gòu)成PID控制器的微分部分D。PID控制器的參數(shù)Kp、Ki、Kd需要經(jīng)過(guò)反復(fù)實(shí)驗(yàn)，選取合適的值以滿足較為穩(wěn)定、準(zhǔn)確的姿態(tài)解算。

3實(shí)驗(yàn)驗(yàn)證及結(jié)果分析

本文是基于Arduino101開(kāi)發(fā)板構(gòu)建的四軸硬件平臺(tái)，Arduino101開(kāi)發(fā)板搭載有BMI160六軸慣性測(cè)量單元，包括三軸加速度計(jì)與三軸陀螺儀，磁力計(jì)采用的是MAG3110。測(cè)試數(shù)據(jù)通過(guò)串口傳輸給PC機(jī)，在PC機(jī)上完成數(shù)據(jù)分析，采樣頻率為100 Hz。

陀螺儀與加速度計(jì)的原始數(shù)據(jù)，則如圖4所示。圖4是在小幅連續(xù)振動(dòng)的情況下在陀螺儀與加速度計(jì)上測(cè)得的原始數(shù)據(jù)。由圖4可見(jiàn)，在前3 s連續(xù)振動(dòng)的情況下，陀螺儀數(shù)據(jù)能實(shí)時(shí)反映四軸平臺(tái)的運(yùn)動(dòng)狀態(tài)，動(dòng)態(tài)數(shù)據(jù)較為敏感，可以實(shí)時(shí)跟蹤四軸平臺(tái)的姿態(tài)變化。而加速度計(jì)數(shù)據(jù)不能切實(shí)反映四軸平臺(tái)的運(yùn)動(dòng)狀態(tài)，動(dòng)態(tài)反應(yīng)失真。四軸平臺(tái)回穩(wěn)后，陀螺儀數(shù)據(jù)仍有波動(dòng)，靜態(tài)條件下不穩(wěn)定，有較為明顯的誤差。而加速度計(jì)數(shù)據(jù)靜態(tài)數(shù)據(jù)非常準(zhǔn)確，變化較小。

將本文的四軸測(cè)試平臺(tái)經(jīng)過(guò)多次實(shí)驗(yàn)測(cè)試后可知，改進(jìn)算法的PID參數(shù)設(shè)為Kp=2.0，Ki=0.005，Kd=0.05時(shí)，能獲得較為理想的姿態(tài)解算結(jié)果。將該算法與傳統(tǒng)的互補(bǔ)濾波算法、梯度下降算法對(duì)四軸的俯仰角（Pitch）與翻滾角（Roll）分別給出觀測(cè)的效果對(duì)比。

運(yùn)算得到3種算法的解算的姿態(tài)角結(jié)果，即如圖5所示。

圖5是將測(cè)試平臺(tái)的俯仰角與翻滾角各自在約10°、3°的位置保持幾秒后恢復(fù)到水平靜止?fàn)顟B(tài)所測(cè)得的數(shù)據(jù)。分析后發(fā)現(xiàn)：梯度下降法數(shù)據(jù)波動(dòng)較大，所需回穩(wěn)時(shí)間較長(zhǎng)；互補(bǔ)濾波算法與改進(jìn)型算法波動(dòng)幅度較小，回穩(wěn)速度快，能夠更快地趨于穩(wěn)定。

PID改進(jìn)型算法是在互補(bǔ)濾波算法的基礎(chǔ)上改進(jìn)得來(lái)。為了驗(yàn)證改進(jìn)算法在抑制振動(dòng)上的優(yōu)勢(shì)，分別在俯仰方向與翻滾方向作連續(xù)振動(dòng)測(cè)試，并對(duì)2種算法的運(yùn)算結(jié)果進(jìn)行數(shù)據(jù)采集，在PC上位機(jī)調(diào)取相應(yīng)數(shù)據(jù)分析，可得結(jié)果如圖6所示。

圖6是分別截取2種姿態(tài)角連續(xù)10 s的數(shù)據(jù)。可以看出，改進(jìn)型互補(bǔ)濾波算法運(yùn)行解算的姿態(tài)數(shù)據(jù)波動(dòng)更小，所需趨于平穩(wěn)的時(shí)間更短，而且在回穩(wěn)后數(shù)據(jù)較為穩(wěn)定。相比于互補(bǔ)濾波算法，改進(jìn)型算法的穩(wěn)定特性使得四軸更容易控制。

為了測(cè)試3種算法的靜態(tài)狀況下的解算準(zhǔn)確性，對(duì)3種算法在靜止水平（俯仰角、翻滾角均為零度）狀態(tài)下的姿態(tài)角各取50組求出平均值，最終得到測(cè)量的平均誤差可見(jiàn)表1。從表1可以得出改進(jìn)型互補(bǔ)濾波算法在俯仰角上平均誤差相對(duì)于梯度下降法和互補(bǔ)濾波算法分別降低了45.4%、52.4%；在翻滾角上，平均誤差分別下降了26.3%、62.2%。

結(jié)束語(yǔ)

針對(duì)互補(bǔ)濾波算法存在的數(shù)據(jù)波動(dòng)大等問(wèn)題，本文提出了四元數(shù)PID改進(jìn)型互補(bǔ)濾波算法。該算法通過(guò)四元數(shù)描述四軸姿態(tài)，避免了歐拉角的奇點(diǎn)問(wèn)題，實(shí)現(xiàn)了全姿態(tài)解算，在此基礎(chǔ)上引入了PID控制器調(diào)節(jié)姿態(tài)誤差，使得解算的姿態(tài)角更加穩(wěn)定。后續(xù)又通過(guò)仿真實(shí)驗(yàn)對(duì)比展示了互補(bǔ)濾波算法、梯度下降算法與改進(jìn)型算法的效果及性能。實(shí)驗(yàn)結(jié)果表明，運(yùn)動(dòng)狀態(tài)下改進(jìn)型算法解算出的俯仰角與翻滾角波動(dòng)較小，且所需穩(wěn)定時(shí)間更短；在水平靜止的狀態(tài)下，改進(jìn)型算法的解算誤差最小。由此可以得知，改進(jìn)型算法有利于四軸的控制，能使四軸的飛行更加平穩(wěn)。

參考文獻(xiàn)

[1] 萬(wàn)曉鳳， 康利平， 余運(yùn)俊， 等. 互補(bǔ)濾波算法在四旋翼飛行器姿態(tài)解算中的應(yīng)用[J]. 測(cè)控技術(shù)， 2015，34（2）：8-11.

[2] 陳亮， 楊柳慶， 肖前貴. 基于梯度下降法和互補(bǔ)濾波的航向姿態(tài)參考系統(tǒng)[J]. 電子設(shè)計(jì)工程， 2016，24（24）：38-41，45.

[3] MADGWICK S. An efficient orientation filter for inertial and inertial/magnetic sensor arrays[R]. UK：University of Bristol，2010.

[4] BALDWIN G， MAHONY R， TRUMPF J， et al. Complementary filter design on the Special Euclidean group SE（3）[C]// Proceedings of the European Control Conference. Kos， Greece： IEEE，2007：3763-3770.

[5] 季元揚(yáng)， 陳躍東， 陳孟元. 基于PID改進(jìn)型互補(bǔ)濾波[J]. 控制工程， 2017，24（5）：952-957.

[6] 徐云川. 四軸飛行器姿態(tài)解算算法設(shè)計(jì)與仿真[J]. 科技視界， 2016（23）：17-18.

[7] 陳孟元， 謝義建， 陳躍東. 基于四元數(shù)改進(jìn)型互補(bǔ)濾波的MEMS姿態(tài)解算[J]. 電子測(cè)量與儀器學(xué)報(bào)， 2015，29（9）：1391-1397.

[8] 秦永元. 慣性導(dǎo)航[M]. 2版. 北京：科學(xué)出版社有限責(zé)任公司，2014.

[9] MAHONY R， HAMEL T， PFLIMLIN J. Nonlinear complementary filters on the special orthogonal group[J]. IEEE Transactions on Automatic Control， 2008，53（5）：1203-1218.

[10]呂印新， 肖前貴， 胡壽松. 基于四元數(shù)互補(bǔ)濾波的無(wú)人機(jī)姿態(tài)解算[J]. 燕山大學(xué)學(xué)報(bào)， 2014，38（2）：175-180.

[11]杜杉杉， 吳昊， 張繼文， 等. 一種面向AHRS的改進(jìn)互補(bǔ)濾波融合算法[J]. 國(guó)外電子測(cè)量技術(shù)， 2015，34（3）：13-18.

[12]陳雷， 卞鴻巍. 基于模糊自適應(yīng)互補(bǔ)濾波的姿態(tài)解算算法[J]. 兵器裝備工程學(xué)報(bào)， 2017，38（4）：183-187.（上接第126頁(yè)）

[4] Yang Xinshe. Flower pollination algorithm for global optimization[C]//LNCS 7445： Proceedings of the 11th Internation Conference on Unconventional Computation and Natural Computation， Orléan， France， Sep 3-7， 2012. Berlin， Heidelberg： Springer， 2012： 240-249.

[5] Cheng S， Zhang Q， Qin Q. Big data analytics with swarm intelligence[J]. Industrial Management & Data Systems， 2016， 116（4）：646-666.

[6] 馮翔，張進(jìn)文，虞慧群. 仿生蚊子追蹤算法[J]. 計(jì)算機(jī)學(xué)報(bào)，2014，37（8）：1794-1808.

[7] 喬瑩瑩，宋威，馬偉. 基于GA優(yōu)化QPSO算法的文本聚類[J]. 計(jì)算機(jī)應(yīng)用研究，2014，31（10）：2912-2915.

[8] Yang X S. A New Meta-heuristic Bat-Inspired Algorithm[J]. Computer Knowledge & Technology， 2010， 284：65-74.

[9] Wang G G， Chu H C E， Mirjalili S. Three-dimensional path planning for UCAV using an improved bat algorithm[J]. Aerospace Science & Technology， 2016， 49：231-238.

[10]Karri C， Jena U. Fast vector quantization using a Bat algorithm for image compression[J]. Engineering Science & Technology An International Journal， 2016， 19（2）：769-781.

[11]姚妮，李紅嬋. 基于混合蝙蝠算法的多目標(biāo)柔性作業(yè)車間調(diào)度問(wèn)題[J]. 微電子學(xué)與計(jì)算機(jī)，2017，34（3）：25-29，34.

[12]肖輝輝，段艷明. 基于DE算法改進(jìn)的蝙蝠算法的研究及應(yīng)用[J]. 計(jì)算機(jī)仿真，2014，31（1）：272-277.

[13]He X S， Ding W J， Yang X S. Bat algorithm based on simulated annealing and Gaussian perturbations[J]. Neural Computing & Applications， 2014， 25（2）：459-468.

[14]裴宇航，劉景森，李煜. 一種動(dòng)態(tài)調(diào)整慣性權(quán)重的自適應(yīng)蝙蝠算法[J]. 計(jì)算機(jī)科學(xué)，2017，44（6）：240-244.

[15]Mirjalili S. SCA： A Sine Cosine Algorithm for solving optimization problems[J]. Knowledge-Based Systems， 2016， 96：120-133.

[16]王文， 王勇， 王曉偉. 采用機(jī)動(dòng)飛行的蝙蝠算法[J]. 計(jì)算機(jī)應(yīng)用研究， 2014， 31（10）：2962-2964.