胡康，周江華，*，張曉軍

（1.中國科學院空天信息創(chuàng)新研究院，北京 100094； 2.中國科學院大學光電學院，北京 100049）

高空氣球作為平流層科學實驗的主要載體，具有有效載荷大、發(fā)放成本低等顯著優(yōu)點。但是作為一種無動力飛行器，高空氣球在駐空階段只能通過調(diào)節(jié)浮力和重力來控制垂直方向的運動，而無法主動調(diào)節(jié)飛行的水平方向，因此，平流層風場水平方向風速分量的數(shù)據(jù)精度和密度將直接影響高空氣球的飛行軌跡［1］。目前，高空氣球飛行仿真和實際實驗使用的全球風場數(shù)據(jù)主要來自于歐洲中期天氣預報中心（ECWMF）的ERA系列全球氣象數(shù)據(jù)模型，其中最新的ERA5模型提供了全球氣象再分析資料和未來5 d內(nèi)的氣象預報數(shù)據(jù)，每一時刻的數(shù)據(jù)水平分辨率為0.25°×0.25°（經(jīng)緯度），包涵從地面到距地面80 km處的37個等壓面風場模型［2］。然而在使用ERA5數(shù)據(jù)對高空氣球進行飛行仿真和實際操控的過程中，通過控制氣球高度處于合適的風速帶來達成區(qū)域駐留目標，其有效駐留控制半徑在100 km以上，控制精度差，對高空氣球?qū)嶒灥能壽E預測、設備回收和空域申請等工作都帶來很大程度的影響，因此，需要研究一種根據(jù)現(xiàn)有數(shù)據(jù)生成大范圍、高精度、高密度風場數(shù)據(jù)的插值方法來提高高空氣球的控制精度。

現(xiàn)有的風場插值方法可以分成2類：直接插值類方法和融合插值類方法。直接插值類方法是對風場的預報數(shù)據(jù)或者觀測數(shù)據(jù)采用各種插值方法進行插值，從而獲得最終結(jié)果。Jardin和Erzberger［3］利用全局多項式插值方法和風場預報對風場格點數(shù)據(jù)進行了插值。Tao等［4］研究了Hermite插值在風場光譜表示中的誤差傳遞和數(shù)值模型。湯新民等［5］研究了3種不同變異模型下的克里格插值法在風場插值中的應用效果。董志南等［6］對地面區(qū)域內(nèi)的實測風場數(shù)據(jù)進行交叉驗證，比較了反距離權(quán)重插值、全局多項式插值、局部多項式插值、徑向基函數(shù)插值、最近鄰插值及普通克里金插值6種方法。上述直接插值類方法［3-6］是對風場的預報數(shù)據(jù)或者觀測數(shù)據(jù)直接進行插值，具有計算速度快、數(shù)據(jù)質(zhì)量要求低等優(yōu)點。這些方法主要面對地面及近地風場，觀測的時間和空間尺度范圍小且采樣難度低，通過密集的風場實測數(shù)據(jù)，補插出與實際風場高度擬合的插值函數(shù)。但是平流層風場時間和空間觀測尺度大、采樣成本高、數(shù)據(jù)稀疏且缺乏連續(xù)性特征［7］，直接插值類方法難以給出有效的插值結(jié)果。

