王式太 姜新偉 殷 敏 王文貫 蘭小艷 楊可心

1 桂林理工大學(xué)測(cè)繪地理信息學(xué)院，桂林市雁山街319號(hào)，541006 2 廣西空間信息與測(cè)繪重點(diǎn)實(shí)驗(yàn)室，桂林市雁山街319號(hào)，541006 3 廣西建設(shè)職業(yè)技術(shù)學(xué)院土木工程學(xué)院，南寧市羅文大道33號(hào)，530007

土壤濕度又稱土壤含水率，是全球碳水循環(huán)中的重要環(huán)節(jié)，對(duì)于水文和氣象研究、農(nóng)業(yè)發(fā)展具有重要價(jià)值[1]。傳統(tǒng)的土壤濕度測(cè)量方法(如中子儀法、烘干稱重法、頻域反射儀法、時(shí)域反射儀法等)操作過程較為復(fù)雜、探測(cè)范圍小、時(shí)空分辨率較低[2]。而全球?qū)Ш叫l(wèi)星系統(tǒng)反射(GNSS-R)技術(shù)可利用GNSS反射信號(hào)估計(jì)反射面相關(guān)參數(shù)[3]。Martin-Neira[4]通過GPS反射信號(hào)首次實(shí)現(xiàn)海洋高度測(cè)量；Larson等[5]根據(jù)海水潮汐變化率提出動(dòng)態(tài)改正方法，可提高海洋高度變化情況下的測(cè)量精度。在土壤濕度方面，Larson等[6]認(rèn)為土壤濕度與GNSS接收機(jī)記錄的信噪比(SNR)震蕩幅度有關(guān)，并表明該特征是由衛(wèi)星直射信號(hào)與反射信號(hào)干涉形成。在此基礎(chǔ)上，孫波等[7]考慮信噪比序列的相位、頻率和振幅，通過GA-SVM構(gòu)建該特征與土壤濕度的非線性回歸模型；荊麗麗等[8]提出使用熵值法對(duì)2個(gè)頻點(diǎn)的數(shù)據(jù)進(jìn)行融合，從而提高土壤濕度的反演精度；文獻(xiàn)[9-10]通過小波變換方法去除信噪比序列中的趨勢(shì)項(xiàng)和高頻誤差，從而提高反演精度。

以上對(duì)于土壤濕度的反演研究均采用標(biāo)準(zhǔn)余弦函數(shù)擬合信噪比序列求解特征參數(shù)，未考慮信噪比振幅隨時(shí)間和高度角的變化情況。因此，本文選擇顧及阻尼因子的信噪比干涉模型，考慮信噪比振幅隨時(shí)間和衛(wèi)星仰角的變化，并通過遺傳算法(GA)輔助非線性最小二乘方法(NLS)求解干涉模型中的反射信號(hào)參數(shù)，提高特征參數(shù)的求解精度。通過仿真實(shí)驗(yàn)分析該方法的有效性，并基于60 d的實(shí)驗(yàn)數(shù)據(jù)對(duì)2種不同計(jì)算方案進(jìn)行對(duì)比分析。

1 GNSS-IR基本原理

信噪比(SNR)通常用于衡量信號(hào)的質(zhì)量，定義為信號(hào)接收功率Pr與噪聲功率N0之比，即

(1)

由于接收機(jī)實(shí)際接收到的為干涉信號(hào)，因此可將信噪比公式近似表達(dá)為：

(2)

式中，Ad為直射信號(hào)振幅；Am為反射信號(hào)振幅；ψ為直射信號(hào)與反射信號(hào)的相位差，由直射信號(hào)與反射信號(hào)傳播路程存在差異(光程差)所致，可表示為：

(3)

式中，λ為衛(wèi)星載波波長(zhǎng)，θ為天線處衛(wèi)星高度角[11]。將ψ代入式(2)，則SNR可表示為：

(4)

式中，φ為相對(duì)延遲相位。

接收機(jī)天線增益會(huì)隨衛(wèi)星高度角的增大而增大，接收的直射信號(hào)強(qiáng)度會(huì)隨衛(wèi)星高度角的增大而提高，而反射信號(hào)強(qiáng)度會(huì)隨衛(wèi)星高度角的增大而下降，導(dǎo)致整個(gè)信噪比序列呈現(xiàn)二次曲線變化趨勢(shì)。通常使用低階多項(xiàng)式擬合來表示該趨勢(shì)，通過將原始信噪比序列減去趨勢(shì)項(xiàng)來構(gòu)建多路徑信噪比SNRMP。低高度角范圍下多路徑信噪比可表示為余弦函數(shù)形式：

(5)

式中，AMP和φMP分別為多路徑信噪比序列的振幅和延遲相位。若將sinθ視作變量，則：

