嵌入式大氣數(shù)據(jù)傳感系統(tǒng)的模糊邏輯建模方法

王 坤,黃 達(dá)

(南京航空航天大學(xué) 航空學(xué)院,南京 210016)

0 引 言

現(xiàn)代飛行器的發(fā)展對(duì)高速、隱身性能的要求越來越高。傳統(tǒng)的安裝在飛行器前端的空速管等探針式大氣數(shù)據(jù)傳感器對(duì)飛行器的隱身、機(jī)動(dòng)性能的不利影響較大。嵌入式大氣數(shù)據(jù)傳感系統(tǒng)(Flush Air Data Sensing System,FADS System)是依靠嵌入在飛行器頭部的若干壓力傳感器來獲得飛行器表面壓力分布，并由此壓力分布間接得出飛行參數(shù)的系統(tǒng)，曾用于X-33、X-34和X-38等飛行器上[1]。由于壓力傳感器分布與飛行器表面平齊，相比傳統(tǒng)的探針式傳感器系統(tǒng)，這種嵌入式傳感系統(tǒng)在大馬赫數(shù)、大迎角下依然能夠很好地工作，且有利于提高飛行器隱身性能。

早期對(duì)FADS系統(tǒng)的研究主要針對(duì)球形頭部的飛行器，提出了基于圓球繞流理論解的三點(diǎn)法[1-3]。后來發(fā)展了針對(duì)該系統(tǒng)的神經(jīng)網(wǎng)絡(luò)算法[4-9]，通過大量的試驗(yàn)數(shù)據(jù)訓(xùn)練，適用于非線性系統(tǒng)建模。

模糊邏輯建模方法以往多用于對(duì)非定常氣動(dòng)力的建模[10-12]或者工業(yè)控制領(lǐng)域。本文研究了模糊邏輯對(duì)FADS系統(tǒng)的建模。以模型表面若干測(cè)壓點(diǎn)的壓力或壓力系數(shù)作為模糊系統(tǒng)的輸入，以迎角、側(cè)滑角、來流速度和海拔高度作為輸出，混合使用梯度下降法和最小二乘估計(jì)法來識(shí)別模糊系統(tǒng)的參數(shù)，從而建立針對(duì)該系統(tǒng)的模糊邏輯模型。該模型建立完成之后，將它運(yùn)用于飛行器的傳感系統(tǒng)，可以實(shí)時(shí)地根據(jù)若干測(cè)壓點(diǎn)的壓力得到飛行器的飛行狀態(tài)參數(shù)，然后傳遞給控制系統(tǒng)。

1 風(fēng)洞試驗(yàn)及CFD仿真

1.1 風(fēng)洞試驗(yàn)部分

試驗(yàn)在南京航空航天大學(xué)NH-2低速風(fēng)洞進(jìn)行。風(fēng)洞試驗(yàn)段截面為3.0 m×2.5 m的矩形，試驗(yàn)?zāi)Ｐ蜑槟筹w翼布局飛行器，表面共布置62個(gè)測(cè)壓點(diǎn)，頭部分布46個(gè)測(cè)壓點(diǎn)，內(nèi)翼左右各分布8個(gè)測(cè)壓點(diǎn)，位置基本對(duì)稱(圖1、圖2)。測(cè)壓采用內(nèi)埋測(cè)壓管的方式，測(cè)壓管內(nèi)徑0.3 mm。每個(gè)測(cè)壓管均與壓力測(cè)量模塊相連。

圖1 風(fēng)洞試驗(yàn)?zāi)Ｐ虵ig.1 Model in the wind tunnel

圖2 模型幾何外形及測(cè)壓點(diǎn)分布Fig.2 Model’s geometric shape and pressure tap positions

圖3中給出了半側(cè)機(jī)身上的測(cè)壓點(diǎn)分布，除1～8號(hào)測(cè)壓孔外，另半側(cè)機(jī)身上的測(cè)壓孔分布與圖中“Z”開頭編號(hào)的測(cè)壓孔對(duì)稱，相對(duì)應(yīng)的編號(hào)以“Y”開頭，例如“Z18”對(duì)應(yīng)“Y18”。

