孫宇 陳書(shū)馳 邵勝利 蔡俊武 何億強(qiáng) 王芳棟 史薈燕

(1 北京航天飛行控制中心，北京 100094； 2 中國(guó)人民解放軍32021部隊(duì)，北京 100094)

引言

臨近空間氣象火箭探空測(cè)風(fēng)模式為探空儀與火箭在上升段分離后，繼續(xù)慣性運(yùn)動(dòng)到頂點(diǎn)，在下落過(guò)程中，降落傘逐漸張開(kāi)，達(dá)到穩(wěn)定狀態(tài)后，在探空儀牽引下下落并隨風(fēng)飄移，降落傘成為高空風(fēng)場(chǎng)的示蹤物。根據(jù)探空儀記錄的位置變化，可以反演得到臨近空間20～60 km大氣風(fēng)場(chǎng)[1-7]。由此可知，相比于激光雷達(dá)或無(wú)線(xiàn)電雷達(dá)等遙感探測(cè)設(shè)備，火箭探空儀測(cè)風(fēng)屬于原位測(cè)量，其精度較高，對(duì)于工程應(yīng)用和比對(duì)試驗(yàn)具有重要意義[8-9]。

氣象火箭測(cè)風(fēng)資料處理與常規(guī)氣球探空不同。在50～60 km高度，由于空氣密度很低，對(duì)降落傘—探空儀系統(tǒng)的拖洩力較小，根據(jù)降落傘位置信息算出的水平運(yùn)動(dòng)速度并不等于該高度空氣塊的水平移動(dòng)速度，因此要對(duì)降落傘—探空儀系統(tǒng)的水平運(yùn)動(dòng)速度進(jìn)行修正，從而得到該高度空氣塊的水平運(yùn)動(dòng)速度(即該高度上的水平風(fēng)速)。由于風(fēng)速修正公式中要用到降落傘—探空儀系統(tǒng)的移動(dòng)速度和加速度，因此如何從所測(cè)的降落傘空間位置坐標(biāo)計(jì)算獲得這些速度和加速度？其誤差是多少，以及因此導(dǎo)致風(fēng)場(chǎng)反演結(jié)果的不確定度如何評(píng)估？這就是本文主要研究的問(wèn)題。

1 風(fēng)場(chǎng)反演數(shù)學(xué)模型

根據(jù)牛頓第二定律且忽略浮力和科氏力的影響，將降落傘—探空儀系統(tǒng)作為1個(gè)質(zhì)點(diǎn)，建立平衡方程如下：

(1)

將式(1)改寫(xiě)為三維分量表達(dá)式：

(2)

式中，ax、ay、az為降落傘-探空儀系統(tǒng)運(yùn)動(dòng)加速度的三維分量；vx、vy、vz為降落傘-探空儀系統(tǒng)運(yùn)動(dòng)速度的三維分量；Vx、Vy為風(fēng)速的南北、東西分量；Vz為垂直氣流速度；gx、gy、gz為重力加速度的三維分量。

由于重力加速度在X、Y方向上的分量很小，故忽略gx、gy，令gz=g；通常情況下垂直氣流速度Vz遠(yuǎn)小于降落傘-探空儀系統(tǒng)的落速，即vz-Vz≈vz。于是式(2)可簡(jiǎn)化為：

(3)

式(3)即為世界氣象組織推薦的風(fēng)速修正公式。

風(fēng)向G可由下式表達(dá)(具體方向由Vx、Vy的正負(fù)號(hào)確定)：

(4)

2 風(fēng)場(chǎng)反演不確定度

根據(jù)誤差理論[10]，若間接測(cè)量量y=f(x1、x2、…、xn)，則系統(tǒng)誤差Δy有如下近似表達(dá)式：

(5)

隨機(jī)誤差是用表征其取值分散程度的標(biāo)準(zhǔn)差來(lái)評(píng)定的，對(duì)于函數(shù)y的隨機(jī)誤差，可用函數(shù)的標(biāo)準(zhǔn)差來(lái)評(píng)定。若各測(cè)量值的隨機(jī)誤差相互獨(dú)立，則隨機(jī)誤差σy表達(dá)式如下：

(6)

式中，σx1，σx2，…，σxn為各變量的隨機(jī)誤差。

誤差合成采用方和根合成法，則測(cè)量結(jié)果的總標(biāo)準(zhǔn)差Δ總為：

(7)

風(fēng)場(chǎng)反演不確定度評(píng)估采用A類(lèi)評(píng)定法，其標(biāo)準(zhǔn)不確定度uc等同于由一系列觀測(cè)值獲得的標(biāo)準(zhǔn)差Δ總，即u=Δ總，測(cè)量結(jié)果不確定度Y表達(dá)式為：

