基于地心坐標(biāo)系的衛(wèi)星導(dǎo)航測(cè)風(fēng)平滑算法

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

(1 天津市氣象探測(cè)中心，天津 300061；2 中國(guó)氣象局氣象探測(cè)中心，北京 100081；3 中國(guó)洛陽電子裝備試驗(yàn)中心，濟(jì)源 459000)

引言

20世紀(jì)80年代以來，出現(xiàn)了奧米伽、羅蘭-C以及全球定位系統(tǒng)(Global Positioning System, GPS)、北斗等導(dǎo)航測(cè)風(fēng)系統(tǒng)[1]，隨著北斗衛(wèi)星導(dǎo)航系統(tǒng)的發(fā)展以及考慮探空經(jīng)濟(jì)成本和加密觀測(cè)的進(jìn)一步需求，我國(guó)探空業(yè)務(wù)將從二次雷達(dá)測(cè)風(fēng)逐步過渡到衛(wèi)星導(dǎo)航測(cè)風(fēng)[2]。與此同時(shí)，國(guó)產(chǎn)衛(wèi)星導(dǎo)航探空系統(tǒng)的數(shù)據(jù)采集率、準(zhǔn)確度以及系統(tǒng)的自動(dòng)化程度有很大提升，其數(shù)據(jù)采集率已由原59-701探空系統(tǒng)的每分鐘1組提高到每秒鐘1組[3-4]，而現(xiàn)行《常規(guī)高空氣象觀測(cè)業(yè)務(wù)規(guī)范》(以下簡(jiǎn)稱《規(guī)范》)規(guī)定的業(yè)務(wù)應(yīng)用數(shù)據(jù)和文件計(jì)算數(shù)學(xué)模型還停留在以分鐘間隔處理的基礎(chǔ)上[5]。分鐘風(fēng)的計(jì)算方式浪費(fèi)了大量數(shù)據(jù)資源[6]，不能完全發(fā)揮衛(wèi)星導(dǎo)航探空系統(tǒng)的效益，但是利用秒間隔定位數(shù)據(jù)直接計(jì)算的原始高空風(fēng)脈動(dòng)太大，不能用于生成業(yè)務(wù)應(yīng)用文件。

關(guān)于原始高空風(fēng)數(shù)據(jù)的處理，在目前業(yè)務(wù)應(yīng)用的L波段探空系統(tǒng)中，已有很多專家論證采用矢量滑動(dòng)平均法，并對(duì)滑動(dòng)平均窗口提出了建議：王緬等[7]建議采用窗口為1 min的滑動(dòng)平均方式或者前20 min采用30 s，以后采用1 min的分段滑動(dòng)平均方法；梁建平等[8]提出使用30～45 s的滑動(dòng)平均窗口；陳磊等[9]認(rèn)為在全程使用50～60 s窗口或前50 min使用30～40 s窗口，50 min以后使用50～60 s窗口結(jié)果較好。在衛(wèi)星導(dǎo)航測(cè)風(fēng)相關(guān)研究中，也有不少學(xué)者使用矢量滑動(dòng)平均法處理原始高空風(fēng)數(shù)據(jù)：張國(guó)舫等[10]在比較GPS探空儀多普勒效應(yīng)方法和探空氣球定位方法計(jì)算高空風(fēng)時(shí)，采用20 s作為矢量滑動(dòng)平均窗口；姚雯等[11]在評(píng)估國(guó)產(chǎn)GPS探空儀秒級(jí)探空數(shù)據(jù)隨機(jī)誤差時(shí)，對(duì)秒間隔高空風(fēng)數(shù)據(jù)進(jìn)行了30點(diǎn)和15點(diǎn)滑動(dòng)平均，并證明平滑程度的差異對(duì)間接計(jì)算風(fēng)的標(biāo)準(zhǔn)差影響較大。

