李 海 張彬蕾 孟凡旺 孫 研 韓雁飛

(1.中國(guó)民航大學(xué)天津市智能信號(hào)與圖像處理重點(diǎn)實(shí)驗(yàn)室 天津 300300;2.中國(guó)航空工業(yè)集團(tuán)公司雷華電子技術(shù)研究所 江蘇無錫 214063;3.中國(guó)電子科技集團(tuán)航空電子有限公司 成都 611731)

0 引言

當(dāng)民航飛機(jī)在巡航階段(一般在22000英尺以上)飛行時(shí),高空冰晶是對(duì)飛機(jī)飛行安全威脅較大的危險(xiǎn)天氣現(xiàn)象之一。在高空冰晶聚集區(qū)域,發(fā)動(dòng)機(jī)葉片表面的冰晶凝結(jié)會(huì)導(dǎo)致發(fā)動(dòng)機(jī)失控與喘振,燃燒室熄火或發(fā)動(dòng)機(jī)葉片損壞等[1]。高空冰晶極易造成機(jī)載傳感器探頭堵塞,使得探測(cè)數(shù)據(jù)異常,造成飛機(jī)失控[2]。近年來,高濃度冰晶導(dǎo)致噴氣式飛機(jī)發(fā)動(dòng)機(jī)動(dòng)力損失的事件已經(jīng)引起了發(fā)動(dòng)機(jī)制造商、安全監(jiān)管部門和氣象研究界的廣泛注意[3]。據(jù)波音公司統(tǒng)計(jì),自20世紀(jì)90年代初以來,在通勤飛機(jī)和大型運(yùn)輸飛機(jī)中,已有超過150起飛機(jī)動(dòng)力損失或發(fā)動(dòng)機(jī)回滾(發(fā)動(dòng)機(jī)失去控制)等不安全事件被歸因于高空冰晶[4]。這些統(tǒng)計(jì)報(bào)告表明:高空冰晶對(duì)航空運(yùn)輸業(yè)和旅客飛行安全造成了嚴(yán)重影響。因此高空冰晶檢測(cè)及危險(xiǎn)等級(jí)區(qū)域告警成為保障航空安全的一項(xiàng)重要課題。高空冰晶檢測(cè)及其危險(xiǎn)級(jí)別量化通過估計(jì)雷達(dá)所探測(cè)空間的冰晶密度,即冰水含量(Ice Water Content,IWC)實(shí)現(xiàn),單位為g/m3[5]。冰水含量的準(zhǔn)確估計(jì)是檢測(cè)高空冰晶的前提。