機(jī)頭為坐標(biāo)系原點(diǎn)。x軸位于飛機(jī)對(duì)稱面指向向后，與機(jī)身水平線平行；y軸向上為正；z軸指向左機(jī)翼為正。模型展長(zhǎng)3.5 m，部分關(guān)鍵測(cè)壓點(diǎn)坐標(biāo)如表1所示。

圖3 測(cè)壓孔編號(hào)及位置分布(前視圖視角)Fig.3 Pressure tap NOs.and positions (front view)

表1 部分測(cè)壓點(diǎn)位置(單位：mm)Table 1 Some pressure tap positions (unit:mm)

試驗(yàn)采用DTC Initium電子掃描系統(tǒng)采集表面壓力。該系統(tǒng)有一臺(tái)主控制器，可帶8個(gè)測(cè)壓模塊，每個(gè)模塊可連接64個(gè)測(cè)壓通道，數(shù)據(jù)采樣頻率330 Hz，系統(tǒng)精度0.1%。試驗(yàn)采用1個(gè)量程為1 psi的測(cè)壓模塊測(cè)量模型表面壓力。

試驗(yàn)風(fēng)速為35 m/s、55 m/s和70 m/s，試驗(yàn)迎角-18°～24°，每隔2°進(jìn)行試驗(yàn)，側(cè)滑角0°～16°，每隔4°進(jìn)行試驗(yàn)。最終得到在不同迎角、側(cè)滑角和風(fēng)速情況下的模型表面62個(gè)測(cè)壓點(diǎn)的壓力和壓力系數(shù)值。

由于地面風(fēng)洞試驗(yàn)不能模擬不同海拔高度的情況，所以輔助使用CFD仿真來計(jì)算在不同海拔高度下的流場(chǎng)。

1.2 CFD仿真部分

流場(chǎng)網(wǎng)格使用四面體/棱柱體混合網(wǎng)格，使用Fluent求解器，湍流模型采用S-A模型。來流速度為35 m/s、70 m/s，迎角為-10°、0°、10°，側(cè)滑角為0°、8°，海拔高度為500 m、1000 m、1500 m、2000 m、3000 m、4000 m、5000 m和6000 m共8個(gè)高度。一共96個(gè)計(jì)算狀態(tài)。CFD計(jì)算得到的壓力數(shù)據(jù)用于對(duì)海拔高度的建模。

2 模糊邏輯建模

2.1 模糊結(jié)構(gòu)

(1)

下面以2個(gè)輸入量為例，說明該模糊結(jié)構(gòu)的建立。

圖4以5層為例展示模糊結(jié)構(gòu)。

圖4 模糊結(jié)構(gòu)示意圖Fig.4 Fuzzy model structure

第1層表示對(duì)于兩個(gè)輸入量x1、x2各分配若干個(gè)隸屬函數(shù)，分別是A1、A2、B1、B2、B3。

第2層表示對(duì)分別屬于各輸入變量的隸屬函數(shù)進(jìn)行交叉乘積計(jì)算，得到內(nèi)部函數(shù)的權(quán)重值，分別是w1、w2、w3、w4、w5、w6。

(2)

第4層表示將權(quán)重系數(shù)與內(nèi)部函數(shù)相乘。

第5層表示將第4層的輸出進(jìn)行求和運(yùn)算，得到最終輸出y。

(3)

2.2 隸屬函數(shù)

識(shí)別隸屬函數(shù)參數(shù)時(shí)，根據(jù)不同的識(shí)別方法，分別使用了三角形[10]、梯形和鐘形[14]隸屬函數(shù)。本文使用到的3種不同形狀的隸屬函數(shù)示意圖如圖5～圖7所示。

鐘形隸屬函數(shù)表達(dá)式：

(4)

