基于地心坐标系的卫星导航测风平滑算法

2022-04-28李宛桐姜明史静杨荣康黄威

气象科技 2022年2期

李宛桐姜明史静杨荣康黄威

(1 天津市气象探测中心，天津 300061；2 中国气象局气象探测中心，北京 100081；3 中国洛阳电子装备试验中心，济源 459000)

引言

20世纪80年代以来，出现了奥米伽、罗兰-C以及全球定位系统(Global Positioning System, GPS)、北斗等导航测风系统[1]，随着北斗卫星导航系统的发展以及考虑探空经济成本和加密观测的进一步需求，我国探空业务将从二次雷达测风逐步过渡到卫星导航测风[2]。与此同时，国产卫星导航探空系统的数据采集率、准确度以及系统的自动化程度有很大提升，其数据采集率已由原59-701探空系统的每分钟1组提高到每秒钟1组[3-4]，而现行《常规高空气象观测业务规范》(以下简称《规范》)规定的业务应用数据和文件计算数学模型还停留在以分钟间隔处理的基础上[5]。分钟风的计算方式浪费了大量数据资源[6]，不能完全发挥卫星导航探空系统的效益，但是利用秒间隔定位数据直接计算的原始高空风脉动太大，不能用于生成业务应用文件。

关于原始高空风数据的处理，在目前业务应用的L波段探空系统中，已有很多专家论证采用矢量滑动平均法，并对滑动平均窗口提出了建议：王缅等[7]建议采用窗口为1 min的滑动平均方式或者前20 min采用30 s，以后采用1 min的分段滑动平均方法；梁建平等[8]提出使用30～45 s的滑动平均窗口；陈磊等[9]认为在全程使用50～60 s窗口或前50 min使用30～40 s窗口，50 min以后使用50～60 s窗口结果较好。在卫星导航测风相关研究中，也有不少学者使用矢量滑动平均法处理原始高空风数据：张国舫等[10]在比较GPS探空仪多普勒效应方法和探空气球定位方法计算高空风时，采用20 s作为矢量滑动平均窗口；姚雯等[11]在评估国产GPS探空仪秒级探空数据随机误差时，对秒间隔高空风数据进行了30点和15点滑动平均，并证明平滑程度的差异对间接计算风的标准差影响较大。

上述研究在利用秒间隔定位数据计算原始高空风时，均参照《规范》规定的量得风层计算方法，即基于站心坐标系进行计算：首先将气球的大地坐标垂直投影到站心坐标平面上，然后用该平面上两相邻投影点的坐标数据计算风向风速，没有考虑地球曲率的影响。风是平行于地面的空气流动，气球测风是将气球看作随气流移动的质点，气球在升空过程中同时受上升力和地球引力的影响，基于两种力的作用，其实际的飞行轨迹，是在上升和随空气的水平运动向远方漂移同时，对于地心坐标系不断下降，始终沿地球曲面飞行，即气球的运动是基于地心坐标系的[12]。此外，根据天气学的要求，规定高度风、规定标准气压层风和最大风层等需要以量得风层为基础内插计算的高空风资料，其高度是以垂直于地球表面的位势高度表示的[13]，即业务上分析量得风层时所使用的站心坐标系与使用高空风资料的坐标系并不一致，造成这种不一致的原因与59-701系统下风向风速以人工计算为主有关[14]。如今计算机广泛应用，卫星导航探空系统提供的气球定位又是大地坐标，应在地心坐标系中计算高空风[15]。

本文采用2013年12月中国气象局在阳江探空站进行北斗-GPS双模式探空仪试验的实际施放数据，分析基于地心坐标系的卫星导航测风矢量滑动平均窗口选择，为探空业务高空风计算方法提供参考。

1 高空风数据处理

1.1 坐标系统

在卫星导航系统中，坐标系是描述卫星运动、处理观测数据和求解接收机位置的数学和物理计算基础[16]。利用坐标系可以方便地为卫星导航建立数学公式，当参考坐标系确立后，卫星和接收机的状态就很容易用数学模型表示出来[17]。根据应用场合的不同，选用的坐标系也不相同。

1.1.1 地心大地坐标系

目前，GPS所采用的坐标系是世界大地坐标系(World Geodetic System of 1984, WGS-84)，GPS所有的星历参数就是基于此坐标系[18]。WGS-84坐标系的原点位于地球质心，其Z轴指向国际时间局(BIH)1984年00:00定义的BIH1984.0协议地球极(Conventional Terrestrial Pole, CTP)方向，X轴指向BIH1984.0的零子午面与CTP赤道的交点，X轴、Y轴、Z轴构成右手坐标系，如图1所示。