Y=y±uc

(8)

2.1 風(fēng)速誤差

根據(jù)誤差理論，由式(3)，可得風(fēng)速系統(tǒng)誤差表達(dá)式如下：

(9)

由式(6)，可得風(fēng)速隨機(jī)誤差表達(dá)式如下：

(10)

隨機(jī)誤差與系統(tǒng)誤差合成采用方和根法。風(fēng)速誤差在X、Y方向上的分量為：

(11)

由誤差合成原理可得風(fēng)速誤差為：

(12)

2.2 風(fēng)向誤差

根據(jù)誤差理論，由式(4)，可得風(fēng)向誤差表達(dá)式為：

(13)

在實(shí)際數(shù)據(jù)處理時(shí)，由于速度和加速度分量是對(duì)離散位置擬合后獲得的，因此速度和加速度分量系統(tǒng)誤差和隨機(jī)誤差需由擬合過(guò)程得出。下面通過(guò)推導(dǎo)得到采用最小二乘法進(jìn)行多項(xiàng)式擬合的風(fēng)場(chǎng)系統(tǒng)誤差和隨機(jī)誤差表達(dá)式。

3 最小二乘法多項(xiàng)式擬合的系統(tǒng)誤差和隨機(jī)誤差

3.1 速度擬合的系統(tǒng)誤差

設(shè)探空儀的在下落過(guò)程中的真實(shí)位置可以通過(guò)一個(gè)k次多項(xiàng)式給出：

x(t)=a0+a1t+a2t2+…+aktk

(14)

(15)

式中，Pk(i)為正交多項(xiàng)式的k次項(xiàng)；P′k(i)為關(guān)于Pk(i)中對(duì)i的一階導(dǎo)數(shù)；Δt為數(shù)據(jù)點(diǎn)之間的時(shí)間間隔；N為用于擬合的數(shù)據(jù)點(diǎn)個(gè)數(shù)；k為擬合多項(xiàng)式次數(shù)。

根據(jù)正交多項(xiàng)式性質(zhì)可得:

(16)

(17)

下面對(duì)最小二乘法進(jìn)行1～4次多項(xiàng)式擬合的系統(tǒng)誤差公式進(jìn)行推導(dǎo)。

3.1.1 線(xiàn)性擬合

線(xiàn)性擬合條件下，擬合段中點(diǎn)的速度表達(dá)式為：

(18)

用式(14)代換為：

(19)

(20)

式(20)中大括號(hào)中的項(xiàng)可利用正交多項(xiàng)式的恒等式特點(diǎn)求解：

(21)

將式(21)代入式(20)可得：

(22)

由式(14)，擬合段的中點(diǎn)速度可表示為：

(23)