此函數(shù)由3個(gè)參數(shù)ai、bi、ci決定，在進(jìn)行參數(shù)識(shí)別的時(shí)候利用梯度下降法動(dòng)態(tài)調(diào)整這3個(gè)參數(shù)。

圖5 三角隸屬函數(shù)示意圖Fig.5 Triangular membership function

圖6 梯形隸屬函數(shù)示意圖Fig.6 Trapezoid membership function

圖7 鐘形隸屬函數(shù)示意圖Fig.7 Bell-shaped membership function

2.3 參數(shù)識(shí)別

2.3.1 隸屬函數(shù)參數(shù)識(shí)別

(5)

式中a,b,c,…表示該節(jié)點(diǎn)包含的參數(shù)。假設(shè)共有m組樣本數(shù)據(jù)，對(duì)于第j組樣本數(shù)據(jù)，定義對(duì)某組樣本數(shù)據(jù)的平方差Ej：

(6)

為求解E對(duì)隸屬函數(shù)參數(shù)的偏導(dǎo)數(shù)，首先求解?Ej/?O。對(duì)于第L層結(jié)構(gòu)可以根據(jù)Ej的表達(dá)式直接推導(dǎo)出來。

(7)

(8)

根據(jù)式(7)和式(8)，可以求出對(duì)任意節(jié)點(diǎn)處的?Ej/?O。

現(xiàn)在來求解Ej對(duì)隸屬函數(shù)參數(shù)的偏導(dǎo)數(shù)。假設(shè)α是其中一個(gè)隸屬函數(shù)的某個(gè)參數(shù)。

(9)

S表示節(jié)點(diǎn)輸出依賴于α的節(jié)點(diǎn)集合。因?yàn)橛衜組樣本數(shù)據(jù)，則：

(10)

(11)

k表示在參數(shù)空間內(nèi)梯度變化的步長(zhǎng)。在編程計(jì)算時(shí)可以通過改變k來加快收斂。例如當(dāng)平方差之和E連續(xù)四次迭代降低時(shí)，可以將步長(zhǎng)變?yōu)?.25k，當(dāng)E出現(xiàn)增加減少來回波動(dòng)的情況時(shí)，將步長(zhǎng)變?yōu)?.8k。

2.3.2 內(nèi)部函數(shù)參數(shù)識(shí)別

可以使用梯度下降法或者最小二乘估計(jì)法來識(shí)別內(nèi)部函數(shù)參數(shù)。

模糊輸出：

(12)

2.3.2.1 梯度下降法

模糊輸出與樣本輸出的平方差之和：

(13)

m是樣本個(gè)數(shù)。定義相關(guān)系數(shù)：

(14)

待定的內(nèi)部函數(shù)參數(shù)的個(gè)數(shù)是N1×N2×…×Nr×…×Nk×(k+1)。

目標(biāo)模糊邏輯結(jié)構(gòu)要使得E最小，使用梯度下降法更新內(nèi)部函數(shù)參數(shù)。

(15)

αr是收斂因子，在(0,1)之間。當(dāng)αr選取的較小時(shí)，建模速度慢，當(dāng)αr選取過大時(shí)，計(jì)算出現(xiàn)溢出。需要不斷調(diào)整選取合適的值，最大化建模速度。

將E的表達(dá)式代入展開，得到：

當(dāng)r=0時(shí)，

(16)

當(dāng)r=1,2,3,…,k時(shí)，

(17)

2.3.2.2 最小二乘估計(jì)法

模糊結(jié)構(gòu)可以表示為AP=B。

P表示內(nèi)部函數(shù)參數(shù)構(gòu)成的列向量矩陣，元素個(gè)數(shù)為(k+1)×N1×…×Nr×…×Nk。A矩陣由樣本輸入和隸屬函數(shù)計(jì)算得到，行數(shù)為樣本數(shù)m，列數(shù)為(k+1)×N1×…×Nr×…×Nk。B矩陣表示樣本輸出，元素個(gè)數(shù)為樣本數(shù)m。對(duì)此一階線性方程組，通常使用最小二乘法求解。