上述研究在利用秒間隔定位數(shù)據(jù)計(jì)算原始高空風(fēng)時(shí)，均參照《規(guī)范》規(guī)定的量得風(fēng)層計(jì)算方法，即基于站心坐標(biāo)系進(jìn)行計(jì)算：首先將氣球的大地坐標(biāo)垂直投影到站心坐標(biāo)平面上，然后用該平面上兩相鄰?fù)队包c(diǎn)的坐標(biāo)數(shù)據(jù)計(jì)算風(fēng)向風(fēng)速，沒有考慮地球曲率的影響。風(fēng)是平行于地面的空氣流動(dòng)，氣球測(cè)風(fēng)是將氣球看作隨氣流移動(dòng)的質(zhì)點(diǎn)，氣球在升空過程中同時(shí)受上升力和地球引力的影響，基于兩種力的作用，其實(shí)際的飛行軌跡，是在上升和隨空氣的水平運(yùn)動(dòng)向遠(yuǎn)方漂移同時(shí)，對(duì)于地心坐標(biāo)系不斷下降，始終沿地球曲面飛行，即氣球的運(yùn)動(dòng)是基于地心坐標(biāo)系的[12]。此外，根據(jù)天氣學(xué)的要求，規(guī)定高度風(fēng)、規(guī)定標(biāo)準(zhǔn)氣壓層風(fēng)和最大風(fēng)層等需要以量得風(fēng)層為基礎(chǔ)內(nèi)插計(jì)算的高空風(fēng)資料，其高度是以垂直于地球表面的位勢(shì)高度表示的[13]，即業(yè)務(wù)上分析量得風(fēng)層時(shí)所使用的站心坐標(biāo)系與使用高空風(fēng)資料的坐標(biāo)系并不一致，造成這種不一致的原因與59-701系統(tǒng)下風(fēng)向風(fēng)速以人工計(jì)算為主有關(guān)[14]。如今計(jì)算機(jī)廣泛應(yīng)用，衛(wèi)星導(dǎo)航探空系統(tǒng)提供的氣球定位又是大地坐標(biāo)，應(yīng)在地心坐標(biāo)系中計(jì)算高空風(fēng)[15]。

本文采用2013年12月中國(guó)氣象局在陽江探空站進(jìn)行北斗-GPS雙模式探空儀試驗(yàn)的實(shí)際施放數(shù)據(jù)，分析基于地心坐標(biāo)系的衛(wèi)星導(dǎo)航測(cè)風(fēng)矢量滑動(dòng)平均窗口選擇，為探空業(yè)務(wù)高空風(fēng)計(jì)算方法提供參考。

1 高空風(fēng)數(shù)據(jù)處理

1.1 坐標(biāo)系統(tǒng)

在衛(wèi)星導(dǎo)航系統(tǒng)中，坐標(biāo)系是描述衛(wèi)星運(yùn)動(dòng)、處理觀測(cè)數(shù)據(jù)和求解接收機(jī)位置的數(shù)學(xué)和物理計(jì)算基礎(chǔ)[16]。利用坐標(biāo)系可以方便地為衛(wèi)星導(dǎo)航建立數(shù)學(xué)公式，當(dāng)參考坐標(biāo)系確立后，衛(wèi)星和接收機(jī)的狀態(tài)就很容易用數(shù)學(xué)模型表示出來[17]。根據(jù)應(yīng)用場(chǎng)合的不同，選用的坐標(biāo)系也不相同。

1.1.1 地心大地坐標(biāo)系

目前，GPS所采用的坐標(biāo)系是世界大地坐標(biāo)系(World Geodetic System of 1984, WGS-84)，GPS所有的星歷參數(shù)就是基于此坐標(biāo)系[18]。WGS-84坐標(biāo)系的原點(diǎn)位于地球質(zhì)心，其Z軸指向國(guó)際時(shí)間局(BIH)1984年00:00定義的BIH1984.0協(xié)議地球極(Conventional Terrestrial Pole, CTP)方向，X軸指向BIH1984.0的零子午面與CTP赤道的交點(diǎn)，X軸、Y軸、Z軸構(gòu)成右手坐標(biāo)系，如圖1所示。

圖1 世界大地坐標(biāo)系(WGS-84)

北斗所使用的坐標(biāo)系是2000國(guó)家大地坐標(biāo)系(China Geodetic Coordinate System 2000, CGCS2000)，自2008年7月1日開始啟用，是全球地心坐標(biāo)系在我國(guó)的具體體現(xiàn)。CGCS坐標(biāo)系的原點(diǎn)為地球質(zhì)心，初始定向由1984.0的BIH定向給定，CGCS2000的參考?xì)v元為2000.0[19]。

CGCS2000與WGS-84坐標(biāo)系在原點(diǎn)、尺度、定向及定向的定義基本相同，參考橢球非常相近，只有扁率f有微小的差異，鑒于在坐標(biāo)系定義和實(shí)現(xiàn)上的比較，可以認(rèn)為CGCS2000與WGS-84是相容的；在坐標(biāo)系的實(shí)現(xiàn)精度范圍內(nèi)，CGCS2000與WGS-84是一致的[20-22]。

在地心大地坐標(biāo)系中，點(diǎn)用(L,B,H)來表示。其中，L表示經(jīng)度，B表示緯度，H表示高程。