融合插值類方法則先對風場數(shù)據(jù)進行融合，在融合過程中分析風場特征并設計插值函數(shù)，從而對融合結(jié)果進行插值，獲得最終結(jié)果。Gandin［8］采用最優(yōu)插值法進行了風場融合插值，以最小方差作為約束條件，對觀測風場和預報風場進行數(shù)據(jù)融合，從而獲得對風場插值結(jié)果的最優(yōu)估計。朱成陣等［9］通過風場的實測數(shù)據(jù)和氣象預報數(shù)據(jù)，用內(nèi)插外推法預測特定位置的風場信息。上述融合插值類方法［8-9］將小范圍的實測數(shù)據(jù)作為觀測場，將風場預報數(shù)據(jù)作為背景場，進行數(shù)據(jù)融合，分析風場特征，再根據(jù)特征進一步設計風場插值模型。但是統(tǒng)計數(shù)據(jù)表明，風場在時間和空間維度上有顯著的差異性［10］，只采用小范圍的實測數(shù)據(jù)進行數(shù)據(jù)融合獲得的插值方法面對更大的區(qū)域難以具有良好的泛化性能。同時，計算高精度的風場預報作為模型的背景場需要巨大的算力和時間開銷。以國內(nèi)全球／區(qū)域同化預測系統(tǒng)（global／regional assimilation prediction system，GRAPES）中的全球版本為例，作為序貫模型，對某一時刻0.5°×0.5°（經(jīng)緯度）網(wǎng)格31層模型的積分計算需要200個以上的CPU節(jié)點計算10 d［11］，面對需要更高數(shù)據(jù)密度的高空氣球?qū)嶒灒瑐鹘y(tǒng)的融合插值類方法很難符合其時效性要求。

為了解決平流層風場插值面臨的插值精度低、計算時間長和算力開銷大等問題，進行速度快、精度高的插值計算，本文提出一種基于地轉(zhuǎn)風（based on geostrophic wind，BGW）模型的平流層水平方向風場插值方法，用2019年7月ERA5全球氣象再分析數(shù)據(jù)對平流層的全球風場進行插值計算實驗驗證。相比于直接插值類方法，BGW 插值對平流層水平分量的插值效果有明顯提升；相比于融合插值類方法，BGW 插值大幅度降低了計算復雜度，有效提高了平流層風場插值的計算速度。

1 基于地轉(zhuǎn)風模型的風場插值

本文提出的BGW 風場插值方法是一種改進的融合插值類方法，主要包含兩部分：①為了提升計算效率，用重力位勢數(shù)據(jù)計算地轉(zhuǎn)風模型作為插值背景場；②為了減小背景場誤差和觀測誤差對插值結(jié)果的影響，將已有的風速數(shù)據(jù)作為觀測場，用迭代的二元線性回歸對背景場數(shù)據(jù)進行修正，獲得最終的插值結(jié)果。

1.1 地轉(zhuǎn)風模型改進

一般的融合插值類方法需要使用氣象預報模型計算的目標精度風場預報數(shù)據(jù)作為背景場。氣象預報模型一般是基于大氣運動方程組的數(shù)值模型，能夠高精度地預測風速的東向、北向及垂直方向分量，對中短期風場都有著良好的預測效果。然而數(shù)值預報模型的實際計算過程需要大量的歷史數(shù)據(jù)、算力和時間作為支撐。因此，當前的融合類插值方法需要大量的時間來進行計算，或依賴于外部輸入的目標密度風場預報數(shù)據(jù)，難以在實驗現(xiàn)場根據(jù)需求進行實時的離線計算。為了解決這一問題，本文提出一種用地轉(zhuǎn)風模型替代風場預報作為背景場的方法。

地轉(zhuǎn)風模型由白貝羅于1857年提出，是大氣運動方程組在大尺度風場條件下推導而出，當大氣運動滿足科里奧利力和氣壓梯度力平衡，即地轉(zhuǎn)平衡條件時的理論風［12］，在局地直角坐標系下，其公式如下：

式中：Ω 為地球的旋轉(zhuǎn)角速度，值為7.29×10-5rad／s；φ為緯度。

地轉(zhuǎn)風雖然是根據(jù)氣壓分布計算出的估計值，并非實際存在的風，但是根據(jù)實驗驗證，以地轉(zhuǎn)平衡作為條件計算的中高緯度自由大氣風場與真實風場十分接近［13］，因此用地轉(zhuǎn)風模型作為風場插值的背景場能夠在提高計算速度的同時，準確地顯示中高緯度地區(qū)的風速分布特征。

