基于多重多元回歸的高空冰晶冰水含量估計方法

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

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

0 引言

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

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

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

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

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

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

1.1 建立特征矩陣和標簽矩陣

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

(1)

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

(2)

(3)

1.2 建立回歸模型

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

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

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

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

(4)

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

(5)

分別對p1和q1求偏導,并令之為0,得:

(6)

(7)

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

(8)

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

(9)

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

(10)

(11)

(12)

(13)

(14)

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

(15)

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

(16)

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

(17)

迭代計算過程中可得:

(18)

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

(19)

其中:

(20)

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

(21)

記為

(22)

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

Y=XD+Fn

(23)

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

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

1.3 高空冰晶冰水含量估計

(24)

(25)

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

2 方法流程與實驗步驟

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

圖2 算法流程

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

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

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

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

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

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

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

3 實驗結(jié)果分析

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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié)束語

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

附錄：

(A-1)

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

(A-2)

(A-3)

即可得到λ1=θ1。

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

(A-4)