閆文輝，黃興友，趙鈺錦，張瑩，楊濤

(1. 中國飛行試驗研究院氣象臺，陜西 西安 710089；2. 南京信息工程大學氣象災害預報預警與評估協同創(chuàng)新中心，江蘇 南京 210044)

1 引 言

陣風載荷指飛行器穿越陣風時，升力面產生的載荷增量，是大氣擾動的客觀表現[1]。氣象學上把陣風載荷稱為飛機顛簸，顛簸特別嚴重時會破壞飛機結構，形成安全隱患，甚至造成飛行事故。2017年 6月 18日，東航空客 330-200E 客機在俄羅斯境內遭遇極重度顛簸，導致28 人受傷。Sharman 等[2]估計美國民航飛機每年遭遇6 300 次中到重度顛簸，5 000 次重度到極重度顛簸。Williams等[3]估計每年飛機顛簸造成全美約2億美元的經濟損失。2019年中國飛行試驗研究院試飛任務超5 000架次，試飛任務的不斷增長致使飛行器遭遇飛機顛簸的風險也不斷增加，嚴重影響試飛安全與效率。

飛機顛簸與大氣湍流有著密切的聯系，Kelvin-Helmholtz (K-H)不穩(wěn)定是最為公認的大氣湍流產生機制[4]。大氣湍流區(qū)域的中值長度約60 km，中值厚度約1 km[5]。大氣湍流區(qū)域存在著尺度不同的湍渦，飛機顛簸是由那些與飛機尺度相當的、混亂的湍渦造成的，這種湍渦的直徑一般在15～150 m[6]。Meneguz 等[7]研究表明只有約 14%飛機顛簸是發(fā)生在對流天氣附近，這意味著造成飛機顛簸的主要原因是切變重力波破碎和地形因素。Jaeger等[8]發(fā)現約66.7%的飛機顛簸發(fā)生在高空急流附近，因此高空急流是誘發(fā)飛機顛簸最主要的原因。

俞飛等[9]發(fā)現高空急流、高空槽、切變線以及高空脊是造成華北地區(qū)飛機顛簸的主要天氣形勢。這種基于天氣形勢分析的方法主觀經驗性太強，易造成很高的漏報率和空報率。因此，運用數值模式輸出量計算飛機顛簸指數進行顛簸診斷和預報目前已經成為主流。長期研究建立了大量飛機顛簸單一指數，主要有Brown 指數、Dutton 指數、L-P指數、Ri 數以及 Ellrod 指數等[10-19]。徐佳男等[20]選取診斷效果較好的3種指數，統(tǒng)計了中國地區(qū)飛機顛簸的氣候特征。但不同指數重點考慮的因素有所差異，因此在不同區(qū)域、不同季節(jié)及不同成因的飛機顛簸中，指數自身及指數之間的預報能力有所不同。Sharman 等[2]提出集成指數預報算法（Graphical Turbulence Guidance, GTG），將幾種顛簸指數組合使用，其總體預報性能比使用單一指數診斷效果要好。Kim 等[21]將GTG 方法應用于東亞晴空顛簸的研究，并估計了當前晴空顛簸數值預報水平。沈強等[22]選取了9個預報效果較好的顛簸指數，根據預報得分計算權重，構建了飛機顛簸綜合指數。

為拓展民機遭遇飛機顛簸時的飛機包線、提升飛機品質、完善民航試航取證大綱，中國飛行試驗研究院以ARJ21 支線客機作為試驗機，承接民機陣風載荷測量試飛型號研發(fā)項目。陣風載荷測量試飛在國內是首次進行的，具有技術新、找風難、風險高及工作量大等特點。為準確預報潛在飛機顛簸區(qū)域，解決“找風難”問題，同時避免其他航空危險天氣對試飛安全的影響，本文選取對流層高層（7 500 m 以上）誘發(fā)飛機顛簸的常見天氣系統(tǒng)（高空急流），利用歐洲中期天氣預報中心（European Centre for Medium - Range Weather Forecast, ECMWF）提供的 ERA5 再分析資料，從飛機顛簸飛行員天氣報告（Pilot Reports, PIREPs）入手，基于10個應用成熟的飛機顛簸指數，通過加權集成給出了顛簸潛勢綜合指數，同時對顛簸潛勢綜合指數的性能進行了評估。