在緯度逼近赤道的過程中，由于科里奧利頻率f的值逐漸趨近于0，導致在低緯度地區(qū)地轉(zhuǎn)風計算值趨于無窮，因此對f的計算公式進行改進。對于一般意義上的赤道附近區(qū)域，即南北緯10°之間的區(qū)域，選取邊界對科里奧利頻率公式進行截斷，截斷邊界的選擇范圍為（0°，10°］，越接近0則赤道附近區(qū)域風場誤差越大，越接近10則初值場和準地轉(zhuǎn)運動偏離越多。綜合考慮選取中值5°為截斷邊界，將北緯5°和南緯5°之間的f值固定為2Ωsin 5°，改進后的公式如下：

經(jīng)過實驗驗證，相比于使用原始的地轉(zhuǎn)風模型，導致最終插值精度較差，用式（3）的BGW 插值在低緯度地區(qū)也能獲得良好的插值效果。

1.2 基于二維卷積的地轉(zhuǎn)風快速計算方法

式（1）的地轉(zhuǎn)風模型需要大氣密度ρ作為輸入，而在實際實驗中無法實時獲得大氣密度數(shù)據(jù)，因此對式（1）進行改進。另外，為了進一步提升背景場的計算速度，結(jié)合GPU的矩陣運算優(yōu)勢改進了背景場計算方法。

1.2.1 基于重力位勢的地轉(zhuǎn)風模型

由于氣壓p是局地直角坐標系（x，y，z，t）的單值函數(shù)，可以將坐標系中的高度z用氣壓p替代，從而對于任意坐標系（x，y，z，t）下的函數(shù)F（x，y，z，t）有

式（9）把式（1）中的大氣密度和氣壓這2個變量轉(zhuǎn)換為只有重力位勢這一個變量，在降低計算復雜度的同時，解決了地轉(zhuǎn)風模型需要大氣密度參與計算的問題。

1.2.2 基于二維卷積的優(yōu)化計算

式中：“*”表示矩陣的二維卷積運算；A為提取的卷積核。

由于輸入是全球重力位勢數(shù)據(jù)的平面展開投影，而卷積計算對矩陣邊緣有計算損失，考慮到投影過程是將地球沿180°經(jīng)線東西向分開，南北兩極將緯度90°處展開，在進行卷積計算時將重力位勢矩陣左右進行循環(huán)擴展，上下進行復制擴展，這種改進有效抑制了卷積計算的邊緣損失。

相比于傳統(tǒng)模型沿等壓線根據(jù)大氣密度進行逐格計算來求解地轉(zhuǎn)風［14］，BGW 插值通過提取模型卷積核，用二維卷積進行快速的矩陣計算，有效提升了地轉(zhuǎn)風模型的計算速度。

1.3 迭代二元線性回歸確定插值函數(shù)

BGW 插值方法將地轉(zhuǎn)風模型計算的目標密度風場格點數(shù)據(jù)作為背景場，用部分格點的實際風速觀測數(shù)據(jù)作為觀測場，用觀測場對背景場進行修正以獲得最終插值結(jié)果。每一個格點上的插值結(jié)果是由背景場值和修正值用迭代二元線性回歸確定，修正值由一定范圍內(nèi)N個格點上的觀測值與背景場值的偏差加權(quán)獲得，最終獲得如下線性回歸方程式（12）和其約束條件式（13）：

在地理信息領(lǐng)域，傳統(tǒng)的距離相關(guān)權(quán)重矩陣一般采用反距離權(quán)重矩陣。反距離權(quán)重受到觀測點和待插值點的距離及冪次參數(shù)控制，需要根據(jù)數(shù)據(jù)密度、搜索范圍和插值效果不斷調(diào)整冪次參數(shù)［15］。在權(quán)重wi中引入搜索半徑參數(shù)r并進行重新設計，使得權(quán)重可以根據(jù)搜索范圍進行自適應調(diào)整，解決了對不同搜索范圍都要重新設計權(quán)重的問題，其公式如下：