因此，對(duì)于線(xiàn)性擬合而言，速度的系統(tǒng)誤差(實(shí)際速度和擬合速度之間的差(式(22)-式(23))如下：

(24)

由于系數(shù)a0，a1…是與t0=0的坐標(biāo)相聯(lián)系，則式(24)可簡(jiǎn)化為：

(25)

3.1.2 三次多項(xiàng)式擬合

與線(xiàn)性擬合推導(dǎo)方法相同，同理可得到三次多項(xiàng)式速度擬合系統(tǒng)誤差公式：

Δvx=a5項(xiàng)+a7項(xiàng)+a9項(xiàng)+…

(26)

由式(25)、(26)可知，二次多項(xiàng)式擬合與線(xiàn)性擬合中點(diǎn)斜率相同(只有3次項(xiàng)系數(shù)以后的項(xiàng))，四次多項(xiàng)式擬合與三次多項(xiàng)式擬合中點(diǎn)斜率相同(只有5次項(xiàng)系數(shù)以后的項(xiàng))。因此，二次多項(xiàng)式擬合的速度系統(tǒng)誤差與線(xiàn)性擬合相同，四次多項(xiàng)式擬合的速度系統(tǒng)誤差與三次多項(xiàng)式擬合相同。

3.2 加速度擬合的系統(tǒng)誤差

與速度擬合求系統(tǒng)誤差推導(dǎo)方法相同，二次多項(xiàng)式擬合求二階導(dǎo)數(shù)獲取的加速度系統(tǒng)誤差表達(dá)式為：

(27)

如果坐標(biāo)的初始值t0=0，那么式(27)可簡(jiǎn)化為：

(28)

四次擬合k次多項(xiàng)式求二階導(dǎo)數(shù)的加速度系統(tǒng)誤差推導(dǎo)方法與之類(lèi)似，如果初始坐標(biāo)系是t0=0，那么系統(tǒng)誤差公式含有的項(xiàng)就只包括a6項(xiàng)，a8項(xiàng)……。也就是說(shuō)，二次多項(xiàng)式擬合求加速度時(shí)，其系統(tǒng)誤差與k次多項(xiàng)式a4項(xiàng)以后的偶數(shù)項(xiàng)相關(guān)；四次多項(xiàng)式擬合求加速度時(shí)，其系統(tǒng)誤差與k次多項(xiàng)式a6項(xiàng)以后的偶數(shù)項(xiàng)相關(guān)。

3.3 速度擬合的隨機(jī)誤差

(29)

(30)

整理可得：

(31)

(32)

(33)

因此，式(33)進(jìn)一步簡(jiǎn)化為：

(34)

從式(34)中誤差的系數(shù)項(xiàng)可以看出，二次多項(xiàng)式擬合與一次多項(xiàng)式擬合的速度隨機(jī)誤差相同，四次多項(xiàng)式擬合與三次多項(xiàng)式擬合的速度隨機(jī)誤差相同。

式(34)中右邊括弧里的項(xiàng)可表示為N的函數(shù)，由正交多項(xiàng)式遞推公式和βk+1聯(lián)合求解可得：

(35)

3.4 加速度擬合的隨機(jī)誤差

與速度擬合隨機(jī)誤差推導(dǎo)方法相同，同理可得k次多項(xiàng)式擬合段中點(diǎn)加速度表達(dá)式：

(36)

由式(36)中誤差的系數(shù)項(xiàng)可以看出，三次多項(xiàng)式擬合與二次多項(xiàng)式擬合的加速度隨機(jī)誤差相同，五次多項(xiàng)式擬合與四次多項(xiàng)式擬合的加速度隨機(jī)誤差相同。

式(36)括弧內(nèi)的前兩項(xiàng)可表示為N的函數(shù)：

(37)

通過(guò)以上推導(dǎo)，得到了速度和加速度1～4次多項(xiàng)式擬合的系統(tǒng)誤差和隨機(jī)誤差公式。結(jié)論總結(jié)如下：①對(duì)于速度擬合而言，二次多項(xiàng)式擬合的速度系統(tǒng)誤差、隨機(jī)誤差均與線(xiàn)性擬合相同；四次多項(xiàng)式擬合的速度系統(tǒng)誤差、隨機(jī)誤差與三次多項(xiàng)式擬合相同。②對(duì)于加速度擬合而言，二次多項(xiàng)式擬合的速度系統(tǒng)誤差、隨機(jī)誤差均與三次多項(xiàng)式擬合相同；四次多項(xiàng)式擬合的速度系統(tǒng)誤差、隨機(jī)誤差與五次多項(xiàng)式擬合相同。

下面以2004年11月16日11:59在酒泉地區(qū)獲取的1次探測(cè)數(shù)據(jù)為例，根據(jù)上述推導(dǎo)結(jié)論，計(jì)算風(fēng)場(chǎng)系統(tǒng)誤差、隨機(jī)誤差及其不確定度。

4 個(gè)例分析

從實(shí)際探測(cè)獲取的位移—時(shí)間曲線(xiàn)(圖1)看，位置滑動(dòng)擬合時(shí)(滑動(dòng)擬合點(diǎn)數(shù)根據(jù)分層高度選取)，采用3次多項(xiàng)式即可，即：

f(t)=a0+a1t+a2t2+a3t3

(38)

式中，f為緯向、經(jīng)向和垂向分量。

圖1 2004年11月16日11:59酒泉探空實(shí)測(cè)位移—時(shí)間曲線(xiàn): (a)緯向,(b)徑向,(c)垂向

4.1 雷達(dá)跟蹤定位誤差

雷達(dá)跟蹤體制應(yīng)答定位探空儀位置在站心坐標(biāo)系中坐標(biāo)計(jì)算公式為：

(39)

式中,x、y、z為探空儀位置點(diǎn)p在X、Y、Z軸上對(duì)應(yīng)的坐標(biāo)，單位:m；Rd為探空儀位置點(diǎn)p對(duì)應(yīng)的雷達(dá)跟蹤斜距，單位:m；α為探空儀位置點(diǎn)p對(duì)應(yīng)的雷達(dá)跟蹤仰角，單位:rad；δ為探空儀位置點(diǎn)p對(duì)應(yīng)的雷達(dá)跟蹤方位角，單位:rad；zh為探空儀的海拔高度，單位:m；h0為站心海拔高度，單位:m。

根據(jù)誤差傳遞原理，雷達(dá)跟蹤定位體制探空儀的定位誤差為：

(40)

式中，σx、σy、σz為探空儀定位誤差在站心坐標(biāo)系X、Y、Z方向上的分量，單位:m；σR為雷達(dá)測(cè)距誤差的標(biāo)準(zhǔn)差，單位:m；σγ為雷達(dá)測(cè)角誤差的標(biāo)準(zhǔn)差，單位:rad。

采用常規(guī)高空氣象702D雷達(dá)跟蹤定位，其測(cè)距精度約為20 m，測(cè)角精度約為0.1°。由式(40)計(jì)算可得探空儀位置分量的定位誤差曲線(xiàn)(圖2)?？梢钥吹?，x方向(緯向)全程定位誤差<60 m，y方向(經(jīng)向)全程定位誤差≤100 m，z方向(垂向)全程定位誤差≤80 m。

圖2 702D氣象雷達(dá)跟蹤定位誤差曲線(xiàn)

4.2 速度和加速度擬合

與常規(guī)高空氣象探測(cè)氣球勻速上升不同，火箭探空儀在下落過(guò)程中，落速逐漸減小，探測(cè)數(shù)據(jù)點(diǎn)分布呈上疏下密，因此需根據(jù)垂直分層要求，將探測(cè)高度每10 km間隔，劃分為50～60 km、40～50 km，30～40 km和20～30 km 4個(gè)高度區(qū)間，采用最小二乘法逐點(diǎn)對(duì)離散位置數(shù)據(jù)(x(t)、y(t)、z(t))進(jìn)行滑動(dòng)擬合。首先對(duì)位置數(shù)據(jù)進(jìn)行線(xiàn)性擬合，分別求一階導(dǎo)數(shù)獲得速度(vx(t)、vy(t)、vz(t))；然后對(duì)速度分量進(jìn)行2次多項(xiàng)式擬合，求一階導(dǎo)數(shù)獲得加速度(ax(t)、ay(t)、az(t))。

根據(jù)實(shí)測(cè)數(shù)據(jù)分析(表1)，50～60 km區(qū)間探空儀系統(tǒng)平均落速為87.6 m/s，擺動(dòng)周期約為7 s，采用7點(diǎn)擬合(采樣時(shí)間間隔Δt為1.05 s)，擬合高度分辨率約為610 m；40～50 km區(qū)間宜取11點(diǎn)平滑，擬合高度分辨率約為570 m；30～40 km區(qū)間宜取21點(diǎn)平滑，擬合高度分辨率約為410 m；30 km以下區(qū)間宜取35點(diǎn)平滑，擬合高度分辨率約為275 m。根據(jù)對(duì)擬合誤差的評(píng)估，通常采用線(xiàn)性、二次擬合平滑較為合適，即通過(guò)位置線(xiàn)性擬合求導(dǎo)獲得速度，再對(duì)速度進(jìn)行2次多項(xiàng)式擬合求導(dǎo)獲得加速度。根據(jù)雷達(dá)定位數(shù)據(jù)擬合后的速度分量見(jiàn)圖3。

表1 2004年11月16日11:59酒泉實(shí)測(cè)數(shù)據(jù)落速統(tǒng)計(jì)

圖3 速度分量擬合曲線(xiàn):(a)vx,(b)vy,(c)vz

4.3 擬合量系統(tǒng)誤差和隨機(jī)誤差計(jì)算

速度采用線(xiàn)性擬合時(shí)，由式(25)，其系統(tǒng)誤差為：

(41)

圖4 擬合的速度分量系統(tǒng)誤差:(a)Δvx,(b)Δvy,(c)Δvz

由于位置多項(xiàng)式方程的3次項(xiàng)(a3、b3、c3)未知，采用三次多項(xiàng)式最小二乘擬合系數(shù)代替。位置分量擬合系統(tǒng)誤差見(jiàn)圖4。可以看到，x方向(緯向)擬合系統(tǒng)誤差基本在±1 m/s以?xún)?nèi)；y方向(經(jīng)向)擬合系統(tǒng)誤差50～60 km基本在±1.5 m/s以?xún)?nèi)，50 km以下在±0.5 m/s以?xún)?nèi)；z方向(垂向)擬合系統(tǒng)誤差50～60 km約在±4 m/s以?xún)?nèi)，40～50 km約在±2 m/s以?xún)?nèi)，40 km以下約在±1 m/s以?xún)?nèi)。