機(jī)載氣象雷達(dá)是飛機(jī)探測(cè)危險(xiǎn)氣象保障飛機(jī)安全運(yùn)行的重要監(jiān)視設(shè)備。目前,實(shí)現(xiàn)冰晶冰水含量估計(jì)的方法主要是依據(jù)冰水含量與雷達(dá)反射率因子之間的經(jīng)驗(yàn)關(guān)系。1954年,Atlas通過研究飛機(jī)實(shí)測(cè)的譜參數(shù)數(shù)據(jù),得到雷達(dá)反射率因子-粒子分布-冰水含量之間的經(jīng)驗(yàn)關(guān)系[6]。1992年,Liao等人探討了Ka波段雷達(dá)和W波段雷達(dá)的冰晶冰水含量與雷達(dá)反射率因子之間的關(guān)系[7]。1995年,Brown用毫米波機(jī)載雷達(dá)探測(cè)了冰晶云,并給出了冰水含量與雷達(dá)反射率因子經(jīng)驗(yàn)公式[9]。2004年,Matrosov研究了利用雷達(dá)反射率因子以估測(cè)海洋性層云含水量的方法,比較了在不同的降水云與非降水云強(qiáng)度界限時(shí)云內(nèi)含水量的變化[10]。2015年,Protat給出了在35GHz和95GHz毫米波雷達(dá)下的冰水含量與雷達(dá)反射率因子統(tǒng)計(jì)關(guān)系[11]。上述研究分別針對(duì)不同場(chǎng)景下不同波段的雷達(dá)估計(jì)冰水含量時(shí)雷達(dá)反射率因子與冰水含量之間的關(guān)系進(jìn)行統(tǒng)計(jì)分析,得出了不同情況下雷達(dá)反射率因子與冰水含量之間的經(jīng)驗(yàn)關(guān)系。但由于每個(gè)研究人員擁有不同的數(shù)據(jù)集,所得經(jīng)驗(yàn)關(guān)系具有較大的局限性,且僅利用雷達(dá)反射率因子對(duì)冰水含量進(jìn)行估計(jì)的準(zhǔn)確性較差,難以依據(jù)高空冰晶濃度的分布進(jìn)行危險(xiǎn)性評(píng)判及區(qū)域告警,必須尋求更精準(zhǔn)的冰水含量估計(jì)方法。因此,開展高空冰晶冰水含量估計(jì)的研究具有非常重要的意義。2016年,美國(guó)霍尼韋爾國(guó)際公司的研究團(tuán)隊(duì)提出了一種機(jī)載X波段氣象雷達(dá)的高空冰晶檢測(cè)方法,并成功進(jìn)行飛行測(cè)試[12]。該算法基于機(jī)器學(xué)習(xí)的方法,利用冰晶的雷達(dá)反射率因子、高度和溫度等特征建立特征向量估計(jì)高空冰晶冰水含量,理論上能夠?qū)崿F(xiàn)冰晶檢測(cè)。雖然關(guān)于該功能的數(shù)據(jù)和算法原理以及相關(guān)雷達(dá)結(jié)構(gòu)等信息均未被公開,但利用可用信息實(shí)現(xiàn)對(duì)于冰晶冰水含量參數(shù)的預(yù)測(cè)、估計(jì)這種思想,為準(zhǔn)確實(shí)現(xiàn)冰水含量估計(jì)提供了研究思路。

回歸算法是常應(yīng)用于氣象預(yù)測(cè)、降水量估計(jì)等領(lǐng)域的一種以數(shù)據(jù)特征提取為基礎(chǔ),用數(shù)據(jù)分析的相關(guān)模型和數(shù)據(jù)挖掘技術(shù)進(jìn)行處理分析,用于估計(jì)或預(yù)測(cè)未知參數(shù)的數(shù)據(jù)處理技術(shù)[13]。這種特性為回歸算法應(yīng)用于高空冰晶冰水含量估計(jì)提供了理論基礎(chǔ)。目前,將回歸算法應(yīng)用于高空冰晶冰水含量估計(jì)的研究尚屬空缺。

本文提出了一種基于多重多元回歸的高空冰晶冰水含量估計(jì)方法,分析驗(yàn)證了將回歸算法引入高空冰晶冰水含量估計(jì)的可行性。該方法將冰晶冰水含量單一的目標(biāo)值估計(jì)問題轉(zhuǎn)化為多重多元回歸的冰水含量標(biāo)記向量預(yù)測(cè)問題,利用訓(xùn)練好的回歸模型最終實(shí)現(xiàn)高空冰晶的冰水含量估計(jì)。針對(duì)典型深對(duì)流區(qū)域仿真數(shù)據(jù)進(jìn)行實(shí)驗(yàn)分析表明,該方法可有效地估計(jì)高空冰晶的冰水含量,估計(jì)結(jié)果較為準(zhǔn)確。

1 基于多重多元回歸的高空冰晶冰水含量估計(jì)

基于多重多元回歸的高空冰晶冰水含量估計(jì)方法需要通過帶標(biāo)簽的訓(xùn)練數(shù)據(jù)建立回歸模型,建立多重多元回歸模型首先需要建立訓(xùn)練樣本的特征矩陣及標(biāo)簽矩陣,然后利用偏最小二乘回歸(Partial Least Squares Regression,PLSR)算法建立起冰晶樣本多維關(guān)鍵特征向量與該樣本冰水含量之間的映射關(guān)系。得到訓(xùn)練好的多重多元回歸模型后,將待測(cè)數(shù)據(jù)輸入模型得到冰晶冰水含量估計(jì)結(jié)果。下面分別對(duì)建立樣本特征矩陣和標(biāo)簽矩陣,建立多重多元回歸模型以及利用回歸模型完成高空冰晶冰水含量估計(jì)的過程進(jìn)行介紹。