式中：d為觀測點到待插值點的距離。同時為了加快線性回歸的收斂速度，對wi進行歸一化處理。

為了進一步提高插值精度，在完成一輪線性回歸計算后，將插值結(jié)果矩陣作為背景場數(shù)據(jù)再一次輸入回歸模型進行迭代，重復迭代過程，直到插值精度收斂。經(jīng)過驗證，在兩輪迭代時模型的插值精度已經(jīng)收斂，取兩輪迭代的模型作為最終的插值函數(shù)。

2 實驗驗證和對比分析

2.1 數(shù)據(jù)集

平流層風場插值方法的目的是為了獲得盡可能貼近真實風場情況的插值結(jié)果，ECWMF的ERA系列氣象再分析模型使用了全球氣象資料融合系統(tǒng)和完善的氣象數(shù)據(jù)數(shù)據(jù)庫，從1970年至今，對全球各地來自地面雷達測站、船載氣象測站、海上浮標、探空氣球、飛機和衛(wèi)星的多源觀測資料進行質(zhì)量評價和數(shù)據(jù)融合，從而生成高度接近真實風場的完整數(shù)據(jù)集［2］。以往的高空氣球飛行試驗也以ERA系列模型作為參考，歷史數(shù)據(jù)也驗證了ERA系列模型的可靠性。因此，選擇ERA系列最新的ERA5模型2019年7月的氣象再分析數(shù)據(jù)進行插值計算和實驗驗證。由于高空氣球的駐空階段一般駐留在平流層的底部區(qū)域，距地面20 km左右，位于50 hPa大氣壓等壓面附近，選用50 hPa等壓面模型的0.5°×0.5°（經(jīng)緯度）網(wǎng)格的重力位勢數(shù)據(jù)作為輸入，用經(jīng)緯度網(wǎng)格0.25°×0.25°的風速數(shù)據(jù)進行驗證。

2.2 評價指標

通過實驗對BGW 插值方法中，背景場對實際風場特征的估計能力、插值模型的優(yōu)劣程度和插值精度進行評價，分別使用相關(guān)系數(shù)矩陣、R2檢驗和平均絕對誤差（mean absolute error，MAE）作為評價指標。

2.2.1 相關(guān)系數(shù)矩陣

背景場對實際風速分布特征的表達能力用背景場數(shù)據(jù)與再分析資料的相關(guān)系數(shù)矩陣表示，矩陣X和Y的相關(guān)系數(shù)矩陣ρXY的公式為

2.3 BGW 插值背景場風場特征估計能力實驗

為驗證BGW 插值的背景場模型對實際風場特征的估計能力，實驗選擇ERA 模型2019年7月第一個星期的50 hPa等壓面氣象數(shù)據(jù)進行實驗，對一周內(nèi)每日UTC時間8：00的平流層數(shù)據(jù)，計算其BGW 插值的背景場和再分析風場的相關(guān)系數(shù)矩陣，計算時將風場數(shù)據(jù)沿緯度線剪開，每一緯度的數(shù)據(jù)作為一維特征。由于每一個完整的相關(guān)系數(shù)矩陣包涵721個相關(guān)系數(shù)，難以直接分析，取相關(guān)系數(shù)矩陣的中位數(shù)進行分析，結(jié)果保留6位有效數(shù)字，數(shù)據(jù)如表1所示。

表1 相關(guān)系數(shù)矩陣中位數(shù)Table 1 Median of correlation coefficient matrix

從表1中可以看出，西風、南風和合風的背景場與再分析風場的相關(guān)系數(shù)中位數(shù)都在0.4以上，屬于中等相關(guān)和強相關(guān)范疇，而且合風的相關(guān)系數(shù)都在0.6以上，屬于強相關(guān)，說明背景場數(shù)據(jù)能夠一定程度上表達實際風場的風速分布特征。