2 資料介紹

2.1 ERA再分析資料

ERA5 再分析資料是 ECMWF 利用 4DVar 同化產生的第5代大氣再分析全球氣候數據，該數據提供1979年至今的逐小時大氣、陸地和海洋的氣候變量估計值，其中大氣數據的空間分辯率精確到了0.25 °×0.25 °網格，等壓面數據垂直方向包括1～1 000 hPa的 37 層等壓面[23]。數據包主要包括位勢高度、溫度、經向風、緯向風、垂直風、相對濕度、比濕、散度及渦度等參量的網格數據。

大氣再分析數據在飛機顛簸指數的診斷和預報中具體廣泛的應用，為了進一步說明ERA5再分析資料的準確性，本文隨機下載了2020年9月21日08 時（北京時間，下同）300 hPa 天氣圖(圖1)，圖1a 為 ERA5 再分析資料，圖 1b 為 ECMWF 細網格預報資料中08 時的天氣實況。對比可見，二者基本一致，歐亞中高緯度均成兩槽一脊型，中緯度為緯向平直氣流，中低緯度受反氣旋控制，并且風場、等溫線及等高線的分布、大小基本一致。

圖1 2020年9月21日08時300 hPa天氣圖

2.2 飛行員天氣報告（PIREPs）

在飛行過程中若遇到飛機顛簸、飛機積冰等現象時，機組人員會以話音通信方式將現象、時間、位置、高度及強度等信息傳輸給管制單位，管制人員再將該信息發(fā)送給氣象人員，氣象人員將其記錄并通過專業(yè)技術分析增添造成原因、備注等信息，最后進行統(tǒng)計保存，這種資料稱為飛行員天氣報告（PIREPs）[24]。為了準確評估各單一顛簸指數的預報性能，并推算顛簸潛勢綜合指數，本文篩選了 2017年 10月—2020年 12月期間，中國范圍內7 500 m 高度以上高空急流背景下的170 例飛機顛簸報告及170 例無顛簸樣本，其中280 例（140 例發(fā)生顛簸、140 無顛簸）樣本用于構建顛簸潛勢綜合指數，剩余的60 例（30 例發(fā)生顛簸、30 無顛簸）樣本用于評估顛簸潛勢綜合指數。具體報告位置分布如圖2所示，顛簸產本基本分布于我國中緯度地區(qū)，與溫帶急流和副熱帶急流位置對應較為良好。需要注意的是，無顛簸樣本選取閻良機場附近9 000 m左右的高度上一個固定位置，試飛期間當飛機未發(fā)生顛簸時取樣一次，所以只有一個點。

圖2 飛機顛簸報告位置分布圖

3 單一顛簸指數介紹

飛機顛簸和大氣湍流活動強弱有密切的關系，兩者皆發(fā)生在風場垂直切變區(qū)、水平切變區(qū)、輻合或輻散區(qū)、水平形變區(qū)以及強的水平溫度梯度區(qū)[10]。由于不同指數側重因素不同，其預報效果也不同，綜合考慮這五個方面，參考前人研究[2,21,26]，選取反映上述特性的，預報效果較為理想的10個飛機顛簸單一指數來構建綜合指數。

(1)L-P指數。

前蘇聯水文氣象中心綜合考慮了溫度切變、水平風場的水平切變和垂直切變情況，構建了L-P指數來預報潛在飛機顛簸區(qū)[10]，L-P指數采用如下計算方式：