1.1 建立特征矩陣和標(biāo)簽矩陣

建立回歸模型首先需要建立訓(xùn)練樣本的特征矩陣及標(biāo)簽矩陣。建立樣本特征矩陣時(shí)考慮高空冰晶的濃度分布受眾多因素的影響,除雷達(dá)反射率因子外,與冰晶結(jié)冰事件相關(guān)的氣象環(huán)境參數(shù)包括溫度、高度、氣壓、空氣密度和對(duì)流風(fēng)暴類型等也會(huì)限制冰晶的產(chǎn)生[5]。從冰晶的產(chǎn)生條件來看,利用氣象環(huán)境參數(shù)高度、溫度確定冰晶可能存在的區(qū)域是可行的。因此,本文依據(jù)影響高空冰晶冰水含量的關(guān)鍵特征參數(shù)建立特征矩陣。除雷達(dá)反射率因子外,再引入高度、溫度信息作為冰晶樣本的特征,建立的特征矩陣X為

(1)

其中,特征矩陣X的每一個(gè)行向量表示一個(gè)樣本的特征參數(shù)向量。假設(shè)有k個(gè)冰晶樣本數(shù)據(jù),第i個(gè)冰晶樣本xi(i=1,2,…,k);有3個(gè)特征參數(shù),表示為xi=(Zi,Hi,Ti),Z代表雷達(dá)反射率因子(單位:dBZ),H代表高度(單位:km),T代表溫度(單位:℃)。

(2)

(3)

1.2 建立回歸模型

在1.1節(jié)基礎(chǔ)上,本文采用多重多元回歸模型來求解由特征矩陣X到標(biāo)簽矩陣Y的映射問題,模型的求解目標(biāo)為利用訓(xùn)練數(shù)據(jù),得到冰晶樣本的特征向量與真實(shí)冰水含量值之間的最優(yōu)映射關(guān)系。將基于經(jīng)驗(yàn)公式的冰水含量值計(jì)算問題轉(zhuǎn)換為利用多自變量對(duì)因變量進(jìn)行預(yù)測(cè)的多重多元回歸分析問題,并采用偏最小二乘回歸方法對(duì)冰水含量進(jìn)行求解。偏最小二乘回歸模型結(jié)構(gòu)如圖1所示。

圖1 偏最小二乘回歸模型結(jié)構(gòu)圖

偏最小二乘回歸模型分別由外部模型和內(nèi)部模型組成,首先由外部模型對(duì)特征矩陣X和標(biāo)簽矩陣Y進(jìn)行成分提取,然后由內(nèi)部模型建立提取出的n對(duì)成分之間的關(guān)系,最后構(gòu)建特征矩陣X和標(biāo)簽矩陣Y的回歸模型。具體過程是首先對(duì)特征矩陣X和標(biāo)簽矩陣Y提取第一對(duì)主成分v1和u1,然后偏最小二乘回歸分別建立特征矩陣X與v1的回歸方程以及標(biāo)簽矩陣Y與v1的回歸方程。如果建立的回歸方程交叉有效性驗(yàn)證指標(biāo)小于門限值,則算法終止;否則,接下來利用特征矩陣X與v1建立回歸方程產(chǎn)生的殘差矩陣以及標(biāo)簽矩陣Y與u1建立回歸方程產(chǎn)生的殘差矩陣進(jìn)行第2輪的成分提取。如此反復(fù)迭代,直到回歸方程的交叉有效性指標(biāo)達(dá)到門限值時(shí)停止迭代。若最終共提取了n個(gè)成分,則利用上述成分v1,v2,…,vn構(gòu)建的X與利用上述成分v1,v2,…,vn構(gòu)建的Y建立標(biāo)簽矩陣Y關(guān)于特征矩陣X的回歸方程。以下對(duì)這一過程的原理及實(shí)現(xiàn)進(jìn)行詳細(xì)闡述。

