邵乙倫，趙娟，劉承禹，王巖，陳新

(1.上海交通大學(xué) 電子信息與電氣工程學(xué)院，上海 200240；2.96901部隊25分隊，北京 100094)

0 引言

無人機、智能車等移動載體利用各式傳感器收集環(huán)境反射信號來感知周邊環(huán)境.該環(huán)境信息可應(yīng)用于場景感知、避障、路徑規(guī)劃和定位增強等領(lǐng)域[1-2].在全球衛(wèi)星導(dǎo)航系統(tǒng)（GNSS）多徑反射信號利用方面，全球衛(wèi)星導(dǎo)航系統(tǒng)反射（GNSS-R）技術(shù)是一種新型的微波遙感技術(shù).該類技術(shù)的核心是通過觀測目標(biāo)GNSS反射信號特征參數(shù)來估算目標(biāo)參數(shù)，如測算海面高度、海面風(fēng)場、土壤濕度等[3-4].然而上述研究無法直接應(yīng)用于接收機周邊環(huán)境感知及導(dǎo)航定位.在城市峽谷環(huán)境下，GNSS信號被城市建筑等障礙物阻擋并產(chǎn)生大量多徑信號.目前，該類多徑信號主要被歸類為有害干擾信號并被GNSS接收機加以抑制或剔除[5-7].若能有效利用多徑信號中的空間環(huán)境信息，則可以在不增加額外開銷的情況下提高接收機的環(huán)境感知能力.

針對上述問題，本文提出了一種基于GNSS多徑信號的反射面參數(shù)估計算法.在GNSS多徑信號感知和特征參數(shù)提取方面，碼相位延遲幅度聯(lián)合跟蹤算法（CADLL）對接收機跟蹤環(huán)路的時域相關(guān)峰進(jìn)行檢測，在多徑信號存在時利用標(biāo)準(zhǔn)參考波形對發(fā)生畸變的時域相關(guān)峰進(jìn)行扣減，將多徑信號成分從時域相關(guān)峰中剝離[8].該算法在抑制多徑誤差的同時能夠?qū)崿F(xiàn)多徑信號特征參數(shù)的提取.在分析了GNSS多徑信號結(jié)構(gòu)及其幾何傳播過程后，給出了多徑信號傳播模型.由于相關(guān)模型含有大量的非線性運算，本文選用粒子濾波算法作為反射面定位算法.粒子濾波算法是一類貝葉斯估計方法，該算法利用蒙特卡洛法對復(fù)雜的統(tǒng)計積分公式進(jìn)行數(shù)值估計，能夠?qū)Ψ蔷€性及非高斯噪聲系統(tǒng)參數(shù)進(jìn)行估計[9-10].GNSS真實信號實驗結(jié)果表明，本文所提算法能夠有效估計GNSS多徑信號反射面位置參數(shù).

1 系統(tǒng)模型

1.1 GNSS多徑信號特征參數(shù)

城市峽谷等復(fù)雜場景下，GNSS信號經(jīng)一次或多次反射到達(dá)接收機天線，產(chǎn)生多徑效應(yīng).接收機接收到衛(wèi)星直達(dá)徑（LOS）信號和多徑信號.根據(jù)GNSS信號結(jié)構(gòu)，接收到的GNSS信號s(t)可表示為

式中：i=0表示LOS信號；A為LOS信號的幅值；D(t)為導(dǎo)航電文；αi為多徑信號相對LOS信號的衰減系數(shù)；C(t?τi)為擴頻碼；τi為多徑碼相位延遲；f0為載波頻率；fdi為載波多普勒頻率；φi為多徑信號相對LOS信號的相位變化.易知α0=1，τ0=0，φ0=0.多徑信號特征參數(shù)，即多徑衰減系數(shù)αi、碼相位延遲 τi、多普勒頻率fdi和相位偏移 φi，唯一確定了GNSS多徑信號.

