胡康，周江华，*，张晓军

（1.中国科学院空天信息创新研究院，北京 100094； 2.中国科学院大学光电学院，北京 100049）

高空气球作为平流层科学实验的主要载体，具有有效载荷大、发放成本低等显著优点。但是作为一种无动力飞行器，高空气球在驻空阶段只能通过调节浮力和重力来控制垂直方向的运动，而无法主动调节飞行的水平方向，因此，平流层风场水平方向风速分量的数据精度和密度将直接影响高空气球的飞行轨迹［1］。目前，高空气球飞行仿真和实际实验使用的全球风场数据主要来自于欧洲中期天气预报中心（ECWMF）的ERA系列全球气象数据模型，其中最新的ERA5模型提供了全球气象再分析资料和未来5 d内的气象预报数据，每一时刻的数据水平分辨率为0.25°×0.25°（经纬度），包涵从地面到距地面80 km处的37个等压面风场模型［2］。然而在使用ERA5数据对高空气球进行飞行仿真和实际操控的过程中，通过控制气球高度处于合适的风速带来达成区域驻留目标，其有效驻留控制半径在100 km以上，控制精度差，对高空气球实验的轨迹预测、设备回收和空域申请等工作都带来很大程度的影响，因此，需要研究一种根据现有数据生成大范围、高精度、高密度风场数据的插值方法来提高高空气球的控制精度。

现有的风场插值方法可以分成2类：直接插值类方法和融合插值类方法。直接插值类方法是对风场的预报数据或者观测数据采用各种插值方法进行插值，从而获得最终结果。Jardin和Erzberger［3］利用全局多项式插值方法和风场预报对风场格点数据进行了插值。Tao等［4］研究了Hermite插值在风场光谱表示中的误差传递和数值模型。汤新民等［5］研究了3种不同变异模型下的克里格插值法在风场插值中的应用效果。董志南等［6］对地面区域内的实测风场数据进行交叉验证，比较了反距离权重插值、全局多项式插值、局部多项式插值、径向基函数插值、最近邻插值及普通克里金插值6种方法。上述直接插值类方法［3-6］是对风场的预报数据或者观测数据直接进行插值，具有计算速度快、数据质量要求低等优点。这些方法主要面对地面及近地风场，观测的时间和空间尺度范围小且采样难度低，通过密集的风场实测数据，补插出与实际风场高度拟合的插值函数。但是平流层风场时间和空间观测尺度大、采样成本高、数据稀疏且缺乏连续性特征［7］，直接插值类方法难以给出有效的插值结果。

融合插值类方法则先对风场数据进行融合，在融合过程中分析风场特征并设计插值函数，从而对融合结果进行插值，获得最终结果。Gandin［8］采用最优插值法进行了风场融合插值，以最小方差作为约束条件，对观测风场和预报风场进行数据融合，从而获得对风场插值结果的最优估计。朱成阵等［9］通过风场的实测数据和气象预报数据，用内插外推法预测特定位置的风场信息。上述融合插值类方法［8-9］将小范围的实测数据作为观测场，将风场预报数据作为背景场，进行数据融合，分析风场特征，再根据特征进一步设计风场插值模型。但是统计数据表明，风场在时间和空间维度上有显著的差异性［10］，只采用小范围的实测数据进行数据融合获得的插值方法面对更大的区域难以具有良好的泛化性能。同时，计算高精度的风场预报作为模型的背景场需要巨大的算力和时间开销。以国内全球／区域同化预测系统（global／regional assimilation prediction system，GRAPES）中的全球版本为例，作为序贯模型，对某一时刻0.5°×0.5°（经纬度）网格31层模型的积分计算需要200个以上的CPU节点计算10 d［11］，面对需要更高数据密度的高空气球实验，传统的融合插值类方法很难符合其时效性要求。

为了解决平流层风场插值面临的插值精度低、计算时间长和算力开销大等问题，进行速度快、精度高的插值计算，本文提出一种基于地转风（based on geostrophic wind，BGW）模型的平流层水平方向风场插值方法，用2019年7月ERA5全球气象再分析数据对平流层的全球风场进行插值计算实验验证。相比于直接插值类方法，BGW 插值对平流层水平分量的插值效果有明显提升；相比于融合插值类方法，BGW 插值大幅度降低了计算复杂度，有效提高了平流层风场插值的计算速度。

1 基于地转风模型的风场插值