訓(xùn)練樣本數(shù)據(jù)含有多維自變量和多維因變量,偏最小二乘算法首先提取訓(xùn)練樣本數(shù)據(jù)的特征矩陣X和標(biāo)簽矩陣Y中的第一對(duì)主成分v1和u1,其中,主成分v1和u1可以作為原始特征矩陣X和標(biāo)簽矩陣Y的線性變換,權(quán)重系數(shù)向量分別設(shè)為p1和q1,即v1=Xp1,u1=Yq1。為了滿足回歸分析的需要,應(yīng)使v1和u1盡可能多地?cái)y帶樣本特征矩陣X和標(biāo)簽矩陣Y的信息,即v1和u1方差最大,且滿足v1和u1相關(guān)程度達(dá)到最大,即相關(guān)性系數(shù)達(dá)到最大值,綜上當(dāng)v1和u1的協(xié)方差達(dá)到最大即可滿足要求[13]。為了得到成分v1和u1首先需要求得p1和q1,即求解滿足式(4)目標(biāo)函數(shù)的p1和q1:

(4)

引入拉格朗日乘子算法可得:

(5)

分別對(duì)p1和q1求偏導(dǎo),并令之為0,得:

(6)

(7)

由式(7)可知,p1是XTYYTX的最大特征值對(duì)應(yīng)的單位特征向量,q1是YTXXTY最大特征值對(duì)應(yīng)的單位特征向量。求得p1和q1后,即可得到成分:

(8)

求解出第一對(duì)成分v1和u1后,建立回歸方程:

(9)

式(9)中,c1=XTv1/‖v1‖2;r1=YTv1/‖v1‖2是回歸系數(shù)向量即投影向量;E1,F1是回歸方程的殘差矩陣。

(10)

(11)

(12)

(13)

(14)

繼續(xù)建立回歸方程為

(15)

將式(15)第二個(gè)公式帶入式(9)第二個(gè)公式可得:

(16)

以此類推,不斷對(duì)殘差矩陣建立回歸方程,對(duì)上述過程進(jìn)行迭代,迭代過程中通過交叉有效性驗(yàn)證判斷是否停止迭代,若最終對(duì)X和Y共提取了n對(duì)成分,最終可得到的回歸模型為

(17)

迭代計(jì)算過程中可得:

(18)

由于式(17)中v1,v2,…,vn中任一成分vh(h=1,2,…,n)是原始特征矩陣X的線性組合[15],即

(19)

其中:

(20)

式(20)中,I為單位向量,p和c均為迭代過程中計(jì)算得到。因此可將公式(17)可轉(zhuǎn)化為標(biāo)簽矩陣Y關(guān)于特征矩陣X的回歸方程:

(21)

記為

(22)

D是偏最小二乘回歸方程的回歸系數(shù)向量,則有

Y=XD+Fn

(23)

其中,Fn為最終達(dá)到迭代要求時(shí)的殘差矩陣,在后續(xù)使用測(cè)試數(shù)據(jù)進(jìn)行冰水含量估計(jì)時(shí)可忽略[15]。

至此,針對(duì)冰晶冰水含量估計(jì)的多重多元回歸模型訓(xùn)練完成。接下來利用訓(xùn)練好的回歸模型對(duì)待測(cè)氣象數(shù)據(jù)樣本進(jìn)行冰晶冰水含量估計(jì)。

1.3 高空冰晶冰水含量估計(jì)

(24)

(25)

冰晶冰水含量估計(jì)結(jié)果大于0,則說明存在冰晶,冰水含量越大說明該空域內(nèi)冰晶含量越高,根據(jù)冰水含量的估計(jì)結(jié)果可判斷冰晶是否存在,且冰水含量估計(jì)結(jié)果可作為后續(xù)冰晶的危險(xiǎn)程度判決依據(jù)。

2 方法流程與實(shí)驗(yàn)步驟

實(shí)驗(yàn)數(shù)據(jù)是利用WRF(Weather Research and Forecasting,WRF)天氣模擬與預(yù)報(bào)軟件模擬得到的氣象場(chǎng)景數(shù)據(jù)以及利用氣象雷達(dá)回波仿真得到的雷達(dá)回波數(shù)據(jù)。仿真數(shù)據(jù)包含建立特征向量需要的雷達(dá)反射率因子,計(jì)算溫度、高度信息所需要的氣壓、海拔高度、擾動(dòng)位溫等數(shù)據(jù)和計(jì)算冰水含量真值所需要的冰水混合比、空氣密度等數(shù)據(jù)。多重多元回歸的高空冰晶冰水含量估計(jì)方法流程圖如圖2所示。冰晶冰水含量估計(jì)的計(jì)算步驟如下:

圖2 算法流程

1)步驟一:對(duì)冰晶樣本數(shù)據(jù)進(jìn)行質(zhì)量控制;

2)步驟二:生成樣本特征矩陣和標(biāo)簽矩陣;

3)步驟三:將所有冰晶樣本數(shù)據(jù)劃分為訓(xùn)練數(shù)據(jù)集和測(cè)試數(shù)據(jù)集;

4)步驟四:將訓(xùn)練數(shù)據(jù)集的特征矩陣和標(biāo)簽矩陣帶入回歸模型,利用偏最小二乘回歸算法建立回歸模型;

5)步驟五:將標(biāo)準(zhǔn)化的回歸變量還原成原始變量,確定回歸模型的最終形式;

6)步驟六:將測(cè)試數(shù)據(jù)集的特征矩陣輸入確定好的回歸模型中,計(jì)算測(cè)試樣本的冰水含量值;

7)步驟七:對(duì)冰晶樣本冰水含量估計(jì)結(jié)果與樣本真實(shí)冰水含量進(jìn)行對(duì)比統(tǒng)計(jì)分析,判斷估計(jì)結(jié)果的準(zhǔn)確性,驗(yàn)證算法性能。

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

本文針對(duì)墨西哥灣地區(qū)緯度北緯28.8°～30.4°,經(jīng)度西經(jīng)85.9°～87.7°區(qū)域內(nèi)的含有高空冰晶天氣的氣象仿真數(shù)據(jù)進(jìn)行實(shí)驗(yàn)分析。該區(qū)域位于沿海地區(qū)屬于熱帶氣候,仿真氣象場(chǎng)景時(shí)間選取2015年8月16日,當(dāng)天該地區(qū)有典型深對(duì)流云產(chǎn)生。高空冰晶主要出現(xiàn)在對(duì)流云頂層的冰云中,其高度主要分布在6～15km范圍內(nèi)。為了驗(yàn)證本文方法的有效性,分別對(duì)不同高度層的氣象場(chǎng)景數(shù)據(jù)進(jìn)行冰晶冰水含量估計(jì)實(shí)驗(yàn)。如圖3所示為針對(duì)該場(chǎng)景不同高度層的氣象數(shù)據(jù)分別采用本文方法與Atlas、Aydin和Mace等人[10]分別給出的經(jīng)驗(yàn)公式法進(jìn)行冰晶冰水含量估計(jì)的結(jié)果對(duì)比圖。

圖3 不同高度層數(shù)據(jù)冰晶冰水含量估計(jì)結(jié)果對(duì)比圖

由圖3所示對(duì)比結(jié)果可以看出,傳統(tǒng)經(jīng)驗(yàn)公式法僅使用雷達(dá)反射率因子對(duì)冰水含量進(jìn)行估計(jì)且受經(jīng)驗(yàn)取值限制,泛用性差,冰晶冰水含量估計(jì)結(jié)果均與真實(shí)值之間存在較大的偏差。而采用本文所述多重多元回回歸算法進(jìn)行冰晶冰水含量估計(jì),冰晶冰水含量估計(jì)值與真實(shí)值基本一致,誤差較小,說明采用多重多元回歸算法實(shí)現(xiàn)的冰晶冰水含量估計(jì)效果優(yōu)于經(jīng)驗(yàn)公式僅利用雷達(dá)反射率因子進(jìn)行冰晶冰水含量估計(jì)的方法。

分別對(duì)不同高度層數(shù)據(jù)進(jìn)行實(shí)驗(yàn),對(duì)冰水含量估計(jì)結(jié)果統(tǒng)計(jì)分析如下,冰晶冰水含量估計(jì)結(jié)果與真值之間相關(guān)系數(shù)如表1所示。

表1 冰晶冰水含量估值與真值之間的相關(guān)系數(shù)