1.2 GNSS多徑信號傳播模型

根據(jù)衛(wèi)星、接收機、反射面間的幾何關(guān)系，建立GNSS多徑信號傳播模型.該模型假設(shè)多徑信號經(jīng)反射面單次反射后到達(dá)接收機天線.以接收機初始時刻位置為原點建立東、北、天（E、N、U）地理坐標(biāo)系.下文所述的位置、速度參數(shù)均使用該地理坐標(biāo)系表示.傳播模型的幾何關(guān)系如圖1所示.衛(wèi)星位置為psat，衛(wèi)星速度為vsat，反射面鏡像衛(wèi)星位置pmirsat，鏡像衛(wèi)星速度為vmirsat，接收機位置為precv，接收機速度為vrecv，反射面法向量為np，反射面距離為d，反射面方位角φazi定義為反射面法向量在地理坐標(biāo)系水平面內(nèi)的投影向量與N方向間的夾角（北偏東為正），反射面俯仰角θele定義為反射面法向量與法向量水平面投影間的夾角（法向量指向天向時為正）.反射面法向量、反射面距離、反射面方位角、俯仰角、距離唯一確定了反射面位置參數(shù).

圖1 GNSS多徑信號傳播模型幾何關(guān)系

反射面法向量可以由反射面方位角、俯仰角表示：

反射面解析方程表達(dá)式如下：

鏡像衛(wèi)星的位置速度參數(shù)可以利用反射面解析方程計算得出，也可以采用類似于姿態(tài)求解過程中等效旋轉(zhuǎn)矢量的思路來求解[11].基于等效旋轉(zhuǎn)矢量方法，公式（4）、（5）給出了鏡像衛(wèi)星位置速度求解公式.

式中：I3為三維單位向量，(np×)為反射面法向量的反對稱矩陣.

由于多徑反射信號的幅值與反射面的介質(zhì)相關(guān)，因此很難建立幅值與反射面位置參數(shù)的數(shù)學(xué)模型.本文主要對多徑碼相位延遲τ和多普勒頻率漂移fdfading進(jìn)行建模分析.

多徑信號的碼相位延遲 τ是由LOS信號和多徑信號的傳播路徑長度差異造成，因此可以利用式（7）建立模型.

式中，“ ∥∥”為向量取模運算.

多徑多普勒頻率漂移fdfading反映了多徑信號多普勒頻率相對于LOS信號多普勒頻率的偏移，與衛(wèi)星和接收機的運動有關(guān)，可以建模為式（8）.

式中：λ為載波波長；“·”為向量點積運算.

2 粒子濾波多徑反射面參數(shù)估計

在接收機跟蹤環(huán)路，CADLL算法實現(xiàn)了GNSS多徑信號的感知和特征參數(shù)提取，為反射面參數(shù)估計提供了多徑碼相位延遲 τ和多普勒頻率漂移fdfading原始觀測量.接收機PVT解算模塊提供的接收機、衛(wèi)星的位置、速度信息為GNSS多徑信號傳播模型提供參數(shù)信息.由傳播模型可以看出，反射面位置參數(shù)與多徑信號觀測量之間存在大量的非線性運算，常規(guī)的卡爾曼濾波已無法滿足此類非線性觀測需求.因此，本文采用適用于非線性系統(tǒng)的粒子濾波算法來估計反射面位置參數(shù)信息，計算過程分為以下幾個步驟：

1)初始化

由于算法中所討論的GNSS多徑信號反射面主要指城市建筑外墻，因此反射面俯仰角初始化值的取值范圍限制在[?10,10]°；反射面方位角初始化值的取值范圍是[φazisat?90,φazisat+90]°，其中 φazisat為衛(wèi)星方位角；反射面距離初始化值的取值范圍是 [0,τ0]m，其中 τ0為初始多徑碼相位延遲.

2)系統(tǒng)狀態(tài)更新