(6)

由于sinθ非均勻變化，因此使用L-S譜(Lomb-Scargle periodogram, LSP)分析法求取頻率f[12]。但AMP會(huì)隨衛(wèi)星高度角的增大而逐漸減小，僅使用余弦函數(shù)擬合會(huì)遺漏振幅衰減因素，進(jìn)而影響特征參數(shù)的求解。振幅衰減主要與接收機(jī)天線增益和地面反射粗糙程度有關(guān)，因此Strandberg等[13]提出具有阻尼因子的多路徑信噪比模型：

SNRMP=

(7)

式中，Λ為阻尼因子。

2 GA輔助NLS的原理與仿真

2.1 方案原理

遺傳算法(GA)是基于模擬生物遺傳進(jìn)化的優(yōu)化算法，主要包括編碼、適應(yīng)度函數(shù)、遺傳算子(包括種群選擇、交叉、變異)、運(yùn)行參數(shù)4部分。相比于傳統(tǒng)優(yōu)化算法，GA具有2個(gè)特點(diǎn)：1)可全局最優(yōu)解搜索，避免陷入局部最優(yōu)解；2)適應(yīng)度函數(shù)適用范圍廣，約束條件小。

信賴域算法的原理是首先確定一個(gè)信賴域半徑，在該半徑內(nèi)計(jì)算目標(biāo)函數(shù)的二階近似極小值，該近似函數(shù)可表示為：

(8)

式中，q(k)(s)為目標(biāo)函數(shù)f(x)的二次逼近式，s=x-xk；gk為目標(biāo)函數(shù)在xk處的梯度；Bk為xk處Hesse陣。通過二次逼近式與目標(biāo)函數(shù)的比值rk來確定下次迭代的信賴半徑，rk可表示為：

(9)

rk越接近1說明二次逼近式越接近目標(biāo)函數(shù)，則下一次迭代可考慮增大信賴域半徑；當(dāng)‖gk‖≤ε時(shí)，迭代停止。

在進(jìn)行多路徑信噪比延遲相位求解時(shí)，可以通過最小二乘法構(gòu)建參數(shù)方程，其中待定參數(shù)分別為AMP、h、φMP、Λ。其參數(shù)方程可表示如下。

(10)

由于誤差方程為非線性方程組，非線性最小二乘法在對(duì)誤差平方和求其最小值對(duì)應(yīng)的特征參數(shù)時(shí)，使用傳統(tǒng)的信賴域算法容易因初值設(shè)置不當(dāng)而陷入局部最優(yōu)(極小值點(diǎn)而非最小值點(diǎn))，從而影響求解結(jié)果，而僅使用GA時(shí)求解精度不高。因此本文提出采用GA輔助NLS求解最優(yōu)參數(shù)，具體方案為：1)使用二階多項(xiàng)式擬合信噪比趨勢(shì)項(xiàng)，將原始信噪比減去趨勢(shì)項(xiàng)計(jì)算SNRMP；2)將SNRMP和sinθ輸入到具有阻尼因子的函數(shù)模型構(gòu)建誤差方程中；3)使用GA求解誤差方程，計(jì)算初始參數(shù)AMP、h、φMP、Λ；4)將初始參數(shù)賦值給信賴域算法進(jìn)行迭代，獲得最終相位參數(shù)φMP。

2.2 實(shí)驗(yàn)仿真

由于多路徑信噪比延遲相位無法獲得真值，本文利用MATLAB人為設(shè)置參數(shù)來進(jìn)行仿真實(shí)驗(yàn)，對(duì)比2種算法的求解能力。其中，AMP設(shè)置為2，h設(shè)置為190.5 cm，φMP設(shè)置為2.452 5，Λ設(shè)置為46，高度角設(shè)置為5°～20°，均勻生成100組數(shù)據(jù)，然后加入標(biāo)準(zhǔn)差為0.2的隨機(jī)誤差。設(shè)置2種方案：1)方案1，采用傳統(tǒng)余弦函數(shù)模型，利用L-S譜分析方法求解SNRMP和sinθ存在的頻率關(guān)系，賦值初始參數(shù)h，統(tǒng)計(jì)多路徑信噪比最大值并賦值初始參數(shù)AMP，解算相位參數(shù)；2)方案2，采用顧及阻尼因子的信噪比干涉模型，通過GA輔助NLS進(jìn)行求解。

由圖1可見，方案2可以較好地?cái)M合SNRMP曲線，而方案1在高度角過高或過低情況下均無法很好地?cái)M合SNRMP曲線。表1(單位cm)為2種方案100次仿真實(shí)驗(yàn)的相位求解結(jié)果。

圖1 2種方案仿真結(jié)果Fig.1 Simulation result of the two schemes