本文提出的BGW 风场插值方法是一种改进的融合插值类方法，主要包含两部分：①为了提升计算效率，用重力位势数据计算地转风模型作为插值背景场；②为了减小背景场误差和观测误差对插值结果的影响，将已有的风速数据作为观测场，用迭代的二元线性回归对背景场数据进行修正，获得最终的插值结果。

1.1 地转风模型改进

一般的融合插值类方法需要使用气象预报模型计算的目标精度风场预报数据作为背景场。气象预报模型一般是基于大气运动方程组的数值模型，能够高精度地预测风速的东向、北向及垂直方向分量，对中短期风场都有着良好的预测效果。然而数值预报模型的实际计算过程需要大量的历史数据、算力和时间作为支撑。因此，当前的融合类插值方法需要大量的时间来进行计算，或依赖于外部输入的目标密度风场预报数据，难以在实验现场根据需求进行实时的离线计算。为了解决这一问题，本文提出一种用地转风模型替代风场预报作为背景场的方法。

地转风模型由白贝罗于1857年提出，是大气运动方程组在大尺度风场条件下推导而出，当大气运动满足科里奥利力和气压梯度力平衡，即地转平衡条件时的理论风［12］，在局地直角坐标系下，其公式如下：

式中：Ω 为地球的旋转角速度，值为7.29×10-5rad／s；φ为纬度。

地转风虽然是根据气压分布计算出的估计值，并非实际存在的风，但是根据实验验证，以地转平衡作为条件计算的中高纬度自由大气风场与真实风场十分接近［13］，因此用地转风模型作为风场插值的背景场能够在提高计算速度的同时，准确地显示中高纬度地区的风速分布特征。

在纬度逼近赤道的过程中，由于科里奥利频率f的值逐渐趋近于0，导致在低纬度地区地转风计算值趋于无穷，因此对f的计算公式进行改进。对于一般意义上的赤道附近区域，即南北纬10°之间的区域，选取边界对科里奥利频率公式进行截断，截断边界的选择范围为（0°，10°］，越接近0则赤道附近区域风场误差越大，越接近10则初值场和准地转运动偏离越多。综合考虑选取中值5°为截断边界，将北纬5°和南纬5°之间的f值固定为2Ωsin 5°，改进后的公式如下：

经过实验验证，相比于使用原始的地转风模型，导致最终插值精度较差，用式（3）的BGW 插值在低纬度地区也能获得良好的插值效果。

1.2 基于二维卷积的地转风快速计算方法

式（1）的地转风模型需要大气密度ρ作为输入，而在实际实验中无法实时获得大气密度数据，因此对式（1）进行改进。另外，为了进一步提升背景场的计算速度，结合GPU的矩阵运算优势改进了背景场计算方法。

1.2.1 基于重力位势的地转风模型

由于气压p是局地直角坐标系（x，y，z，t）的单值函数，可以将坐标系中的高度z用气压p替代，从而对于任意坐标系（x，y，z，t）下的函数F（x，y，z，t）有

式（9）把式（1）中的大气密度和气压这2个变量转换为只有重力位势这一个变量，在降低计算复杂度的同时，解决了地转风模型需要大气密度参与计算的问题。

1.2.2 基于二维卷积的优化计算

式中：“*”表示矩阵的二维卷积运算；A为提取的卷积核。

由于输入是全球重力位势数据的平面展开投影，而卷积计算对矩阵边缘有计算损失，考虑到投影过程是将地球沿180°经线东西向分开，南北两极将纬度90°处展开，在进行卷积计算时将重力位势矩阵左右进行循环扩展，上下进行复制扩展，这种改进有效抑制了卷积计算的边缘损失。

相比于传统模型沿等压线根据大气密度进行逐格计算来求解地转风［14］，BGW 插值通过提取模型卷积核，用二维卷积进行快速的矩阵计算，有效提升了地转风模型的计算速度。

1.3 迭代二元线性回归确定插值函数

BGW 插值方法将地转风模型计算的目标密度风场格点数据作为背景场，用部分格点的实际风速观测数据作为观测场，用观测场对背景场进行修正以获得最终插值结果。每一个格点上的插值结果是由背景场值和修正值用迭代二元线性回归确定，修正值由一定范围内N个格点上的观测值与背景场值的偏差加权获得，最终获得如下线性回归方程式（12）和其约束条件式（13）：