本文假設(shè)反射面為單一平面且保持靜止，因此粒子更新方程如下式所示：

3)計算多徑觀測量預(yù)測值

根據(jù)每個粒子的方位角、俯仰角、距離值及接收機、衛(wèi)星的位置、速度信息，利用多徑信號傳播模型中式（7）、（8）計算預(yù)測多徑碼延遲和預(yù)測多普勒頻率漂移

4)粒子權(quán)重更新及歸一化

本文假設(shè)多徑觀測量服從獨立的二維高斯分布，則粒子權(quán)重更新過程如下式所示：

式中：σ1表示碼延遲觀測噪聲標(biāo)準(zhǔn)差；σ2表示多普勒頻率漂移觀測噪聲標(biāo)準(zhǔn)差.

若粒子俯仰角不滿足[?10,10]°取值范圍或粒子方位角不滿足式（13）所示的取值范圍，

式中，θelesat為衛(wèi)星俯仰角，則認(rèn)為該粒子值不滿足鏡面反射幾何關(guān)系并將該粒子的權(quán)重歸零.

粒子權(quán)重歸一化過程如下式所示：

5)粒子重采樣

根據(jù)更新后的粒子權(quán)重對粒子進(jìn)行重采樣操作，防止粒子退化.根據(jù)更新后的粒子權(quán)重值，計算權(quán)值累積分布函數(shù)（CDF），如下式所示：

6)計算反射面位置參數(shù)估計值

7)計算多徑信號反射點位置

多徑信號反射點坐標(biāo)(prpx,prpy,prpz)可通過求解如下方程組得到：

式中：t為輔助參數(shù)；(precvx,precvy,precvz)為接收機位置坐標(biāo)；nmp為鏡像衛(wèi)星單位視線矢量，坐標(biāo)為(nmpx,nmpy,nmpz)，計算公式如下：

3 實驗結(jié)果與分析

本文分兩次實驗對算法進(jìn)行測試.兩次實驗分別是上海交通大學(xué)微電子樓天臺的標(biāo)定實驗以及上海市陸家嘴市區(qū)環(huán)境實驗.實驗使用GNSS中頻采樣設(shè)備對GNSS原始信號采樣得到中頻原始數(shù)據(jù)，利用實驗室包含CADLL算法的軟件接收機對中頻原始數(shù)據(jù)進(jìn)行基帶信號處理和PVT解算得到多徑信號觀測數(shù)據(jù)及定位解算導(dǎo)航信息.將觀測數(shù)據(jù)代入粒子濾波反射面參數(shù)估計程序，計算得出反射面參數(shù)估計結(jié)果，其中反射面估計頻率為10Hz.

GNSS信號中頻采樣過程中盡量增加采樣帶寬以提高多徑特征參數(shù)提取精度，本文GNSS信號采樣帶寬為32MHz.

3.1 微電子樓天臺標(biāo)定實驗

上海交通大學(xué)微電子樓天臺北側(cè)和西側(cè)天空被建筑外墻遮擋，另外兩側(cè)天空則是開闊的.接收機天線置于天臺地面，方位角接近90°的低仰角衛(wèi)星信號經(jīng)西側(cè)建筑外墻反射產(chǎn)生多徑信號，GNSS天線同時接收LOS信號和多徑信號.微電子樓天臺標(biāo)定實驗中，GNSS信號及接收機參數(shù)如表1所示.

表1 天臺標(biāo)定實驗GNSS信號及接收機參數(shù)

使用高精度實時動態(tài)（RTK）接收機測量墻面標(biāo)定點坐標(biāo)，通過解析幾何計算可以得到地理坐標(biāo)系下的反射面參數(shù)標(biāo)定值.相對接收機位置處的地理坐標(biāo)系，天臺實驗的標(biāo)定結(jié)果為反射面方位角φazicali=69.8°，俯仰角θelecali=0°，距離dcali=19.6m.將多徑觀測數(shù)據(jù)代入粒子濾波反射面估計程序，反射面估計結(jié)果如圖2所示.