1.1.2 地心地固坐標(biāo)系

地心地固坐標(biāo)系(Earth-Centered Earth-Fixed, ECEF)簡(jiǎn)稱地心坐標(biāo)系，是以地球質(zhì)心為原點(diǎn)的笛卡兒坐標(biāo)系。X軸指向0°大地子午線與赤道的交點(diǎn)，Y軸指向東經(jīng)90°大地子午線與赤道的交點(diǎn)，Z軸指向協(xié)議地球北極，組成右手系直角坐標(biāo)，如圖2所示。

圖2 地心坐標(biāo)系(ECEF)

在用戶接收機(jī)中，因?yàn)槔昧薊CEF坐標(biāo)系表示衛(wèi)星位置(作為定位參考點(diǎn))，所以計(jì)算得到的用戶位置一般也表征在ECEF坐標(biāo)系中。

1.1.3 坐標(biāo)轉(zhuǎn)換

衛(wèi)星導(dǎo)航探空系統(tǒng)測(cè)風(fēng)時(shí)，設(shè)某一時(shí)刻氣球的定位數(shù)據(jù)Pi(Li,Bi,Hi)。在ECEF坐標(biāo)系中，Li為ZOX平面與ZOPi平面的夾角，自ZOX平面起算右旋為正；Bi為過Pi點(diǎn)的橢球面法線與XOY平面的夾角，自XOY平面向z軸方向量取為正；Hi為過Pi點(diǎn)的橢球面法線自橢球面至氣球的距離，以遠(yuǎn)離橢球面中心方向?yàn)檎?，此法線和z軸的交點(diǎn)與地心通常并不重合。據(jù)此，氣球某一時(shí)刻在ECEF坐標(biāo)系中的位置Pi(Xi,Yi,Zi)可表示如下：

(1)

式中，Ni為卯酉圈曲率半徑；e為地球橢球第一偏心率，用下式計(jì)算：

(2)

(3)

其中，a=6378137 m為地球的長(zhǎng)半軸長(zhǎng)度，b=6356752.3142 m為地球的短半軸長(zhǎng)度。

1.2 基于地心坐標(biāo)系計(jì)算原始高空風(fēng)

參考李浩等的研究成果[15]，基于地心坐標(biāo)系的衛(wèi)星導(dǎo)航探空系統(tǒng)原始高空風(fēng)計(jì)算方法如下：

將相鄰時(shí)刻氣球位置記為Pi-1(Li-1,Bi-1,Hi-1)和Pi(Li,Bi,Hi)，選取高度為(Hi-1+Hi)/2的橢球面作為計(jì)算高空風(fēng)的水平面，把該水平面叫做投影橢球面，如圖3所示。圖3中已標(biāo)出投影橢球面上的經(jīng)線和緯線。

圖3 投影橢球面

把Pi-1、Pi分別沿卯酉圈曲率半徑方向向上、向下投影到該投影橢球面上，形成的投影點(diǎn)分別記作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經(jīng)線上，投影點(diǎn)記作N；把Pro,i-1沿著Li-1經(jīng)線投影到Bi緯線上，投影點(diǎn)記作M。顯然，N、M的坐標(biāo)分別是(Li,Bi-1,(Hi-1+Hi)/2)、(Li-1,Bi,(Hi-1+Hi)/2)。

(4)

注意到dPN≠dMP，定義(dPN+dMP)/2為氣球在緯圈(即東西方向)上移動(dòng)的水平位移，顯然

(5)

將氣球經(jīng)向運(yùn)動(dòng)速度記作vx，即南北風(fēng)速分量；緯向運(yùn)動(dòng)速度記作vy，即東西風(fēng)速分量。有

(6)

對(duì)應(yīng)風(fēng)速V和風(fēng)向D：

(7)

(8)

1.3 原始高空風(fēng)數(shù)據(jù)平滑處理

本文采用矢量滑動(dòng)平均法對(duì)衛(wèi)星導(dǎo)航探空系統(tǒng)原始高空風(fēng)數(shù)據(jù)進(jìn)行處理。矢量滑動(dòng)平均法要求設(shè)置平均時(shí)段(窗口)和步進(jìn)時(shí)間，步進(jìn)時(shí)間通常以原始數(shù)據(jù)的時(shí)間間隔為準(zhǔn)。利用公式(9)對(duì)每一平均時(shí)段風(fēng)速分量進(jìn)行平滑處理，所得平均值賦給該平均時(shí)段的中間時(shí)。

(9)

式中，n為平均時(shí)段內(nèi)風(fēng)速分量的個(gè)數(shù)。