表1 方案1與方案2誤差統(tǒng)計(jì)

由表1可見，方案2均方根誤差、最大誤差、最小誤差均優(yōu)于方案1，其中RMSE減小32.5%。

3 實(shí)驗(yàn)與分析

3.1 實(shí)驗(yàn)概況

本文實(shí)驗(yàn)數(shù)據(jù)來自PBO網(wǎng)站提供的P043測(cè)站2015年doy182～241觀測(cè)數(shù)據(jù)(共60 d)，所使用的接收機(jī)型號(hào)為Trimble NERT9，同時(shí)使用鋼制三腳架安置，天線罩為SCIT，天線型號(hào)為TRM59800.80。該測(cè)站用于土壤濕度分析的年份較早，具有一定代表性。測(cè)站地勢(shì)平坦、地形無明顯起伏，有利于進(jìn)行土壤濕度實(shí)驗(yàn)。土壤濕度參考值來自https:∥cires1.colorado.edu/portal/index.php，采樣間隔為1 d。

3.2 實(shí)驗(yàn)數(shù)據(jù)處理

通過RTKLIB軟件解算各個(gè)衛(wèi)星的高度角和方位角，并提取L1波段(λ=19.05cm)信噪比。選取觀測(cè)時(shí)間較長(zhǎng)且運(yùn)動(dòng)軌跡較好的衛(wèi)星數(shù)據(jù)(高度角最高可達(dá)70°)，由于測(cè)站周邊地形無明顯差異，故對(duì)不同方位角數(shù)據(jù)不作要求。通過低階多項(xiàng)式擬合去除SNR中的趨勢(shì)項(xiàng)構(gòu)建SNRMP，完成數(shù)據(jù)預(yù)處理。通常情況下，當(dāng)衛(wèi)星高度角設(shè)置為5°～30°時(shí)，土壤濕度反演效果較好[12]，但對(duì)P043測(cè)站各個(gè)衛(wèi)星進(jìn)行去趨勢(shì)項(xiàng)處理后，衛(wèi)星高度角在25°以上時(shí)基本無明顯干涉特征。由于本文所選數(shù)據(jù)月份為7～9月，該期間植被生長(zhǎng)茂盛，過低的高度角數(shù)據(jù)容易受到植被遮擋影響，最終選取5°～25°高度角范圍進(jìn)行實(shí)驗(yàn)。同時(shí)為探究不同高度角范圍的反演結(jié)果，本文設(shè)置3個(gè)高度角取值范圍(5°～15°、5°～20°、5°～25°)進(jìn)行實(shí)驗(yàn)。

圖2為采用2種方案對(duì)G28衛(wèi)星2015年doy241的SNR數(shù)據(jù)進(jìn)行擬合，從圖中可以看出，在處理真實(shí)衛(wèi)星觀測(cè)文件時(shí)，方案2的擬合效果優(yōu)于方案1。

圖2 2種方案擬合結(jié)果Fig.2 Fitting result of the two schemes

由于篇幅限制，3種不同高度角的計(jì)算結(jié)果中，僅給出G14衛(wèi)星5°～15°、G23衛(wèi)星5°～20°高度角范圍內(nèi)2種方案的建模結(jié)果(圖3～4)。由圖可見，G23衛(wèi)星方案1和方案2精度基本相當(dāng)，G14衛(wèi)星方案2精度有較為明顯的提升。在對(duì)G14衛(wèi)星進(jìn)行實(shí)驗(yàn)時(shí)，使用方案1余弦函數(shù)解算時(shí)出現(xiàn)部分異常相位(圖3紅圈)。

圖3 5°～15°高度角范圍下 G14衛(wèi)星2種方案的回歸模型Fig.3 Regression model of the two schemes of G14 satellite at 5° to 15° elevation angle

圖4 5°～20° 高度角范圍下G23衛(wèi)星2種方案的回歸模型Fig.4 Regression model of the two schemes of G23 satellite at 5° to 20° elevation angle

本文統(tǒng)計(jì)不同高度角范圍下2種方案對(duì)G14、G22、G23、G28、G30五顆衛(wèi)星的線性回歸結(jié)果，包括相位和土壤濕度的相關(guān)系數(shù)R、決定系數(shù)R2、均方根誤差RMSE、平均絕對(duì)誤差MAE等，結(jié)果如表2～4所示。