在地理信息领域，传统的距离相关权重矩阵一般采用反距离权重矩阵。反距离权重受到观测点和待插值点的距离及幂次参数控制，需要根据数据密度、搜索范围和插值效果不断调整幂次参数［15］。在权重wi中引入搜索半径参数r并进行重新设计，使得权重可以根据搜索范围进行自适应调整，解决了对不同搜索范围都要重新设计权重的问题，其公式如下：

式中：d为观测点到待插值点的距离。同时为了加快线性回归的收敛速度，对wi进行归一化处理。

为了进一步提高插值精度，在完成一轮线性回归计算后，将插值结果矩阵作为背景场数据再一次输入回归模型进行迭代，重复迭代过程，直到插值精度收敛。经过验证，在两轮迭代时模型的插值精度已经收敛，取两轮迭代的模型作为最终的插值函数。

2 实验验证和对比分析

2.1 数据集

平流层风场插值方法的目的是为了获得尽可能贴近真实风场情况的插值结果，ECWMF的ERA系列气象再分析模型使用了全球气象资料融合系统和完善的气象数据数据库，从1970年至今，对全球各地来自地面雷达测站、船载气象测站、海上浮标、探空气球、飞机和卫星的多源观测资料进行质量评价和数据融合，从而生成高度接近真实风场的完整数据集［2］。以往的高空气球飞行试验也以ERA系列模型作为参考，历史数据也验证了ERA系列模型的可靠性。因此，选择ERA系列最新的ERA5模型2019年7月的气象再分析数据进行插值计算和实验验证。由于高空气球的驻空阶段一般驻留在平流层的底部区域，距地面20 km左右，位于50 hPa大气压等压面附近，选用50 hPa等压面模型的0.5°×0.5°（经纬度）网格的重力位势数据作为输入，用经纬度网格0.25°×0.25°的风速数据进行验证。

2.2 评价指标

通过实验对BGW 插值方法中，背景场对实际风场特征的估计能力、插值模型的优劣程度和插值精度进行评价，分别使用相关系数矩阵、R2检验和平均绝对误差（mean absolute error，MAE）作为评价指标。

2.2.1 相关系数矩阵

背景场对实际风速分布特征的表达能力用背景场数据与再分析资料的相关系数矩阵表示，矩阵X和Y的相关系数矩阵ρXY的公式为

2.3 BGW 插值背景场风场特征估计能力实验

为验证BGW 插值的背景场模型对实际风场特征的估计能力，实验选择ERA 模型2019年7月第一个星期的50 hPa等压面气象数据进行实验，对一周内每日UTC时间8：00的平流层数据，计算其BGW 插值的背景场和再分析风场的相关系数矩阵，计算时将风场数据沿纬度线剪开，每一纬度的数据作为一维特征。由于每一个完整的相关系数矩阵包涵721个相关系数，难以直接分析，取相关系数矩阵的中位数进行分析，结果保留6位有效数字，数据如表1所示。

表1 相关系数矩阵中位数Table 1 Median of correlation coefficient matrix

从表1中可以看出，西风、南风和合风的背景场与再分析风场的相关系数中位数都在0.4以上，属于中等相关和强相关范畴，而且合风的相关系数都在0.6以上，属于强相关，说明背景场数据能够一定程度上表达实际风场的风速分布特征。

牛冲含有优质蛋白质及丰富的维生素 B1、维生素 B2、维生素 C、维生素E，以及雄性激素，具有补肾益精，壮腰膝的功效。常食对于改善性功能有一定作用。

由于表1只对风速的数值特征进行分析，为了更直观地分析风场的几何特征，以ERA5模型2019年7月24日UTC时间8：00的50 hPa等压面气象数据为例，用数据中经纬度网格0.5°×0.5°的重力位势再分析数据作为输入计算经纬度网格0.25°×0.25°的风场水平方向合风速的西风分量u和南风分量v背景场数据，风速数据用热图表示，西风分量和南风分量背景场风速热图分别如图1（a）和图1（b）所示。以ERA5气象再分析模型的0.5°精度西风u和南风v数据作为实际风场进行对比，其风速热图分别如图1（c）和图1（d）所示。

图1中每幅子图的横坐标为经度，纵坐标为纬度。在同一幅西风分量图中，颜色越偏向黄色表示西风越大，颜色越偏向蓝色表示东风越大；在同一幅南风分量图中，颜色越偏向黄色表示南风越大，颜色越偏向蓝色表示北风越大。