計(jì)算完成第一個(gè)時(shí)段的兩個(gè)風(fēng)速分量以后，在原來時(shí)段的前面去掉一個(gè)步進(jìn)時(shí)間的數(shù)據(jù)，在其后面增加一個(gè)步進(jìn)時(shí)間的數(shù)據(jù)，再進(jìn)行同樣的計(jì)算，直至最后一組數(shù)據(jù)。這樣就得到規(guī)定窗口的兩個(gè)風(fēng)速分量的滑動(dòng)平均值。由于計(jì)算得到的風(fēng)速分量起始時(shí)間為第一個(gè)平均時(shí)段的中間時(shí)刻，在這一時(shí)刻之前的數(shù)據(jù)可以通過與地面風(fēng)速分量數(shù)據(jù)內(nèi)插獲得，最終得到從地面至球炸的連續(xù)數(shù)據(jù)。

2 實(shí)例數(shù)據(jù)及分析方法

本文采用2013年12月中國(guó)氣象局在陽江探空站進(jìn)行北斗-GPS雙模式探空儀試驗(yàn)的實(shí)際施放數(shù)據(jù)。試驗(yàn)采取同球施放方法，即在一個(gè)探空氣球上同時(shí)懸掛一個(gè)RS92探空儀與兩個(gè)北斗-GPS雙模式探空儀，北斗-GPS雙模式探空儀分別以北斗-GPS綜合模式、單GPS模式和單北斗模式解算秒間隔定位數(shù)據(jù)。

3 矢量滑動(dòng)平均窗口選擇

3.1 窗口對(duì)測(cè)風(fēng)結(jié)果的影響

根據(jù)上述公式計(jì)算全部探空數(shù)據(jù)10 s、20 s、30 s、40 s、50 s、60 s滑動(dòng)平均風(fēng)，并與RS92探空系統(tǒng)給出的風(fēng)速風(fēng)向相比較，其結(jié)果如下。

圖4是比對(duì)結(jié)果的典型例子，圖中藍(lán)色曲線為北斗-GPS探空系統(tǒng)滑動(dòng)平均風(fēng)速風(fēng)向值。從圖中可以看出，平滑后的風(fēng)速、風(fēng)向曲線與RS92探空系統(tǒng)趨勢(shì)基本一致。和RS92測(cè)風(fēng)數(shù)據(jù)相比，10 s滑動(dòng)平均風(fēng)的曲線波動(dòng)最為明顯，隨著窗口的延長(zhǎng)，曲線變化幅度衰減得越大，變化頻率越低。由于窗口較短、波動(dòng)較大，10 s和20 s滑動(dòng)平均風(fēng)的風(fēng)向在升空擾動(dòng)明顯處，如遇過北情況時(shí)偏差較大。當(dāng)窗口為40 s及以上時(shí)，滑動(dòng)平均風(fēng)風(fēng)速、風(fēng)向曲線趨于平滑，擺動(dòng)幅度逐漸減小，同時(shí)產(chǎn)生相位滯后，其滯后時(shí)間隨窗口的延長(zhǎng)而增加。