由表2～4可見，5顆衛(wèi)星在3種高度角范圍內(nèi)的15種情況中，5種情況下方案1和方案2精度相當(dāng)，10種情況下方案2較方案1精度有顯著提高，占比66.7%。由此可知，相較于方案1，方案2可有效提高反演精度，而方案1由于延遲相位求解會(huì)出現(xiàn)部分異常值而導(dǎo)致整體結(jié)果變差。對(duì)于G14、G22等信號(hào)狀況較差的衛(wèi)星，方案2在3種高度角范圍內(nèi)均能減小φMP的粗差值，顯著提升反演精度。對(duì)于G23這種信號(hào)較好的衛(wèi)星，2種方案在3種高度角范圍內(nèi)精度相當(dāng)；對(duì)于G28和G30衛(wèi)星，方案2相較于方案1有不同程度提高。相較于方案1，G28衛(wèi)星方案2相關(guān)性提高1.22%～38.33%(由于G28在5°～15°高度角下存在異常值，故不統(tǒng)計(jì)該數(shù)據(jù))，R2提高3.44%～94.10%，RMSE減小2.13%～30.3%，MAE減小5.13%～33.93%。相較于方案1，G30衛(wèi)星方案2相關(guān)性提高3.75%～33.33%，R2提高5.72%～76.06%，RMSE減小6.12%～24.24%，MAE減小2.7%～28.3%。

表2 5°～15°高度角下各方案相關(guān)系數(shù)及擬合精度統(tǒng)計(jì)

表3 5°～20° 高度角下各方案相關(guān)系數(shù)及擬合精度統(tǒng)計(jì)

表4 5°～25° 高度角下各方案相關(guān)系數(shù)及擬合精度統(tǒng)計(jì)

圖5 各衛(wèi)星不同高度角的相關(guān)性Fig.5 Correlation of different elevation angles of each satellite

本文統(tǒng)計(jì)2種方案5顆衛(wèi)星在不同高度角范圍下求解相位和土壤濕度的相關(guān)系數(shù)(圖5)，由圖可見，相較于方案1，方案2受高度角的影響較小，結(jié)果更加穩(wěn)定。表5為3個(gè)高度角范圍內(nèi)的平均相關(guān)系數(shù)以及最大差值，由表可見，相較于方案1，方案2在3種高度角范圍內(nèi)的整體相關(guān)系數(shù)更高，且整體反演結(jié)果更加穩(wěn)定，受高度角范圍變化的影響較小。

表5 各高度角范圍的平均相關(guān)系數(shù)

3.3 多星線性回歸測(cè)試

文獻(xiàn)[14]表明，多星線性回歸可以有效提高反演精度和改善單星反演的跳變現(xiàn)象。本文利用§3.2兩種方案所求的延遲相位進(jìn)行多星線性回歸，以此驗(yàn)證輸入數(shù)據(jù)的優(yōu)劣。由于5°～15°高度角下G14、G22、G28衛(wèi)星的方案1結(jié)果出現(xiàn)一定的延遲相位跳變，從而影響多元線性回歸模型的精度，因此本文僅使用5°～20°和5°～25°高度角下G23-G28、G23-G30、G28-G30、G23-G28-G30四種衛(wèi)星組合方式進(jìn)行驗(yàn)證，結(jié)果見表6。

表6 多星組合線性回歸測(cè)試

由表6可知，相較于單星線性回歸，多星線性回歸反演精度有較為明顯的提高，同時(shí)不同方案輸入的數(shù)據(jù)對(duì)多星線性回歸也具有較為明顯的影響。相較于方案1，5°～20°高度角范圍內(nèi)方案2延遲相位R2提高1.02%～37.85%，RMSE減小2.63%～25%；5°～25°高度角范圍內(nèi)方案2延遲相位R2提高4%～31.25%，RMSE減小6.38%～32.65%。

4 結(jié) 語(yǔ)

1)本文使用GA輔助NLS方法對(duì)顧及阻尼因子的多路徑信噪比模型進(jìn)行延遲相位求解，結(jié)果表明，除G23衛(wèi)星精度與傳統(tǒng)方法差別較小外，其他衛(wèi)星均有不同程度的提升。忽略傳統(tǒng)方法反演產(chǎn)生的部分異常數(shù)據(jù)后，GA輔助NLS方法比傳統(tǒng)方法R2提高5.72%～76.06%，RMSE減小6.12%～24.24%，MAE減小2.7%～28.3%。

2)計(jì)算3個(gè)高度角范圍內(nèi)的相位數(shù)據(jù)可知，GA輔助NLS受高度角范圍影響較小，高度角變化對(duì)相關(guān)性的影響不大，相關(guān)系數(shù)最大差值為0.07，余弦函數(shù)方法最大差值為0.51。

3)GA輔助NLS所求的延遲相位用于多星線性回歸時(shí)精度也有較為明顯的提升，相較于余弦函數(shù)方法，5°～20°高度角范圍內(nèi)RMSE平均減小13.60%，5°～25°范圍內(nèi)RMSE平均減小15.85%。