牛沖含有優(yōu)質(zhì)蛋白質(zhì)及豐富的維生素 B1、維生素 B2、維生素 C、維生素E，以及雄性激素，具有補腎益精，壯腰膝的功效。常食對于改善性功能有一定作用。

由于表1只對風速的數(shù)值特征進行分析，為了更直觀地分析風場的幾何特征，以ERA5模型2019年7月24日UTC時間8：00的50 hPa等壓面氣象數(shù)據(jù)為例，用數(shù)據(jù)中經(jīng)緯度網(wǎng)格0.5°×0.5°的重力位勢再分析數(shù)據(jù)作為輸入計算經(jīng)緯度網(wǎng)格0.25°×0.25°的風場水平方向合風速的西風分量u和南風分量v背景場數(shù)據(jù)，風速數(shù)據(jù)用熱圖表示，西風分量和南風分量背景場風速熱圖分別如圖1（a）和圖1（b）所示。以ERA5氣象再分析模型的0.5°精度西風u和南風v數(shù)據(jù)作為實際風場進行對比，其風速熱圖分別如圖1（c）和圖1（d）所示。

圖1中每幅子圖的橫坐標為經(jīng)度，縱坐標為緯度。在同一幅西風分量圖中，顏色越偏向黃色表示西風越大，顏色越偏向藍色表示東風越大；在同一幅南風分量圖中，顏色越偏向黃色表示南風越大，顏色越偏向藍色表示北風越大。

從圖1（a）中可以發(fā)現(xiàn)，模型計算的背景場西風數(shù)據(jù)在南緯60°附近有一個沿緯度圈環(huán)繞地球一周、風速40 m／s以上的大風速西風帶，而其他區(qū)域的西風風速在沿緯度方向分布較為均勻。從圖1（b）可以發(fā)現(xiàn)，背景場在南緯60°緯線圈上西經(jīng)60°和西經(jīng)140°附近有2處南風超過20 m／s的集中分布區(qū)域。而從實際風場圖1（c）和圖1（d）中可以看出，背景場估計的西風分量和南風分量的大風速帶和風速集中分布區(qū)域的主要特征與觀測場基本一致，結(jié)合相關(guān)系數(shù)分析結(jié)果說明，用改進的地轉(zhuǎn)風模型作為插值背景場能夠替代風場預報數(shù)據(jù)，對風速分布特征進行估計。

圖1 風速背景場和實際風場數(shù)據(jù)Fig.1 Background field and actual field of wind speed data

2.4 平流層風場插值實驗

用搜索半徑分別為30 km、60 km和90 km的BGW 插值方法與直接插值類方法中使用最多的最近鄰域插值法（nearest neighbor，NN）和反距離權(quán)重插值法（inverse distance weighted，IDW）進行對比，IDW 的冪次參數(shù)取2，比較幾種方法在2019年7月每天24個50 hPa全球風場數(shù)據(jù)上的模型插值性能，取R2檢驗結(jié)果的均值作為評價指標，結(jié)果保留6位有效數(shù)字，數(shù)據(jù)如表2所示。

表2 不同插值方法的R2 檢驗結(jié)果Table 2 R2 square test results of different interpolation methods

對于數(shù)據(jù)網(wǎng)格間距為0.5°（經(jīng)緯度）的輸入數(shù)據(jù)，待插值點與觀測點的最遠平均距離約為20 km，從表2的R2檢驗結(jié)果可以發(fā)現(xiàn)，BGW 插值的搜索半徑為30 km時模型性能最好，該結(jié)果驗證了1.3節(jié)中搜索半徑的最佳參數(shù)范圍應該在待插值點與觀測點最遠距離的一倍到兩倍之間的結(jié)論。之后的實驗BGW 插值默認搜索半徑選擇30 km。

表3 不同插值方法的平均絕對誤差Table 3 Mean absolute error of different interpolation methods