从图1（a）中可以发现，模型计算的背景场西风数据在南纬60°附近有一个沿纬度圈环绕地球一周、风速40 m／s以上的大风速西风带，而其他区域的西风风速在沿纬度方向分布较为均匀。从图1（b）可以发现，背景场在南纬60°纬线圈上西经60°和西经140°附近有2处南风超过20 m／s的集中分布区域。而从实际风场图1（c）和图1（d）中可以看出，背景场估计的西风分量和南风分量的大风速带和风速集中分布区域的主要特征与观测场基本一致，结合相关系数分析结果说明，用改进的地转风模型作为插值背景场能够替代风场预报数据，对风速分布特征进行估计。

图1 风速背景场和实际风场数据Fig.1 Background field and actual field of wind speed data

2.4 平流层风场插值实验

用搜索半径分别为30 km、60 km和90 km的BGW 插值方法与直接插值类方法中使用最多的最近邻域插值法（nearest neighbor，NN）和反距离权重插值法（inverse distance weighted，IDW）进行对比，IDW 的幂次参数取2，比较几种方法在2019年7月每天24个50 hPa全球风场数据上的模型插值性能，取R2检验结果的均值作为评价指标，结果保留6位有效数字，数据如表2所示。

表2 不同插值方法的R2 检验结果Table 2 R2 square test results of different interpolation methods

对于数据网格间距为0.5°（经纬度）的输入数据，待插值点与观测点的最远平均距离约为20 km，从表2的R2检验结果可以发现，BGW 插值的搜索半径为30 km时模型性能最好，该结果验证了1.3节中搜索半径的最佳参数范围应该在待插值点与观测点最远距离的一倍到两倍之间的结论。之后的实验BGW 插值默认搜索半径选择30 km。

表3 不同插值方法的平均绝对误差Table 3 Mean absolute error of different interpolation methods

从表3中可以发现，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系统下的Python代码进行算法实现，用英特尔7代i7CPU进行计算的平均耗时为1.37 s，用英伟达GTX1050GPU进行同样工作的平均耗时为1.26 s。而且在实际的高空气球飞行仿真和外场实验中，需要对风场进行持续、实时的插值计算，GPU在矩阵计算上的速度优势在面对大量矩阵计算时还会进一步扩大。另外，一般的融合类插值背景场采用气象预报模型，无法对每一层等压面风场数据解耦进行单独计算，而完整的预测模型需要大量的算力设备支持和天级的运算时间，BGW 插值有效解决了融合插值类方法背景场计算时间长、算力开销大的问题。

实验结果表明，相比于直接插值类方法，BGW 插值对一维的重力位势矢量进行插值，再用观测风速对背景场进行修正，误差绝对值要小于对三维的观测风速矢量进行直接插值，有效减小了插值方法的误差。而相比于融合插值类方法，BGW 插值是用当前时刻的重力位势数据结合地转风模型来计算风速的估计值作为背景场，与直接引入预报数据的方法相比，BGW 插值无需引入外部目标精度的气象预报数据，可以只通过已有的数据进行离线运算；与其他自主计算风场预报的方法相比，由于预报模型是一个高度耦合的数值系统，无法对风速的水平方向分量单独进行预测，导致其模型的运算时间和算力消耗极大，而BGW 插值能只对单一等压面的背景场进行计算，且计算速度在秒级，有效提高了模型的计算效率。

3 结 论

面对稀疏的平流层风场数据，基于地转风模型和迭代二元线性回归提出一种针对全球平流层风场水平分量数据的BGW 插值方法，有效提升了计算速度和精度。

1）提出了一种适用于全纬度地区地转风模型计算的科里奥利频率改进公式，同时推导出一种适用于高空气球实验且能有效表达真实风场特征的地转风快速计算方法。

2）通过迭代的二元线性回归和改进的自适应观测场偏差权重矩阵减小了观测场误差对插值结果的影响，并解决了面对不同密度的风场数据都要重新设计权重矩阵的问题。

3）实验结果表明，与2种直接插值类方法相比，NN插值的平均绝对误差约为0.21 m／s，IDW插值的平均绝对误差约为0.43 m／s，BGW 插值将风速平均绝对误差减小到约0.13 m／s，有效提升了插值精度；与融合插值类方法相比，一般使用的背景场模型需要10 d以上的计算时间，BGW 插值对单一等压面的风场模型的插值计算速度则提升到1.26 s，有效提升了插值速度。

为使提出的平流层风场插值方法能对高空气球仿真和实验提供进一步支持，下一阶段工作应将目前的风场插值方法结合大气位温方程，研究时间维度上的风场插值方法。