將P矩陣的最小二乘估計(jì)記做P*，則：

P*=(ATA)-1ATB

(18)

ATA是非奇異矩陣時(shí)，它的逆存在。當(dāng)矩陣維數(shù)較大時(shí)，直接計(jì)算逆矩陣的效率很低。如果樣本數(shù)量m<(k+1)×N1×…×Nr×…×Nk，上述公式無法使用。下面使用卡爾曼濾波算法來計(jì)算P*[13-14]。

(19)

(20)

2.4 模糊結(jié)構(gòu)識(shí)別

此模糊結(jié)構(gòu)的輸入是62個(gè)測(cè)壓點(diǎn)的壓力或壓力系數(shù)，選取其中k個(gè)點(diǎn)的值作為輸入，輸出為迎角、側(cè)滑角、速度和高度。

將輸入記做x1,…,xr,…,xk，將輸出記做yα,yβ,yv,yh。以對(duì)迎角建模為例，對(duì)側(cè)滑角、來流速度和海拔高度的建模相同。

(2)初始時(shí)，對(duì)每個(gè)輸入選取Nbegin個(gè)隸屬函數(shù)。則初始隸屬函數(shù)個(gè)數(shù)的結(jié)構(gòu)是(Nbegin,…,Nbegin,…,Nbegin)。然后對(duì)其中某一個(gè)輸入的隸屬函數(shù)個(gè)數(shù)加1，其他輸入的隸屬函數(shù)個(gè)數(shù)保持不變，這樣可以建立k個(gè)隸屬函數(shù)個(gè)數(shù)不同的結(jié)構(gòu)。

(5)得到最優(yōu)隸屬函數(shù)個(gè)數(shù)的結(jié)構(gòu)之后，用所有樣本數(shù)據(jù)重新進(jìn)行參數(shù)識(shí)別，得到最終模糊結(jié)構(gòu)參數(shù)，建模結(jié)束。

3 建模結(jié)果分析

試驗(yàn)?zāi)Ｐ捅砻婀膊贾?2個(gè)測(cè)壓點(diǎn)。選取其中若干個(gè)測(cè)壓點(diǎn)的壓力或壓力系數(shù)作為輸入，以迎角、側(cè)滑角、風(fēng)速和高度作為輸出。

對(duì)迎角和側(cè)滑角建模時(shí)選擇壓力系數(shù)作為輸入。試驗(yàn)風(fēng)速最大為70 m/s。CFD計(jì)算時(shí)選取的最大高度為6000 m。在此區(qū)間內(nèi)在姿態(tài)角不變的情況下，風(fēng)速和高度的改變只引起壓力系數(shù)的微小變化，幾乎可以忽略。所以對(duì)風(fēng)速和高度的建模選取壓力作為輸入。

考慮此模型的非線性復(fù)雜程度，通常選擇不少于3個(gè)點(diǎn)作為建模輸入。建模點(diǎn)數(shù)量增加時(shí)，參數(shù)識(shí)別消耗的時(shí)間以指數(shù)程度增加，平衡考慮建模的精度和時(shí)間效益，建模點(diǎn)數(shù)量一般不多于7個(gè)。

從62個(gè)測(cè)壓點(diǎn)選擇若干個(gè)點(diǎn)建模，選取方式很多，可以適當(dāng)探討建模點(diǎn)位置選取與建模精度的關(guān)系，為選點(diǎn)位置的優(yōu)化設(shè)計(jì)提供參考。

3.1 建模結(jié)果

以選取4個(gè)點(diǎn)(2,6,Z4,Y4)建模為例，結(jié)果如下。2、6編號(hào)的點(diǎn)位于模型頭部上下區(qū)域，壓力變化對(duì)迎角的改變敏感，Z4、Y4編號(hào)的點(diǎn)位于頭部左右區(qū)域，壓力變化對(duì)側(cè)滑角的改變敏感。

圖8中橫坐標(biāo)為測(cè)壓孔2的壓力系數(shù)Cp。