從表3中可以發(fā)現(xiàn)，BGW 插值在西風、南風和合風上的MAE都要小于NN和IDW 插值。在西風上，BGW 插值誤差是 NN 插值誤差的42.32%，是IDW 插值誤差的20.77%；在南風上，BGW 插值誤差是NN插值誤差的48.41%，是IDW 插值誤差的85.68%；在合風上，BGW 插值誤差是NN插值誤差的59.26%，是IDW 插值誤差的29.16%。

同時對BGW 插值方法的計算速度進行實驗，對于單等壓面風場的平面插值問題，用Windows系統(tǒng)下的Python代碼進行算法實現(xiàn)，用英特爾7代i7CPU進行計算的平均耗時為1.37 s，用英偉達GTX1050GPU進行同樣工作的平均耗時為1.26 s。而且在實際的高空氣球飛行仿真和外場實驗中，需要對風場進行持續(xù)、實時的插值計算，GPU在矩陣計算上的速度優(yōu)勢在面對大量矩陣計算時還會進一步擴大。另外，一般的融合類插值背景場采用氣象預報模型，無法對每一層等壓面風場數(shù)據(jù)解耦進行單獨計算，而完整的預測模型需要大量的算力設備支持和天級的運算時間，BGW 插值有效解決了融合插值類方法背景場計算時間長、算力開銷大的問題。

實驗結(jié)果表明，相比于直接插值類方法，BGW 插值對一維的重力位勢矢量進行插值，再用觀測風速對背景場進行修正，誤差絕對值要小于對三維的觀測風速矢量進行直接插值，有效減小了插值方法的誤差。而相比于融合插值類方法，BGW 插值是用當前時刻的重力位勢數(shù)據(jù)結(jié)合地轉(zhuǎn)風模型來計算風速的估計值作為背景場，與直接引入預報數(shù)據(jù)的方法相比，BGW 插值無需引入外部目標精度的氣象預報數(shù)據(jù)，可以只通過已有的數(shù)據(jù)進行離線運算；與其他自主計算風場預報的方法相比，由于預報模型是一個高度耦合的數(shù)值系統(tǒng)，無法對風速的水平方向分量單獨進行預測，導致其模型的運算時間和算力消耗極大，而BGW 插值能只對單一等壓面的背景場進行計算，且計算速度在秒級，有效提高了模型的計算效率。

3 結(jié) 論

面對稀疏的平流層風場數(shù)據(jù)，基于地轉(zhuǎn)風模型和迭代二元線性回歸提出一種針對全球平流層風場水平分量數(shù)據(jù)的BGW 插值方法，有效提升了計算速度和精度。

1）提出了一種適用于全緯度地區(qū)地轉(zhuǎn)風模型計算的科里奧利頻率改進公式，同時推導出一種適用于高空氣球?qū)嶒炃夷苡行П磉_真實風場特征的地轉(zhuǎn)風快速計算方法。

2）通過迭代的二元線性回歸和改進的自適應觀測場偏差權(quán)重矩陣減小了觀測場誤差對插值結(jié)果的影響，并解決了面對不同密度的風場數(shù)據(jù)都要重新設計權(quán)重矩陣的問題。

3）實驗結(jié)果表明，與2種直接插值類方法相比，NN插值的平均絕對誤差約為0.21 m／s，IDW插值的平均絕對誤差約為0.43 m／s，BGW 插值將風速平均絕對誤差減小到約0.13 m／s，有效提升了插值精度；與融合插值類方法相比，一般使用的背景場模型需要10 d以上的計算時間，BGW 插值對單一等壓面的風場模型的插值計算速度則提升到1.26 s，有效提升了插值速度。

為使提出的平流層風場插值方法能對高空氣球仿真和實驗提供進一步支持，下一階段工作應將目前的風場插值方法結(jié)合大氣位溫方程，研究時間維度上的風場插值方法。