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

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

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

0 引言

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

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

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

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

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

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

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

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

(1)

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

(2)

(3)

1.2 建立回歸模型

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

圖1 偏最小二乘回歸模型結構圖

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

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

(4)

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

(5)

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

(6)

(7)

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

(8)

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

(9)

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

(10)

(11)

(12)

(13)

(14)

繼續建立回歸方程為

(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)可轉化為標簽矩陣Y關于特征矩陣X的回歸方程:

(21)

記為

(22)

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

Y=XD+Fn

(23)

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

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

1.3 高空冰晶冰水含量估計

(24)

(25)

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

2 方法流程與實驗步驟

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

圖2 算法流程

1)步驟一:對冰晶樣本數據進行質量控制;

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

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

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

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

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

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

3 實驗結果分析

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

圖3 不同高度層數據冰晶冰水含量估計結果對比圖

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

分別對不同高度層數據進行實驗,對冰水含量估計結果統計分析如下,冰晶冰水含量估計結果與真值之間相關系數如表1所示。

表1 冰晶冰水含量估值與真值之間的相關系數

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

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

圖4 不同高度層數據冰晶冰水含量估計結果相關性散點圖

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

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

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

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

4 結束語

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

附錄：

(A-1)

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

(A-2)

(A-3)

即可得到λ1=θ1。

則公式(6)可轉換為

(A-4)