图1 世界大地坐标系(WGS-84)

北斗所使用的坐标系是2000国家大地坐标系(China Geodetic Coordinate System 2000, CGCS2000)，自2008年7月1日开始启用，是全球地心坐标系在我国的具体体现。CGCS坐标系的原点为地球质心，初始定向由1984.0的BIH定向给定，CGCS2000的参考历元为2000.0[19]。

CGCS2000与WGS-84坐标系在原点、尺度、定向及定向的定义基本相同，参考椭球非常相近，只有扁率f有微小的差异，鉴于在坐标系定义和实现上的比较，可以认为CGCS2000与WGS-84是相容的；在坐标系的实现精度范围内，CGCS2000与WGS-84是一致的[20-22]。

在地心大地坐标系中，点用(L,B,H)来表示。其中，L表示经度，B表示纬度，H表示高程。

1.1.2 地心地固坐标系

地心地固坐标系(Earth-Centered Earth-Fixed, ECEF)简称地心坐标系，是以地球质心为原点的笛卡儿坐标系。X轴指向0°大地子午线与赤道的交点，Y轴指向东经90°大地子午线与赤道的交点，Z轴指向协议地球北极，组成右手系直角坐标，如图2所示。

图2 地心坐标系(ECEF)

在用户接收机中，因为利用了ECEF坐标系表示卫星位置(作为定位参考点)，所以计算得到的用户位置一般也表征在ECEF坐标系中。

1.1.3 坐标转换

卫星导航探空系统测风时，设某一时刻气球的定位数据Pi(Li,Bi,Hi)。在ECEF坐标系中，Li为ZOX平面与ZOPi平面的夹角，自ZOX平面起算右旋为正；Bi为过Pi点的椭球面法线与XOY平面的夹角，自XOY平面向z轴方向量取为正；Hi为过Pi点的椭球面法线自椭球面至气球的距离，以远离椭球面中心方向为正，此法线和z轴的交点与地心通常并不重合。据此，气球某一时刻在ECEF坐标系中的位置Pi(Xi,Yi,Zi)可表示如下：

(1)

式中，Ni为卯酉圈曲率半径；e为地球椭球第一偏心率，用下式计算：

(2)

(3)

其中，a=6378137 m为地球的长半轴长度，b=6356752.3142 m为地球的短半轴长度。

1.2 基于地心坐标系计算原始高空风

参考李浩等的研究成果[15]，基于地心坐标系的卫星导航探空系统原始高空风计算方法如下：

将相邻时刻气球位置记为Pi-1(Li-1,Bi-1,Hi-1)和Pi(Li,Bi,Hi)，选取高度为(Hi-1+Hi)/2的椭球面作为计算高空风的水平面，把该水平面叫做投影椭球面，如图3所示。图3中已标出投影椭球面上的经线和纬线。

图3 投影椭球面

把Pi-1、Pi分别沿卯酉圈曲率半径方向向上、向下投影到该投影椭球面上，形成的投影点分别记作Pro,i-1(Li-1,Bi-1,(Hi-1+Hi)/2)和Pro,i(Li,Bi,(Hi-1+Hi)/2)，把Pro,i-1沿着Bi-1纬线投影到Li经线上，投影点记作N；把Pro,i-1沿着Li-1经线投影到Bi纬线上，投影点记作M。显然，N、M的坐标分别是(Li,Bi-1,(Hi-1+Hi)/2)、(Li-1,Bi,(Hi-1+Hi)/2)。

(4)

注意到dPN≠dMP，定义(dPN+dMP)/2为气球在纬圈(即东西方向)上移动的水平位移，显然

(5)

将气球经向运动速度记作vx，即南北风速分量；纬向运动速度记作vy，即东西风速分量。有

(6)

对应风速V和风向D：

(7)

(8)

1.3 原始高空风数据平滑处理

本文采用矢量滑动平均法对卫星导航探空系统原始高空风数据进行处理。矢量滑动平均法要求设置平均时段(窗口)和步进时间，步进时间通常以原始数据的时间间隔为准。利用公式(9)对每一平均时段风速分量进行平滑处理，所得平均值赋给该平均时段的中间时。

(9)

式中，n为平均时段内风速分量的个数。

计算完成第一个时段的两个风速分量以后，在原来时段的前面去掉一个步进时间的数据，在其后面增加一个步进时间的数据，再进行同样的计算，直至最后一组数据。这样就得到规定窗口的两个风速分量的滑动平均值。由于计算得到的风速分量起始时间为第一个平均时段的中间时刻，在这一时刻之前的数据可以通过与地面风速分量数据内插获得，最终得到从地面至球炸的连续数据。