相關(guān)性系數(shù)越接近1說明算法估計(jì)的越準(zhǔn)確,由表1可見,不同高度層的測(cè)試數(shù)據(jù)冰水含量估值與真值之間的相關(guān)系數(shù)都大于0.8,平均值達(dá)到了0.8884,具有很強(qiáng)的相關(guān)性,說明誤差較小估計(jì)結(jié)果較為準(zhǔn)確。

圖4為本文方法得到的高空冰晶冰水含量與真實(shí)值的相關(guān)性散點(diǎn)圖。陰影數(shù)據(jù)點(diǎn)橫坐標(biāo)所對(duì)應(yīng)的是采用本文方法得的冰水含量值,縱坐標(biāo)是冰水含量真值,當(dāng)冰晶冰水含量的估計(jì)值完全等于其真實(shí)值時(shí),每個(gè)樣本點(diǎn)的橫縱坐標(biāo)數(shù)值均相等,即分布在對(duì)角線上,則越接近對(duì)角線表示冰水含量估計(jì)值的準(zhǔn)確程度越高。由圖4可見,冰晶樣本點(diǎn)大部分都分布在對(duì)角線附近,說明算法估計(jì)結(jié)果都較為準(zhǔn)確。

圖4 不同高度層數(shù)據(jù)冰晶冰水含量估計(jì)結(jié)果相關(guān)性散點(diǎn)圖

對(duì)于不同高度層測(cè)試數(shù)據(jù)的冰水含量回歸估計(jì)結(jié)果進(jìn)行誤差統(tǒng)計(jì),統(tǒng)計(jì)誤差小于0.05g/m3的樣本數(shù)占總樣本數(shù)的百分比,結(jié)果如表2所示。

表2 誤差小于0.05g/m3的樣本數(shù)占比

由表2可見,對(duì)不同高度層的數(shù)據(jù)進(jìn)行實(shí)驗(yàn),計(jì)算得到誤差小于0.05g/m3的樣本數(shù)占總樣本數(shù)的比例平均值為75.67%,說明大部分樣本的冰晶冰水含量估計(jì)結(jié)果較為準(zhǔn)確。以9.8 km高度層數(shù)據(jù)為例對(duì)冰水含量估計(jì)結(jié)果進(jìn)行統(tǒng)計(jì)分析,表3為冰晶樣本原始數(shù)據(jù)和回歸算法估計(jì)的結(jié)果中不同冰水含量的冰晶樣本占總體數(shù)量的比重,結(jié)合表3結(jié)果分析可知:回歸算法估計(jì)結(jié)果與真實(shí)情況基本一致,大部分樣本都小于1g/m3,大于1.5g/m3的樣本估計(jì)結(jié)果相差較大是因?yàn)榇笥?.5g/m3的樣本數(shù)量較少。

表3 不同濃度范圍的冰晶樣本占數(shù)據(jù)總體的比重

4 結(jié)束語

本文提出了一種基于多重多元回歸的高空冰晶冰水含量估計(jì)方法,分析驗(yàn)證了將回歸算法引入高空冰晶冰水含量估計(jì)的可行性。將冰晶的雷達(dá)反射率因子、高度和溫度等特征作為代表冰晶樣本的特征向量,將冰晶目標(biāo)的冰水含量估計(jì)問題轉(zhuǎn)化為多重多元回歸預(yù)測(cè)問題,利用訓(xùn)練好的回歸模型最終實(shí)現(xiàn)高空冰晶的冰水含量估計(jì)。針對(duì)典型深對(duì)流區(qū)域仿真數(shù)據(jù)進(jìn)行實(shí)驗(yàn)分析表明,該方法可有效地估計(jì)高空冰晶的冰水含量,相比于僅利用雷達(dá)反射率因子的經(jīng)驗(yàn)公式方法,結(jié)果更為準(zhǔn)確。

附錄：

(A-1)

將公式(A-1)第一個(gè)式子轉(zhuǎn)置,第二個(gè)式子保持不變可得

(A-2)

(A-3)

即可得到λ1=θ1。

則公式(6)可轉(zhuǎn)換為

(A-4)