表2和圖8～圖12展示了選擇4個(gè)點(diǎn)(2,6,Z4,Y4)建模的結(jié)果。模糊輸出與試驗(yàn)值的對(duì)比達(dá)到了很高的精度。

圖8 迎角模糊結(jié)構(gòu)輸出與試驗(yàn)值對(duì)比Fig.8 Comparison between modeling results and experimental angle of attack

在對(duì)迎角與側(cè)滑角及風(fēng)速建模時(shí)，此時(shí)有300多組樣本，在樣本較多時(shí)，選用三角隸屬函數(shù)與梯形隸屬函數(shù)差異不大。對(duì)高度進(jìn)行建模時(shí)，樣本較少，有96組。在本文使用的樣本及樣本數(shù)量較少的情況下,梯形隸屬函數(shù)的精度更高。

對(duì)高度建模時(shí)，加入了使用變參數(shù)鐘形隸屬函數(shù)的建模結(jié)果。使用變參數(shù)隸屬函數(shù)建模時(shí)，分別使用了梯度下降法和最小二乘估計(jì)法識(shí)別隸屬函數(shù)參數(shù)與內(nèi)部函數(shù)參數(shù),從表2的數(shù)據(jù)及圖11和圖12的對(duì)比可以看出，精度相比三角隸屬函數(shù)與梯形隸屬函數(shù)高很多，且建模速度很快。使用卡爾曼濾波的最小二乘估計(jì)法更適用于小樣本的情況。圖13和圖14展示了模型訓(xùn)練前后針對(duì)測(cè)壓孔2的隸屬函數(shù)曲線。

圖9 側(cè)滑角模糊結(jié)構(gòu)輸出與 試驗(yàn)值對(duì)比Fig.9 Comparison between modeling results and experimental angle of sideslip

圖10 風(fēng)速模糊結(jié)構(gòu)輸出與 試驗(yàn)值對(duì)比Fig.10 Comparison between modeling results and experimental velocity

圖11 海拔高度模糊結(jié)構(gòu)輸出與 試驗(yàn)值對(duì)比(梯形隸屬函數(shù))Fig.11 Comparison between modeling results and experimental altitude (trapezoid membership function)

圖12 海拔高度模糊結(jié)構(gòu)輸出與 試驗(yàn)值對(duì)比(鐘形變參數(shù)隸屬函數(shù))Fig.12 Comparison between modeling results and experimental altitude (bell-shaped membership function)

圖13 測(cè)壓點(diǎn)2初始隸屬函數(shù)形狀Fig.13 Initial membership function of pressure tap NO.2

圖14 測(cè)壓點(diǎn)2最終隸屬函數(shù)形狀Fig.14 Final membership function of pressure tap NO.2

3.2 建模點(diǎn)數(shù)量與位置的優(yōu)化選擇

以對(duì)迎角建模為例，分別選擇不同位置的3個(gè)點(diǎn)、4個(gè)點(diǎn)或5個(gè)點(diǎn)用于建模。在建模點(diǎn)數(shù)目相同時(shí)，使用相同的收斂因子、迭代次數(shù)和相同的隸屬函數(shù)個(gè)數(shù)和形狀(三角隸屬函數(shù))，記錄精度R2，比較篩選出最優(yōu)的建模位置。

建模結(jié)果如下，表格數(shù)據(jù)以R2降序排序。

表3～表6中選點(diǎn)是編寫程序從62個(gè)測(cè)壓孔中隨機(jī)選取的，精度均很高，相差較小，從R2降序結(jié)果觀察比對(duì)選點(diǎn)的位置區(qū)域，考慮到試驗(yàn)時(shí)使用的測(cè)壓儀器在模型不同區(qū)域的測(cè)量精度不同，通常選取模型

頭部附近的測(cè)壓點(diǎn)進(jìn)行建模。