圖2 天臺標(biāo)定實驗反射面參數(shù)估計結(jié)果

根據(jù)反射面參數(shù)估計結(jié)果和反射面參數(shù)標(biāo)定值，評估反射面估計效果.反射面估計的絕對值誤差如圖3所示.

圖3 天臺標(biāo)定實驗反射面參數(shù)估計絕對值誤差

計算反射面估計的均方根誤差（RMSE），計算公式如式（23）所示，計算結(jié)果如表2所示.

式中：y為反射面參數(shù)估計值；ycali為反射面參數(shù)標(biāo)定值；N為實驗數(shù)據(jù)長度.

表2 天臺標(biāo)定實驗反射面位置估計RMSE

從實驗結(jié)果可以看出，反射面位置參數(shù)估計結(jié)果收斂于標(biāo)定值，估計RMSE較小，算法實現(xiàn)了對微電子樓天臺西側(cè)建筑外墻反射面位置參數(shù)的估計.

3.2 陸家嘴市區(qū)環(huán)境實驗

在陸家嘴市區(qū)環(huán)境進(jìn)行數(shù)據(jù)采集并進(jìn)行粒子濾波反射面參數(shù)估計，GNSS信號及接收機參數(shù)如表3所示，反射面參數(shù)估計結(jié)果如圖4所示.對反射點位置坐標(biāo)求平均值得到反射點位置估計，如圖5所示.

表3 陸家嘴市區(qū)環(huán)境GNSS信號及接收機參數(shù)

圖4 陸家嘴市區(qū)環(huán)境實驗反射面參數(shù)估計結(jié)果

圖5 陸家嘴市區(qū)環(huán)境實驗反射點位置估計結(jié)果

圖6為粒子濾波反射面參數(shù)估計算法多徑碼相位延遲原始觀測值.由圖6可知，碼相位延遲值小于80m.根據(jù)GNSS多徑反射信號傳播幾何關(guān)系，該多徑信號的反射面距離小于80m.由衛(wèi)星方位角為173.9°可知，反射面方位角的取值范圍為 [83.9,263.9]°.利用百度地圖接收機所在位置處的全景地圖，滿足上述方位角和距離要求的反射面僅有如圖5所示的東亞銀行金融大廈.算法實現(xiàn)了市區(qū)環(huán)境GNSS多徑信號反射面感知與參數(shù)估計.

圖6 陸家嘴市區(qū)環(huán)境實驗多徑碼相位延遲觀測值

由反射面估計結(jié)果可知反射面方位角在175°左右.由于本文反射面法向量定義為由反射面指向地理坐標(biāo)系坐標(biāo)原點，而衛(wèi)星方位角或百度全景地圖朝向的法向量定義為地理坐標(biāo)系中由接收機指向目標(biāo)，因此實驗估計的反射面在百度地圖全景功能中的朝向為355°.圖7為利用百度全景地圖得到的多徑反射場景示意圖，其中反射點所在建筑為東亞銀行金融大廈.

圖7 陸家嘴市區(qū)環(huán)境實驗多徑反射場景

4 結(jié)束語

本文對城市峽谷環(huán)境下GNSS信號中的LOS信號和多徑信號的傳播模型進(jìn)行建模分析，基于粒子濾波算法實現(xiàn)了多徑反射面的位置參數(shù)估計.該算法具有較高精度.靜態(tài)環(huán)境下，反射面方位角RMSE小于10°，俯仰角RMSE小于5°，距離估計RMSE小于10m.在不增加傳感器的情況下，本文所提算法提高了接收機的周邊環(huán)境感知能力，相關(guān)感知信息可應(yīng)用于場景感知、避障、路徑規(guī)劃和定位增強等領(lǐng)域.