式中，SV為水平風場的垂直切變(單位：m/(s·100 m))；ST為水平溫度梯度(單位：K/(100 km))；SH為水平風場的水平切變(單位：m/(s·100 km))；u和v為水平風分量(單位：m/s)；T為氣溫（單位：K）；z為垂直高度(單位：位勢米)。當L＜0 時，預報沒有顛簸，反之在此基礎上根據L指數計算概率P：

P的大小反映了飛機顛簸強度。

(2) -Ri數。

Richardson 數（Ri）被定義為表征大氣穩(wěn)定度的 Brunt-V?is?l? 頻率的平方和表征破壞大氣穩(wěn)定度的垂直風切變的平方之比[11]，其表達式為：

式中，θv為虛位溫(單位：K)；SV為水平風場的垂直切變(單位：s-1)；g為重力加速度(單位：m/s2)，由于Ri數越小，顛簸強度越大，為了方便計算，取-Ri數作為本文顛簸指數。對于低靜力穩(wěn)定度和高垂直風切變區(qū)域，例如高空急流或鋒區(qū)，如果晴空湍流將發(fā)展，通常-Ri大于0.25。

(3) Ellrod指數。

Ellrod 等[12]對鋒生函數進行簡化，綜合考慮了水平風場的垂直切變、總形變以及輻合輻散，提出了兩種飛機顛簸指數TI1和TI2，其表達式如下：

式中，SV為水平風場的垂直切變(單位：s-1)；DEF為水平風場的總形變(單位：s-1)，包括拉伸形變和切變形變兩項；Div為水平風場的散度(單位：s-1)。在業(yè)務應用過程中，王洪芳等[25]發(fā)現散度項的指示性不強，所以本文采用TI1指數，TI1高值區(qū)對應強的飛機顛簸區(qū)。

(4) MOSCAT指數。

MOSCAT 指數是由NECP 嵌套網格模式NGM（Nested Grid Model）基 于 MOS（Model Output Statistics）方法統(tǒng)計輸出的晴空湍流指數，又稱概率預報因子指數[13]，其表達式為：

式中，|V|為水平合成風速(單位：m/s)，DEF為水平風場的總形變(單位：s-1)。

(5) Brown指數。

Brown[14]利用熱成風關系簡化了Ri 傾向方程，提出了Brown指數：

式中，ζa=ζ+f為絕對渦度(單位：s-1)，ζ為相對渦度，f=2Ωsinφ為科氏參數。DSH為切變形變項(單位：s-1)，DST為拉伸形變項(單位：s-1)，SV為水平風場的垂直切變(單位：s-1)。由于大氣運動為準水平，所以將渦度簡化為垂直方向的分量。

(6) Dutton指數。

Dutton 指數[15]綜合考慮水平風場的水平切變和垂直切變對飛機顛簸強度的貢獻，得到一個經驗指數，其表達式為：

式中，SH為水平風場的水平切變(單位：m/(s·100 km))；SV為水平風場的垂直切變(單位：m/(s·km))；10.5為經驗常數。

(7) 曲率方法指數。

綜合分析飛機顛簸報告和環(huán)流形勢，發(fā)現高空槽脊的曲率與飛機顛簸強度相關，而曲率與渦度有關[16]，定義曲率指數為：

式中，ζ2單位為 s-2，u和v為水平風分量(單位：m/s)，此處渦度也取垂直方向的分量。

(8) 水平溫度梯度指數。

水平溫度梯度指數是一個從熱成風關系中計算形變和垂直風切變的方法[17]，其表達式為：

式中，HTG單位為K/m，T為氣溫 (單位：K)。

(9) 水平散度指數。

水平風場的強輻合輻散區(qū)與風場的切變區(qū)域對應良好，并且輻合可以對鋒生產生作用[2]。定義水平散度指數為：

式中，DIV單位為s-1，u和v為水平風分量。

(10) 風相關指數。

急流背景下的飛機顛簸與大的風速相關[11]，定義風相關指數如下：

式中，S 單位為 m/s，u和v為水平風分量(單位：m/s)。

4 顛簸潛勢綜合指數算法