2 实例数据及分析方法

本文采用2013年12月中国气象局在阳江探空站进行北斗-GPS双模式探空仪试验的实际施放数据。试验采取同球施放方法，即在一个探空气球上同时悬挂一个RS92探空仪与两个北斗-GPS双模式探空仪，北斗-GPS双模式探空仪分别以北斗-GPS综合模式、单GPS模式和单北斗模式解算秒间隔定位数据。

3 矢量滑动平均窗口选择

3.1 窗口对测风结果的影响

根据上述公式计算全部探空数据10 s、20 s、30 s、40 s、50 s、60 s滑动平均风，并与RS92探空系统给出的风速风向相比较，其结果如下。

图4是比对结果的典型例子，图中蓝色曲线为北斗-GPS探空系统滑动平均风速风向值。从图中可以看出，平滑后的风速、风向曲线与RS92探空系统趋势基本一致。和RS92测风数据相比，10 s滑动平均风的曲线波动最为明显，随着窗口的延长，曲线变化幅度衰减得越大，变化频率越低。由于窗口较短、波动较大，10 s和20 s滑动平均风的风向在升空扰动明显处，如遇过北情况时偏差较大。当窗口为40 s及以上时，滑动平均风风速、风向曲线趋于平滑，摆动幅度逐渐减小，同时产生相位滞后，其滞后时间随窗口的延长而增加。

图4 2013年12月18日不同滑动平均窗口北斗-GPS探空系统滑动平均(蓝色线)与RS92探空系统探测(红色线)风向风速变化:(a1、a2)10 s，(b1、b2)20 s，(c1、c2)30 s，(d1、d2)40 s，(e1、e2)50 s，(f1、f2)60 s

3.2 误差分析

统计窗口分别为10 s、12 s、14 s、…、56 s、58 s、60 s的滑动平均风在每千米高度层上风速分量的误差，绘制风速分量系统误差变化曲线和95%置信概率误差范围变化图，如图5和图6所示。40次施放探空仪达到的最大高度层为34～35 km。

图5 40次施放北斗-GPS双模式探空仪，不同窗口滑动平均风的各高度层系统误差:(a)南北风速分量,(b)东西风速分量

图6 40次施放北斗-GPS双模式探空仪，不同窗口滑动平均风的各高度层误差范围:(a)南北风速分量,(b)东西风速分量

从统计结果看，各窗口风速分量的系统误差较小，其数值在整个探测范围均在±0.3 m/s之间，具有可比较性。在整个探空高度内，10～20 s滑动平均风的风速分量误差范围较大；在0～32 km高度范围内，30～40 s滑动平均风的风速分量误差范围较小；在32 km高度以上接近球炸时，风速分量误差范围有较大幅度增加，此时窗口时间越长，风速分量误差范围越稳定、变化越小。

由于0～32 km高度范围内，30～40 s滑动平均风风速分量误差范围较接近，为获得更为清晰的图像，对30～40 s滑动平均风在0～32 km各高度层的误差范围上限和下限单独作图。从图7～8可以看出，34 s滑动平均风误差范围相对较小。

图7 图6中30～40 s滑动平均风在0～32 km各高度层误差范围上限:(a)南北风速分量,(b)东西风速分量

就上述分析，建议在0～32 km高度范围内取34 s作为滑动平均风窗口、在32 km高度以上取60 s 作为滑动平均风窗口。

图8 图6中30～40 s滑动平均风在0～32 km各高度层误差范围下限:(a)南北风速分量,(b)东西风速分量

4 计算方法应用

将本文建议的基于地心坐标系的平滑算法同文献[10]和[11]中提到的窗口为20 s和30 s、基于站心坐标系的算法进行比较。由图9三次试验结果可见，两种坐标系计算的结果总体一致，尤其是基于站心坐标系计算的30 s滑动平均风与本文算法计算结果差异不明显，基于站心坐标系计算的20 s滑动平均风依然存在着窗口短、曲线变化频率高的情况。

图9 两种坐标系下的北斗-GPS探空系统滑动平均风个例:(a1、a2)2013年12月15日,(b1、b2)2013年12月16日,(c1、c2)2013年12月17日(蓝色线为基于站心坐标系的20 s滑动平均风，绿色线为基于站心坐标系的30 s滑动平均风，红色线为基于地心坐标系、0～32 km窗口设置为34 s、32 km高度以上窗口设置为60 s的滑动平均风)

为分别计算结果间的细微差异，计算两种坐标系下滑动平均风风速分量在各高度层的95%置信概率误差范围，结果如图10所示。