圖4 2013年12月18日不同滑動(dòng)平均窗口北斗-GPS探空系統(tǒng)滑動(dòng)平均(藍(lán)色線)與RS92探空系統(tǒng)探測(cè)(紅色線)風(fēng)向風(fēng)速變化:(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 誤差分析

統(tǒng)計(jì)窗口分別為10 s、12 s、14 s、…、56 s、58 s、60 s的滑動(dòng)平均風(fēng)在每千米高度層上風(fēng)速分量的誤差，繪制風(fēng)速分量系統(tǒng)誤差變化曲線和95%置信概率誤差范圍變化圖，如圖5和圖6所示。40次施放探空儀達(dá)到的最大高度層為34～35 km。

圖5 40次施放北斗-GPS雙模式探空儀，不同窗口滑動(dòng)平均風(fēng)的各高度層系統(tǒng)誤差:(a)南北風(fēng)速分量,(b)東西風(fēng)速分量

圖6 40次施放北斗-GPS雙模式探空儀，不同窗口滑動(dòng)平均風(fēng)的各高度層誤差范圍:(a)南北風(fēng)速分量,(b)東西風(fēng)速分量

從統(tǒng)計(jì)結(jié)果看，各窗口風(fēng)速分量的系統(tǒng)誤差較小，其數(shù)值在整個(gè)探測(cè)范圍均在±0.3 m/s之間，具有可比較性。在整個(gè)探空高度內(nèi)，10～20 s滑動(dòng)平均風(fēng)的風(fēng)速分量誤差范圍較大；在0～32 km高度范圍內(nèi)，30～40 s滑動(dòng)平均風(fēng)的風(fēng)速分量誤差范圍較?。辉?2 km高度以上接近球炸時(shí)，風(fēng)速分量誤差范圍有較大幅度增加，此時(shí)窗口時(shí)間越長(zhǎng)，風(fēng)速分量誤差范圍越穩(wěn)定、變化越小。

由于0～32 km高度范圍內(nèi)，30～40 s滑動(dòng)平均風(fēng)風(fēng)速分量誤差范圍較接近，為獲得更為清晰的圖像，對(duì)30～40 s滑動(dòng)平均風(fēng)在0～32 km各高度層的誤差范圍上限和下限單獨(dú)作圖。從圖7～8可以看出，34 s滑動(dòng)平均風(fēng)誤差范圍相對(duì)較小。

圖7 圖6中30～40 s滑動(dòng)平均風(fēng)在0～32 km各高度層誤差范圍上限:(a)南北風(fēng)速分量,(b)東西風(fēng)速分量

就上述分析，建議在0～32 km高度范圍內(nèi)取34 s作為滑動(dòng)平均風(fēng)窗口、在32 km高度以上取60 s 作為滑動(dòng)平均風(fēng)窗口。

圖8 圖6中30～40 s滑動(dòng)平均風(fēng)在0～32 km各高度層誤差范圍下限:(a)南北風(fēng)速分量,(b)東西風(fēng)速分量

4 計(jì)算方法應(yīng)用

將本文建議的基于地心坐標(biāo)系的平滑算法同文獻(xiàn)[10]和[11]中提到的窗口為20 s和30 s、基于站心坐標(biāo)系的算法進(jìn)行比較。由圖9三次試驗(yàn)結(jié)果可見，兩種坐標(biāo)系計(jì)算的結(jié)果總體一致，尤其是基于站心坐標(biāo)系計(jì)算的30 s滑動(dòng)平均風(fēng)與本文算法計(jì)算結(jié)果差異不明顯，基于站心坐標(biāo)系計(jì)算的20 s滑動(dòng)平均風(fēng)依然存在著窗口短、曲線變化頻率高的情況。

圖9 兩種坐標(biāo)系下的北斗-GPS探空系統(tǒng)滑動(dòng)平均風(fēng)個(gè)例:(a1、a2)2013年12月15日,(b1、b2)2013年12月16日,(c1、c2)2013年12月17日(藍(lán)色線為基于站心坐標(biāo)系的20 s滑動(dòng)平均風(fēng)，綠色線為基于站心坐標(biāo)系的30 s滑動(dòng)平均風(fēng)，紅色線為基于地心坐標(biāo)系、0～32 km窗口設(shè)置為34 s、32 km高度以上窗口設(shè)置為60 s的滑動(dòng)平均風(fēng))

為分別計(jì)算結(jié)果間的細(xì)微差異，計(jì)算兩種坐標(biāo)系下滑動(dòng)平均風(fēng)風(fēng)速分量在各高度層的95%置信概率誤差范圍，結(jié)果如圖10所示。

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

從統(tǒng)計(jì)結(jié)果看，基于站心坐標(biāo)系計(jì)算的20 s滑動(dòng)平均風(fēng)誤差范圍較大；在32 km高度以下，基于站心坐標(biāo)系計(jì)算的30 s滑動(dòng)平均風(fēng)誤差范圍同本文算法誤差范圍相差不大，本文算法略優(yōu)；在32 km高度以上，本文算法計(jì)算結(jié)果誤差范圍明顯小于其他，氣球升空全程誤差在±0.8～±2.0 m/s之間變化。

5 結(jié)論與討論

根據(jù)上述分析可以得到以下初步結(jié)論：

(1)根據(jù)95%置信概率的風(fēng)速分量誤差范圍分析，建議衛(wèi)星導(dǎo)航探空系統(tǒng)在0～32 km高度范圍內(nèi)取34 s作為滑動(dòng)平均窗口，在32 km高度以上取60 s作為滑動(dòng)平均窗口。

(2)通過對(duì)兩種坐標(biāo)系高空風(fēng)平滑結(jié)果比較，基于站心坐標(biāo)系計(jì)算的30 s滑動(dòng)平均風(fēng)與本文建議的基于地心坐標(biāo)系的平滑算法計(jì)算結(jié)果基本相同。就風(fēng)速分量誤差范圍而言，本文算法在32 km高度以下略優(yōu)，在32 km高度以上減小明顯，升空全程誤差范圍更為穩(wěn)定。