由于不同指數重點考慮的因素有所差異，因而指數自身及指數之間的預報能力也有所差異，為確保顛簸指數能夠有效識別多種復雜條件下的飛機顛簸，本文參考美國國家大氣研究中心（National Center for Atmospheric Research,NCAR）Sharman 等[2]提供的 GTG 算法，給出了我國高空急流背景下的顛簸潛勢綜合指數。

4.1 數據歸一化

由于每個顛簸指數的單位和數值不同，計算結果的量級也相差較大，難以進行比較和加權計算，因此需要進行歸一化處理。本文根據各單一顛簸指數的閾值，基于分段線性函數將計算結果進行歸一化處理，將飛機顛簸強度分為無、輕度、中度、重度和極重度5個量級，中值分別對應數值0.0、0.25、0.5、0.75 和1.0。單一顛簸指數的具體閾值選取參考了 Sharman 等[2]、Kim 等[21]及 Willims等[26]的工作，其分類方法和原理為根據大量飛機顛簸報告強度等級與指數計算值之間的對應關系，統(tǒng)計各顛簸指數計算值的中值為該指數在當前量級的強度閾值，具體閾值如表1所示。

表1 高空急流背景下10種顛簸指數強度閾值劃分

4.2 確定顛簸指數權重

在完成各顛簸指數歸一化處理之后，為了定量評估各指數的預報水準，引入基于檢測概率的評估方法POD(Probability of Detection)[27]，包括發(fā)生顛簸的預報準確率PODY、未發(fā)生顛簸的預報準確率PODN 及總體預報準確率PODA，其表達式如下：

式中，YY、YN 分別表示發(fā)生了顛簸，對應的指數預報正確和錯誤，NN、NY 分別表示未發(fā)生顛簸，對應的指數預報正確和錯誤。PODY 值越大，漏報率越低；PODN值越大，空報率越低；PODA值越大，總體預報效果越好。

為了更加精準地評估各指數的預報水準，再引入了 TSS（True Skill Score）評分[28]來評估總體預報準確性，其表達式如下：

TSS 取值越大，表明該指數的診斷效果越好。根據預報準確率計算各顛簸指數預報得分φ[2]，計算公式如下：

式中，C 和p 為常數，用于調節(jié)fMOG的重要性，通常取C = 1，p = 0.25。fMOG為在數值計算過程中得到的顛簸強度為中度及以上的格點數與總網格點數的比值。

根據得到的預報得分，可計算出每個顛簸指數的權重Wn[2]，具體表達式如下：

4.3 綜合顛簸指數

CAT(i, j, k)的取值范圍為[0, 1]，閾值為 0、0.25、0.5、0.75 及1，分別對應顛簸強度為無、輕度、中度、重度和極重度5個量級的中值。

5 顛簸潛勢綜合指數性能評估

選取與飛機報告時間、高度最臨近的ERA5再分析資料作為輸入量，計算各顛簸報告點的單一指數歸一化值，然后利用PODY、PODN、PODA、TSS 及fMOG對各單一指數進行綜合評估，并計算預報得分φ，最終結果如表2 所示?？梢?，在單一指數中，L-P指數發(fā)生顛簸時的預報準確率PODY 最高，達 97.9%，漏報率最低，TI 指數和 Dutton 指數次之，PODY 分別達到了67.9%和67.1%。Dutton指數未發(fā)生顛簸時的預報準確率PODN 最高，達82.1%，空報率最低，Brown 指數和TI1 指數次之，PODN 分別達到了77.1%和69.3%。從總體預報準確率PODA 和TSS診斷評分來看，L-P指數因具有最高的PODA(76.1%)和TSS(0.521)，總體預報準確率最好，Dutton 指數及Brown 指數的總體預報效果次之，PODA分別為74.6%及71.8%，TSS評分分別為0.493 及0.436。同時，從中度及以上顛簸所占面積fMOG來看，Dutton 指數、Brown 指數及TI1 指數的fMOG分別為7.8%、7.9%及12.4%，接近徐佳男等[20]利用AMADR 資料統(tǒng)計出的中度及以上顛簸僅有 10%左右的結果，較為合理。L-P指數識別到的中度及以上顛簸所占面積fMOG為20.4%，遠大于10%，識別面積過大。綜合考慮預報準確率和中度以上顛簸所占面積，Dutton 指數具有最高的φ(1.043)，總體預報效果最好，Brown指數、L-P指數及TI1 指數次之。另外，水平散度指數 DIV、曲率方法指數ζ2及-Ri 指數是這 10個單一指數中預報效果最差的3個。