加速度2次多項(xiàng)式擬合系統(tǒng)誤差分量中含有a4項(xiàng)、a6項(xiàng)等，故加速度擬合系統(tǒng)誤差為0。

由式(34)，速度分量隨機(jī)誤差為：

(42)

加速度分量隨機(jī)誤差為：

(43)

速度和加速度隨機(jī)誤差曲線(xiàn)如圖5、圖6所示，圖中速度和加速度誤差存在階躍變化的原因在于不同高度區(qū)間采取的擬合點(diǎn)數(shù)N不同。從圖5可知，x方向(緯向)速度擬合隨機(jī)誤差在±1 m/s以?xún)?nèi)；y方向(經(jīng)向)速度擬合隨機(jī)誤差基本在±2 m/s以?xún)?nèi)；z方向(垂向)速度擬合隨機(jī)誤差在±1.5 m/s以?xún)?nèi)。

速度擬合合成誤差如圖7所示，x方向(緯向)為±1.7 m/s，y方向(經(jīng)向)為±3.3 m/s。

圖5 擬合的速度分量隨機(jī)誤差:(a)σvx,(b)σvy,(c)σvz

圖6 擬合的加速度分量隨機(jī)誤差:(a)σax,(b)σay,(c)σaz

圖7 擬合的速度分量合成誤差:(a)Vx總,(b)Vy總

