GA輔助NLS的GNSS-IR土壤濕度反演方法

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

1 桂林理工大學測繪地理信息學院，桂林市雁山街319號，541006 2 廣西空間信息與測繪重點實驗室，桂林市雁山街319號，541006 3 廣西建設職業技術學院土木工程學院，南寧市羅文大道33號，530007

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

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

1 GNSS-IR基本原理

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

(1)

由于接收機實際接收到的為干涉信號，因此可將信噪比公式近似表達為：

(2)

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

(3)

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

(4)

式中，φ為相對延遲相位。

接收機天線增益會隨衛星高度角的增大而增大，接收的直射信號強度會隨衛星高度角的增大而提高，而反射信號強度會隨衛星高度角的增大而下降，導致整個信噪比序列呈現二次曲線變化趨勢。通常使用低階多項式擬合來表示該趨勢，通過將原始信噪比序列減去趨勢項來構建多路徑信噪比SNRMP。低高度角范圍下多路徑信噪比可表示為余弦函數形式：

(5)

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

(6)

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

SNRMP=

(7)

式中，Λ為阻尼因子。

2 GA輔助NLS的原理與仿真

2.1 方案原理

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

信賴域算法的原理是首先確定一個信賴域半徑，在該半徑內計算目標函數的二階近似極小值，該近似函數可表示為：

(8)

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

(9)

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

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

(10)

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

2.2 實驗仿真

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

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

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

表1 方案1與方案2誤差統計

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

3 實驗與分析

3.1 實驗概況

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

3.2 實驗數據處理

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

圖2為采用2種方案對G28衛星2015年doy241的SNR數據進行擬合，從圖中可以看出，在處理真實衛星觀測文件時，方案2的擬合效果優于方案1。

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

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

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

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

本文統計不同高度角范圍下2種方案對G14、G22、G23、G28、G30五顆衛星的線性回歸結果，包括相位和土壤濕度的相關系數R、決定系數R2、均方根誤差RMSE、平均絕對誤差MAE等，結果如表2～4所示。

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

表2 5°～15°高度角下各方案相關系數及擬合精度統計

表3 5°～20° 高度角下各方案相關系數及擬合精度統計

表4 5°～25° 高度角下各方案相關系數及擬合精度統計

圖5 各衛星不同高度角的相關性Fig.5 Correlation of different elevation angles of each satellite

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

表5 各高度角范圍的平均相關系數

3.3 多星線性回歸測試

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

表6 多星組合線性回歸測試

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

4 結 語

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

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

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