表2 280例建模樣本的單一顛簸指數性能評估表

在高空急流背景下，Dutton指數、Brown指數、L-P指數及TI1 指數預報顛簸的效果好，分析物理機理可見，水平風場的水平切變、垂直切變以及水平形變的貢獻最大，強的水平溫度梯度和較大的風速也有一定的影響。而水平散度指數DIV、曲率方法指數ζ2及-Ri指數診斷錯誤的原因可能是，對于對流層高層大氣，大氣環(huán)流以緯向型為主，輻合輻散、氣旋反氣旋特征不明顯，且大氣穩(wěn)定度較高。

然后根據第3 節(jié)中介紹的流程對各單一指數進行加權集成，形成最終的顛簸潛勢綜合指數，利用剩余的60例顛簸樣本對顛簸潛勢綜合指數進行驗證。為了更加深入地評價顛簸潛勢綜合指數的優(yōu)越性，本文額外引入了模型性能評價體系中廣泛 使 用 的 ROC 曲 線（Receiver Operating Characteristic Curve）及 AUC 值（Area Under Curve）[29]。ROC 曲線，它是以 1-PODN 為橫軸，PODY 為縱軸繪制的曲線，它越靠近左上角，表明指數模型區(qū)分發(fā)生顛簸和未發(fā)生顛簸的性能越好。AUC 就是ROC 曲線和x軸（1-PODN 軸）之間的面積，反映了模型把發(fā)生顛簸的樣本排在未發(fā)生顛簸樣本前面的比例（如果AUC=1，說明模型100%的將所有發(fā)生顛簸的樣本排在未發(fā)生顛簸樣本前面）。

60 例驗證樣本的顛簸指數性能評價結果和ROC 曲線如表 3 及圖 3 所示，圖 3a 為 10個單一指數的 ROC 曲線，圖 3b 為 9個綜合指數的 ROC 曲線，其中CAT2 到CAT10 分別表示根據預報得分φ的高低，由2個到10個單一指數組成的綜合指數?？梢?，10個單一指數的預報效果均好于隨機預報模型，且綜合指數較好地集成了各個單一指數的優(yōu)點，預報水平始終優(yōu)于各單一指數。CAT5 為這9個顛簸潛勢綜合指數中最好的綜合指數，由預報效果最好的五個單一指數（Dutton 指數、Brown 指數、L-P指數、TI1 指數及風相關指數S）加權集成，其 權 重Wn分 別 為 0.237、0.220、0.205、0.186 及0.152，預報得分φ達 1.225，且 PODY 達 86.7%，PODN 達 93.3%，PODA 達 90%，TSS 達 0.80，中度及以上顛簸所占面積fMOG為9.2%，各項指標均較為優(yōu)秀，總體效果最好。同時，發(fā)現60例驗證樣本統(tǒng)計得到的單一指數性能與280 例建模樣本統(tǒng)計的單一指數性能略有差異，說明樣本量的大小和樣本間的差異性會影響綜合指數建模的好壞，但總體上Dutton指數、Brown指數及L-P指數效果最好，水平散度指數DIV及曲率方法指數ζ2效果最差。

圖3 60例驗證樣本的ROC曲線

表3 60例驗證樣本的單一顛簸指數及最優(yōu)綜合指數性能評估表

