基于地心坐標系的衛星導航測風平滑算法

李宛桐 姜明 史靜 楊榮康 黃威

(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所示。

圖10 40次施放北斗-GPS雙模式探空儀，兩種坐標系下的風速分量各高度層誤差范圍:(a)南北風速分量,(b)東西風速分量

從統計結果看，基于站心坐標系計算的20 s滑動平均風誤差范圍較大；在32 km高度以下，基于站心坐標系計算的30 s滑動平均風誤差范圍同本文算法誤差范圍相差不大，本文算法略優；在32 km高度以上，本文算法計算結果誤差范圍明顯小于其他，氣球升空全程誤差在±0.8～±2.0 m/s之間變化。

5 結論與討論

根據上述分析可以得到以下初步結論：

(1)根據95%置信概率的風速分量誤差范圍分析，建議衛星導航探空系統在0～32 km高度范圍內取34 s作為滑動平均窗口，在32 km高度以上取60 s作為滑動平均窗口。

(2)通過對兩種坐標系高空風平滑結果比較，基于站心坐標系計算的30 s滑動平均風與本文建議的基于地心坐標系的平滑算法計算結果基本相同。就風速分量誤差范圍而言，本文算法在32 km高度以下略優，在32 km高度以上減小明顯，升空全程誤差范圍更為穩定。