從表3～表5還可以看出，選擇5個(gè)點(diǎn)的建模精度稍高于4個(gè)點(diǎn)的建模精度，4個(gè)點(diǎn)的建模精度稍高于3個(gè)點(diǎn)的建模精度。不過隨著建模點(diǎn)數(shù)的增加，建模速度大大降低。表中數(shù)據(jù)表明選擇4個(gè)點(diǎn)或5個(gè)點(diǎn)的建模精度能夠很好地滿足工程需求。

計(jì)算時(shí)發(fā)現(xiàn)存在一些精度非常低的選點(diǎn)方式(即處于同一水平面或鉛垂面的點(diǎn))，在實(shí)際應(yīng)用中應(yīng)當(dāng)避免。

表3 建模點(diǎn)數(shù)為3時(shí)對(duì)迎角的建模結(jié)果 (迭代次數(shù)：20000，收斂因子：0.006)Table 3 Alpha outputs with three input points (iterations:20000,convergence factor:0.006)

表4 建模點(diǎn)數(shù)為4時(shí)對(duì)迎角的建模結(jié)果 (迭代次數(shù)：10000，收斂因子：0.006)Table 4 Alpha outputs with four input points (iterations:10000,convergence factor:0.006)

表5 建模點(diǎn)數(shù)為5時(shí)對(duì)迎角的建模結(jié)果 (迭代次數(shù)：10000，收斂因子：0.006)Table 5 Alpha outputs with five input points (iterations:10000,convergence factor:0.006)

表6 低精度建模點(diǎn)選取示例 (迭代次數(shù)：10000，收斂因子：0.006)Table 6 Low accuracy modeling cases (iterations:10000,convergence factor:0.006)

4 結(jié) 論

本文研究了對(duì)嵌入式大氣數(shù)據(jù)傳感系統(tǒng)的模糊邏輯建模方法。以模型表面若干測(cè)壓點(diǎn)的壓力或壓力系數(shù)作為模糊系統(tǒng)的輸入，以迎角、側(cè)滑角、來流速度和海拔高度作為輸出，探討了在使用不同隸屬函數(shù)和參數(shù)識(shí)別方法情況下的建模結(jié)果，討論分析了建模點(diǎn)數(shù)量與位置的優(yōu)化選擇。結(jié)論如下：

1)針對(duì)本文的模糊邏輯系統(tǒng)，使用自適應(yīng)形狀參數(shù)的隸屬函數(shù)的建模精度要高于使用固定形狀的三角隸屬函數(shù)和梯形隸屬函數(shù)。且在建模時(shí)分別使用梯度下降法和最小二乘估計(jì)法來識(shí)別隸屬函數(shù)參數(shù)和內(nèi)部函數(shù)參數(shù)，不論是從精度還是速度上來說都優(yōu)于單獨(dú)使用梯度下降法來識(shí)別內(nèi)部函數(shù)參數(shù)。

2)在建模點(diǎn)數(shù)量與位置的優(yōu)化選擇方面，考慮傳感器在模型不同區(qū)域的測(cè)量精度，通常選擇頭部附近的測(cè)壓數(shù)據(jù)進(jìn)行建模。對(duì)迎角的建模應(yīng)當(dāng)避免選擇處于同一水平面的測(cè)壓點(diǎn)，對(duì)側(cè)滑角避免選擇處于同一鉛垂面的測(cè)壓點(diǎn)。結(jié)果表明選擇4個(gè)點(diǎn)或5個(gè)點(diǎn)時(shí)的精度足以滿足實(shí)際工程需求，且建模速度較快。對(duì)于復(fù)雜外形飛行器，可以適當(dāng)增加建模點(diǎn)數(shù)和隸屬函數(shù)個(gè)數(shù)。

受限于風(fēng)洞試驗(yàn)的高昂成本，本文建立的模糊邏輯模型主要基于低速情況下的流動(dòng)數(shù)據(jù)。想要拓寬模型的應(yīng)用范圍，需要更多的不同流動(dòng)情況下的數(shù)據(jù)，未來的方向是將模型的研究與驗(yàn)證拓展到更高速的流動(dòng)場(chǎng)景下。