根據預報得分φ的大小逐步增加單一指數數量，驗證顛簸潛勢綜合指數對單一指數數量的敏感性，結果如圖4a所示。可見，當單一指數增加到5個的時候，綜合指數效果最優(yōu)，預報得分φ達1.225，之后隨著單一指數數量的增加，對綜合指數預報效果開始產生負作用。為了進一步驗證顛簸潛勢綜合指數對參與建模樣本數量的敏感性，將參與建模的樣本從40個逐步增加到280個，計算得到的顛簸潛勢綜合指數最優(yōu)組合的預報得分φ如圖4b 所示?？梢?，當建模樣本數量為40個時，綜合指數的預報得分φ最大，但隨機誤差較大，綜合指數預報準確性置信度較低。綜合指數的預報得分φ總體波動性不強，但是也隨樣本量的增加，綜合指數性能逐步提升。

圖4 顛簸潛勢綜合指數預報得分φ對單一指數數量敏感性圖(a)及對參與建模樣本量敏感性圖(b)

本文選取一例樣本來展示諸單一指數和綜合指數的預報差異。2020年1月6日12:52，在蒲城及以東地區(qū)8 500 m 高度上，試飛院參試飛機（C919）出現中到重度顛簸，持續(xù)20 分鐘以上。選取與顛簸報告出現時間、高度最臨近的ERA5再分析資料（2020年 1月 6日 13:00，300 hPa、350 hPa及400 hPa三層等壓面）作為數據源，350 hPa高空天氣圖如圖5a 所示，陜西中南部及以東地區(qū)出現高空急流，急流軸風速大于40 m/s，顛簸報告點位于急流軸北側風速切變最明顯和溫度梯度最大的區(qū)域。10個單一顛簸指數和顛簸潛勢綜合指數預報結果如圖5b～5l所示，圖中填色表示顛簸的區(qū)域和強度?？梢姡孙L相關指數S 外，其余指數位置都診斷正確，飛機顛簸報告點落在預報區(qū)域內。同時，各指數診斷的顛簸區(qū)域面積、強度分布均有較大差異。

圖5 2020年1月6日12:52高空急流背景下天氣形勢圖及各顛簸指數診斷結果

6 結論與討論

本文利用飛機顛簸報告資料和ERA5 再分析資料，基于10個應用成熟的單一顛簸指數，通過加權集成構建了我國高空急流背景下的顛簸潛勢綜合指數。

(1) Dutton 指數因具有最高的預報得分φ(1.043)，是這10個單一指數中總體預報效果最好的指數，Brown指數、L-P指數及TI1指數次之。另外，水平散度指數DIV、曲率方法指數ζ2及-Ri 指數是這10個單一指數中預報效果最差的3個指數。

(2) 隨著單一指數數量的增加，綜合指數的預報效果先增大后減小，當單一指數增加到5個的時候，即由Dutton指數、Brown指數、L-P指數、TI1指數及風相關指數S 加權集成的時候，綜合指數總體效果最優(yōu)，PODY 達 86.7%，PODN 達 93.3%，PODA 達 90%，TSS 達 0.80，中度及以上顛簸所占面積fMOG為9.2%，預報得分φ達1.225，漏報率和空報率都得到了很好的改善，各項指標均較為優(yōu)秀，總體效果最好。

(3) 隨著參與建模的飛機顛簸樣本量增加，綜合指數性能逐步提升，隨機誤差逐步減小，置信度逐步提升。

該顛簸潛勢綜合指數可用于民機陣風載荷測量試飛解決“找風難”問題。飛機顛簸潛勢綜合指數的準確性取決于ERA5再分析資料的精度，同時該綜合指數的評估效果取決于PIREPs 資料的數量，數量越多，評估效果越好。PIREPs資料在反映實際飛機顛簸強度的時候，僅僅依靠機組人員的主觀意識，客觀性弱且機型差異等因素影響較大，這給本文算法閾值的選擇帶來一定干擾，還值得進一步改進和優(yōu)化。