圖8為修正后經(jīng)緯向風(fēng)速曲線(xiàn)，可以看到，在50 km以上經(jīng)緯向風(fēng)修正較大，50 km以下修正較小。

圖8 修正后經(jīng)向風(fēng)(a)和緯向風(fēng)(b)風(fēng)速隨高度變化

圖9為臨近空間20～60 km風(fēng)場(chǎng)反演結(jié)果及不確定度曲線(xiàn)。可以看到，風(fēng)速不確定度隨高度降低而減小。表2、3、4為風(fēng)速和風(fēng)向反演結(jié)果及不確定度?？梢钥吹?，50～60 km為2.8～3.5 m/s,50 km以下在1 m/s以?xún)?nèi)。其中，在58.3 km高度風(fēng)向不確定度異常增大到44.9°,31 km高度異常增大到18.1°，除風(fēng)向0°漸變段(56～59 km，30～32 km)外，其他高度時(shí)基本在10°以?xún)?nèi)。

值得注意的是，風(fēng)向在0°附近擺動(dòng)時(shí)(此時(shí)風(fēng)速也接近于0 m/s)，會(huì)導(dǎo)致風(fēng)向反演不確定度增大(圖9框選部分)。

圖9 反演風(fēng)場(chǎng)及不確定度: (a)合成風(fēng)風(fēng)速,(b)風(fēng)速不確定度,(c)風(fēng)向,(d)風(fēng)向不確定度

表2 風(fēng)場(chǎng)反演不確定度

表3 50～60 km風(fēng)向反演不確定度

表4 30～35 km風(fēng)向反演不確定度

5 結(jié)論

氣象探空火箭探測(cè)臨近空間風(fēng)場(chǎng)是建立在探空儀位置獲取的基礎(chǔ)上，其反演誤差來(lái)自雷達(dá)定位誤差和擬合誤差。雷達(dá)定位誤差與自身定位精度、實(shí)際探測(cè)時(shí)仰角和斜距相關(guān)。因此，風(fēng)場(chǎng)反演不確定的評(píng)估是基于多項(xiàng)式擬合和定位誤差的分析，本文通過(guò)對(duì)多項(xiàng)式最小二乘擬合的誤差分析研究，分別推導(dǎo)給出了系統(tǒng)誤差和隨機(jī)誤差公式。采用一次實(shí)測(cè)數(shù)據(jù)計(jì)算表明，風(fēng)速反演不確定度隨高度降低而減小，50～60 km為2.8～3.5 m/s,50 km以下在1 m/s 以?xún)?nèi)；風(fēng)向在0°附近擺動(dòng)時(shí)(此時(shí)風(fēng)速也接近于0 m/s)，會(huì)導(dǎo)致風(fēng)向反演不確定度異常增大，在其他高度時(shí)風(fēng)向反演不確定度基本在10°以?xún)?nèi)。若采用更高精度的雷達(dá)，在高空大斜距條件下風(fēng)場(chǎng)反演精度提升有限，而且會(huì)增加火箭探空儀結(jié)構(gòu)設(shè)計(jì)的復(fù)雜性。若火箭探空儀采用導(dǎo)航衛(wèi)星定位體制[12]，獲取的探空儀位置和實(shí)時(shí)速度在精度上會(huì)有量級(jí)上的提高，則可以極大地提升風(fēng)場(chǎng)反演精度。