景云鵬，劉 剛，金志坤

（1. 中國農(nóng)業(yè)大學(xué)現(xiàn)代精細農(nóng)業(yè)系統(tǒng)集成研究教育部重點實驗室，北京10083； 2. 中國農(nóng)業(yè)大學(xué)農(nóng)業(yè)農(nóng)村部農(nóng)業(yè)信息獲取技術(shù)重點實驗室，北京 100083）

0 引 言

基于GNSS 的農(nóng)田平整技術(shù)對節(jié)水灌溉、土地利用和保肥增產(chǎn)有重要的作用，是現(xiàn)代精細農(nóng)業(yè)的重要研究內(nèi)容之一[1-2]。地形測量是農(nóng)田平整作業(yè)的必要環(huán)節(jié)，通過對未平整農(nóng)田的地勢進行測量，得到準確的位置和高程信息，可用于農(nóng)田平整路徑規(guī)劃，同時為平整工作中挖填土方的計算和作業(yè)效果的評價提供技術(shù)支持[3-5]。

當前，適用于農(nóng)田平整的地形測量設(shè)備主要包括三維激光掃描儀[6-7]、激光雷達[8]、水準儀[9]和 GNSS（global navigation satellite system，全球?qū)Ш叫l(wèi)星系統(tǒng)）。GNSS 可安置在農(nóng)用車輛上實時獲取位置信息，多用于農(nóng)機導(dǎo)航、農(nóng)田采樣和農(nóng)田地圖測繪。國內(nèi)外學(xué)者在應(yīng)用GNSS 測量農(nóng)田地形方面進行了較深入的研究，Clark 等[10]使用RTK-GPS 快速獲取地形數(shù)據(jù)并插值成圖，在113 m2的土地上測量大約30 m 網(wǎng)格的高程數(shù)據(jù)，高程誤差為9 cm。Yao 等[11]利用亞米級水平位置精度的GPS 對農(nóng)田地形圖進行測量，得到的垂直誤差標準差為10～12 cm。李益農(nóng)等[12]利用GPS 代替人工水準儀測量田間高程，通過對比分析得出，根據(jù)GPS 測量數(shù)據(jù)進行土地平整工程設(shè)計，其平地土方總量相對誤差在8%以內(nèi)。孟志軍等[13-14]設(shè)計了基于RTK-GPS 的機載式農(nóng)田三維地形測量系統(tǒng)，并驗證系統(tǒng)在車輛低速且勻速的行駛狀態(tài)下測量精度更高。李宏鵬[15]將RTK-GPS 定位技術(shù)與農(nóng)田地形測量方法集成，開發(fā)了基于GNSS 的農(nóng)田快速地形測量系統(tǒng)，獲得小于1cm 的靜態(tài)高程測量精度。上述研究多側(cè)重于對GNSS 測量的地形數(shù)據(jù)直接進行插值、模型預(yù)測或聚類等分析，最終得到體現(xiàn)地勢信息的地形圖。目前基于GNSS 的農(nóng)田平整系統(tǒng)多使用單天線進行地形測量，通常將定位天線放置平地鏟的中間。這種方法忽略了農(nóng)田平整過程中由于系統(tǒng)機械振動、平地鏟的姿態(tài)改變以及地形起伏等外界因素引起的數(shù)據(jù)采集誤差，降低了農(nóng)田三維地形測量精度。

本文對當前的GNSS 單天線農(nóng)田測量方法進行改進，采用GNSS 雙天線定位定向技術(shù)獲取農(nóng)田地勢信息，使用AHRS 測量的加速度數(shù)據(jù)校正分析GNSS 得到的農(nóng)田地勢信息，分析了系統(tǒng)結(jié)構(gòu)與地形變化引起的測量誤差。研制了集成地勢信息獲取與誤差處理功能于一體的地形測量系統(tǒng)，并進行了農(nóng)田地形測量試驗。

1 地形測量原理

1.1 地勢數(shù)據(jù)獲取方法

分別將GNSS 定位、定向天線放至平地鏟車的兩端，將姿態(tài)航向傳感器MTI300_AHRS 通過墊片固定在同一水平面，其空間位置分布如圖1 所示。定位天線、定向天線和姿態(tài)航向傳感器位于同一基線O1O2上，雙天線分別位于鏟車2 個支撐輪的中心線與O1O2的垂直交點上，距離為L。地形測量過程中，將測量車輪行駛路徑上的地勢點作為輸入量，由于定向天線無法輸出位置信息，根據(jù)定位天線的空間位置進行坐標轉(zhuǎn)換和平移，獲得另一側(cè)車輪所對應(yīng)定向天線的位置信息。

圖1 地形測量平臺的組成 Fig.1 Composition of topographic survey platform

地形數(shù)據(jù)的動態(tài)采集過程中，平地鏟通過牽引式結(jié)構(gòu)與拖拉機相連接。由于地形變化的影響，2 個衛(wèi)星天線的姿態(tài)會隨著平地鏟一起發(fā)生變化[16-17]。在三維坐標系中通過平地鏟左端定位天線的地勢測量數(shù)據(jù)得到其位置信息矩陣為M'=(XM'，YM'，ZM')T，為求得定位天線在世界坐標系中的位置矩陣M=(XM，YM，ZM)T，根據(jù)AHRS實時測得的翻滾角、俯仰角、航向角信息，建立GNSS天線與AHRS 繞X、Y、Z 軸旋轉(zhuǎn)的矩陣RX、RY、RZ，同時可計算出總旋轉(zhuǎn)矩陣R[18-19]：

式中θR為翻滾角，(°)；θP為俯仰角，(°)；θY為航向角，(°)。

定位天線的位置坐標矩陣M=RM’，為求得另一側(cè)支撐輪所對應(yīng)的定向天線的位置坐標矩陣N=(XN，YN，ZN)T，再通過矩陣平移變換[20]計算：

式中S 為定位天線向定向天線做平移變換時的平移矩陣，XM、YM、ZM為GNSS 定位天線在三維坐標系中的坐標值，m；

GNSS 雙天線也可測得自身的航向角，其角度以正南方向為0°、正西方向為90°、正北方向為180°、正東方向為270°，角度順時針遞增。將經(jīng)緯度轉(zhuǎn)換成平面坐標時，以正東方向為X 軸，正北方向為Y 軸構(gòu)建平面坐標系。結(jié)合式（1）～（3），可得到GNSS 雙天線的位置解算方程組：

式中TX為GNSS 定位天線沿X 軸平移至定向天線的距離，m；TY為GNSS定位天線沿Y軸平移至定向天線的距離，m；TZ為GNSS定位天線沿Z軸平移至定向天線的距離，m；θy為GNSS 雙天線測得的航向角，(°)。

1.2 誤差分析與處理

研究發(fā)現(xiàn)，基于GNSS 的農(nóng)田地勢信息采集誤差主要包括微地形誤差和機械振動誤差。微地形誤差即由于農(nóng)田局部凹凸變化而產(chǎn)生的測量誤差，前人的研究中認為地形測量時的農(nóng)田表面是連續(xù)的，即整體趨勢呈高低起伏[21-22]。但實際上農(nóng)田的表面并非光滑，局部地勢不斷變化。未平整之前存在的凹凸、小坑和小包等，微地形變化會使得采集的數(shù)據(jù)偏離真實值。機械振動誤差是在拖拉機牽引平地鏟行駛過程中產(chǎn)生的，由于拖拉機與平地鏟屬于銷軸式半鋼性連接，且缺少減震裝置，平地鏟會產(chǎn)生無規(guī)則的機械振動。該隨機振動產(chǎn)生的位移偏差疊加至衛(wèi)星天線采集到的位置信息中，進而影響測量數(shù)據(jù)精度。

本文將微地形誤差與機械振動誤差統(tǒng)稱為振動誤差，振動誤差可以通過振動加速度進行衡量。計算分析AHRS 獲取的平地鏟三軸加速度信息，得出的位移變化信號即為振動誤差[23-24]，本文振動誤差的處理步驟如下：

1）振動加速度的測量。將AHRS 測得的平地鏟三軸振動加速度Xa 、Ya 、Za 通過角度與位移變換（坐標變換矩陣詳見1.1 節(jié)），可得到定位天線和定向天線的振動加速度。設(shè)AHRS 獲取的振動加速度為A=(aX，aY，aZ)T，則：

式中aX、aY、aZ分別是在X、Y、Z 軸上的振動加速度；AM和AN分別是GNSS 定位天線和定向天線的振動加速度；A 是AHRS 測量的振動加速度，m/s2。

2）加速度趨勢項預(yù)處理。實際測量過程中由于溫度變化、農(nóng)田濕度、電力磁場、GNSS 信號干擾等因素，會導(dǎo)致AHRS 測量的振動信號產(chǎn)生零點漂移和偏離基線等現(xiàn)象。這些現(xiàn)象會降低振動加速度的測量精度，為消除這些趨勢項的影響，本文采用最小二乘法對振動加速度進行處理。設(shè)采集的振動加速度數(shù)據(jù)為振動擬合系數(shù)為根據(jù)最小二乘法原理可知，將m 階多項式kc 與振動加速度信號進行擬合，以多項式與原始振動信號的誤差平方和最小為目標，計算求得對應(yīng)的擬合系數(shù)ib 。最后將原始振動加速度與對應(yīng)多項式做差，即可得到預(yù)處理后的振動加速度，即：

式中ak為第k 個點的振動加速度，m/s2；ak’為第k 個點消除趨勢項后的振動加速度，m/s2；n 為正整數(shù)；bi為第i階的擬合系數(shù)；ck為第k 個點的擬合多項式；m 為非負整數(shù)。

3）基于頻域的振動位移求解。將消除趨勢項后的振動加速度進行二次積分可獲得振動位移信號。前期研究表明，由于在時域中進行二次積分會產(chǎn)生新的誤差趨勢項從而影響積分效果，本文選用基于快速傅里葉變換（fast fourier transform, FFT）的頻域積分方法對振動位移進行求解。計算步驟如下：

一次積分求得速度信號：

二次積分求得振動位移信號：

式中t 為采樣時間，s；ω 為角頻率，rad/s；a(t)為時域中的加速度，m/s2；A(ω)為頻域中的加速度，m/s2；v 是速度，m/s；s是位移，m；，j為單位虛數(shù)；δ(ω) 為脈沖函數(shù)。

由于振動信號為離散數(shù)據(jù)，若求積分值，需先求得頻域中的角頻率

式中N 為振動信號的采集總數(shù)，f 為采集頻率，Hz；q 為積分次數(shù)，則積分的頻域變換為

積分的相位變換為

式中Fk表示第k 個點在頻域的振動加速度；Fk’為Fk積分變換之后；虛數(shù)單位j 在頻域中表示相移，D 為Fk的實部，B 為Fk的虛部。每次乘以一個j 就逆時針旋轉(zhuǎn)90°，除以一個j 順時針旋轉(zhuǎn)90°。結(jié)合式（8）～（9）可知，在頻域中一次積分順時針旋轉(zhuǎn)90°，二次積分順時針旋轉(zhuǎn)180°。將最終得到的復(fù)數(shù)數(shù)組通過傅里葉反變換至?xí)r域，即可得到所求振動位移。

4）數(shù)據(jù)平滑處理。田間試驗環(huán)境會對測量儀器造成較大影響，長時間采集的高程信號會出現(xiàn)毛刺現(xiàn)象，對最后求得的高程值進行數(shù)據(jù)平滑處理有利于提高高程曲線的信噪比與光滑度[25]。由于系統(tǒng)數(shù)據(jù)為等間距采集，以時間為自變量，則時間間隔為定值。本文選用基于最小二乘的五點三次平滑法進行處理，設(shè)不同時刻所得到的高程值為五點三次平滑法的表達式為

式中sk為第k 個點的高程值，m。

2 GNSS 雙天線結(jié)合AHRS 測量農(nóng)田地形的硬件平臺與軟件算法

2.1 硬件平臺

本文以中國農(nóng)業(yè)大學(xué)精細農(nóng)業(yè)研究中心前期開發(fā)的GNSS 智能農(nóng)田平整系統(tǒng)為試驗平臺，使用John Deere 公司的5-904 拖拉機作為牽引動力。采用北京盛恒天寶公司生產(chǎn)的平地鏟（具有獨立式液壓結(jié)構(gòu)，鏟車寬度為2.5 m），上海華測導(dǎo)航公司生產(chǎn)的GNSS 雙天線導(dǎo)航定位系統(tǒng)（平面測量精度為0.8 cm+1mm，高程測量精度為1.5 cm+1mm，航向精度為0.1°）；AHRS（荷蘭Xens 公司生產(chǎn)型號為MTi-300）橫滾角和俯仰角的靜態(tài)測量精度為0.2°，航向角的測量誤差為1°，加速度的靜態(tài)測量誤差為0.02 m/s2。GNSS 農(nóng)田平地機及地形測量平臺的組成如圖2 所示。

2.2 軟件算法

為保證數(shù)據(jù)運算速度，采用Visual Studio 2008 的C++語言作為系統(tǒng)的編譯環(huán)境，利用Matlab Compiler 編譯器將繪制農(nóng)田三維地形圖的M 函數(shù)文件轉(zhuǎn)換為Visual環(huán)境可調(diào)用的形式，以提高軟件系統(tǒng)的集成度，圖3 為軟件程序的流程圖。

圖2 GNSS 農(nóng)田平地機及地形測量平臺的組成 Fig.2 Composition of GNSS land leveler and topographic survey platform

圖3 農(nóng)田三維地形測量程序流程圖 Fig.3 Flow chart of farmland 3D topographic survey program

3 地形測量試驗與結(jié)果分析

3.1 水泥路面試驗與結(jié)果分析

為檢驗在機械振動影響下，GNSS 雙天線位置解算方程和誤差處理方法的效果，本文于2019 年9 月底在中國農(nóng)業(yè)大學(xué)上莊試驗站進行水泥路面的試驗驗證。分別選取長度為300、900 和1800 m 的水泥路面，拖拉機的行駛速度保持在1.12 m/s，數(shù)據(jù)采集頻率為5 Hz。3 條路徑中定向天線通過位置解算和誤差處理前后的高程測量效果如圖4 所示，表1 為定位、定向天線的高程測量結(jié)果。并以誤差和來評價通過誤差處理后定位天線和定向天線的高程變化情況[26-27]：

式中E 為誤差和，m；Hi為第i 個測點誤差處理前的高程值，m；hi為第i 個測點誤差處理后的高程值，m。

圖4 不同長度路徑下定向天線位置解算和誤差處理前后的高程測量結(jié)果 Fig.4 Elevation measurement result of directional antenna before and after position calculation and error processing for different length road

由圖4 可知，在水泥路面行駛過程中由于拖拉機啟動帶動平地鏟發(fā)生機械振動，定向天線所測得的高程信息上下浮動。結(jié)合表1 可知，300 和900 m 路面長度的地形起伏較小，路面較平坦。定位天線經(jīng)過誤差處理后的誤差和分別為1.986 和2.063 m，數(shù)據(jù)波動范圍分別減小了4.67%和9.57%；定向天線經(jīng)過誤差處理后的誤差和分別下降了1.909 和2.143 m，數(shù)據(jù)波動范圍分別減小了8.42%和8.22%。當路面長度為1 800 m 時，數(shù)據(jù)采集時間更長，系統(tǒng)所處地形起伏較大，經(jīng)過處理后定位天線的誤差和為5.228 m，數(shù)據(jù)波動范圍減小了1.52%；定向天線的誤差和為5.956 m，數(shù)據(jù)波動減小了1.1%。測試結(jié)果證明本文的方法可以有效獲取GNSS 雙天線的位置信息，使定位精度得到較好的改善。

表1 定位、定向天線的高程測量結(jié)果 Table 1 Elevation measurement results of positioning and directional antenna on cement road

3.2 農(nóng)田試驗與結(jié)果分析

為進一步驗證在農(nóng)田地形起伏環(huán)境下，農(nóng)田地形測量和誤差處理方法的效果，選用中國農(nóng)業(yè)大學(xué)上莊試驗站的3 塊農(nóng)田進行標準區(qū)域劃分[28-30]，并手動測量農(nóng)田的地勢信息。將手持GNSS 測量的地勢信息作為真值，單天線測量的地勢信息作為對照，與GNSS 雙天線/AHRS 地形測量系統(tǒng)的測量結(jié)果進行對比。

首先將拖拉機啟動并停至農(nóng)田中，如圖5 所示，AHRS 固定在至鏟車頂部，設(shè)定采集頻率為5 Hz（適合農(nóng)田平整作業(yè)平地鏟的最佳頻率）。分別在5、15 和 30 min 內(nèi)獲取AHRS 測得的垂直方向加速度、航向角、俯仰角和翻滾角等主要信息，以驗證AHRS 在農(nóng)田環(huán)境下的測量精度，具體的試驗結(jié)果如表2。

圖5 AHRS 精度標定 Fig.5 Precision calibration of AHRS

表2 AHRS 精度標定結(jié)果 Table 2 Precision calibration result of AHRS

由表2 中可知，系統(tǒng)在農(nóng)田環(huán)境下雖然受系統(tǒng)結(jié)構(gòu)和地形的影響，隨著測量時間增加至30 min，AHRS 獲取的主要信息波動范圍減小。結(jié)合表2 中數(shù)據(jù)可知，AHR 垂直加速度的均方根誤差（root mean square error, RMSE）小于0.256 m/s2，航向角、俯仰角和翻滾角的均方根誤差分別為0.903°、0.849°和0.531°，將誤差角度帶入式（4）中的空間位置方程，計算可知 |LsinθP| ≤2.9 cm（測量可得雙天線之間的距離L=2 m），對GNSS 天線的測量結(jié)果影響較小，證明文中所選的AHRS可應(yīng)用于復(fù)雜農(nóng)田環(huán)境下的地形測量。

農(nóng)田地形測量試驗分為2 個部分：

1）機載地形測量試驗。選取3 塊未被耕作的農(nóng)田，面積分別為35 m×50 m、35 m×100 m 和35 m×200 m。為保證GNSS 定位數(shù)據(jù)盡可能地遍歷整個農(nóng)田并結(jié)合平臺自身體積影響，如圖6a所示，將每塊試驗田劃分成寬5 m的7 個區(qū)域。試驗人員駕駛拖拉機分別使用GNSS 雙天線/AHRS地形測量系統(tǒng)和單天線GNSS智能農(nóng)田平整系統(tǒng)以1.12 m/s 的速度依次駛過這些區(qū)域，采集頻率5 Hz，系統(tǒng)每次采用地頭轉(zhuǎn)彎的方式駛?cè)胂乱粋€區(qū)域。

2）真實地勢信息采集試驗。如圖6b 所示，真實地勢信息由試驗人員手持GNSS 衛(wèi)星定位天線對試驗農(nóng)田進行測量。通過手持裝有GNSS 定位天線的加長桿，沿著拖拉機行駛的路徑每隔20 cm進行定點測量，3 塊農(nóng)田中分別采集到3 500、7 000 和14 000 個點（與單、雙天線采集頻率相對應(yīng)）。由于機載地形測量平臺的雙天線測量路徑對應(yīng)平地鏟車輪行駛路徑，人工測量時GNSS 的數(shù)據(jù)采集路徑與機載測量的14 條路徑重合，以保證對比試驗的可靠性。

圖6 農(nóng)田地形測量試驗 Fig.6 Field topographic survey test

將不同測量方式得到的地勢信息通過反距離加權(quán)（inverse distance weighting, IDW）插值法制圖，結(jié)果如圖7 所示。采用平整度、最大高程差和高程分布列[31-33]對測量結(jié)果進行評價，結(jié)果如表3 所示。

由圖7 可以看出，原始單天線測量系統(tǒng)測量得到的三維地形圖較為粗糙，雙天線地形測量系統(tǒng)得到的地形圖更為光滑，且更接近手持測量的真實地勢信息。結(jié)合表3可知，經(jīng)過位置方程解算得到的定位天線、定向天線的采集點數(shù)是單天線采集點數(shù)量的2 倍。相比于單天線，在3 塊不同面積農(nóng)田中，雙天線系統(tǒng)測量得到的農(nóng)田平整度分別由3.4 cm 下降到3.0 cm，由5.3 cm 上升到6.2 cm，由10.6 cm 上升到11.8 cm，與真值的平整度相比，準確性分別提高了14.286%、14.063%和10.084%；最大高程差分別由18.1 cm 下降到16.7 cm，由27.6 cm 下降到25.6 cm，由54.8 cm 下降到50.2 cm，與真值的最大高程差相比，準確性分別提高了8.642%、8.333%和8.897%；高差分布列分別由90.532%下降到85.650%，由55.677%下降到50.916%，由27.015%下降到24.024%，與真值的高差分布列相比，準確性分別提高了 1.536%、3.357%和2.991%。試驗結(jié)果表明，雙天線地形測量系統(tǒng)得到的地勢信息更接近真值，地形測量精度高于單天線測量。

圖7 不同測量方式的地形測量結(jié)果對比 Fig. 7 Comparison of topographic survey results with different methods

表3 手持GNSS、單天線和雙天線的方式測量3 塊農(nóng)田的地勢信息 Table 3 Topography information of three farmland measured by hand-held GNSS , single antenna and dual antenna

4 結(jié) 論

1）建立了基于GNSS 農(nóng)田平整系統(tǒng)的雙天線位置解算方程，在農(nóng)田平整作業(yè)過程中可以有效的解算出定位、定向天線的位置信息。

2）提出了農(nóng)田地勢信息誤差處理方法。本文通過分析地形測量時誤差信號的來源和性質(zhì)，以高精度AHRS為載體獲取平地鏟的姿態(tài)和加速度信號。利用最小二乘法對振動加速度進行預(yù)處理、快速傅里葉變換進行誤差位移轉(zhuǎn)換；五點三次法進行數(shù)據(jù)平滑處理，提高了GNSS 天線獲取數(shù)據(jù)的光滑度和準確度。

3）進行了水泥路面和田間測量試驗，試驗過程中拖拉機的行駛速度為1.12 m/s。水泥路面的試驗結(jié)果表明，雙天線位置解算方程在動態(tài)環(huán)境下可以獲得定位天線與定向天線的位置信息。經(jīng)過誤差處理后定位天線的數(shù)據(jù)波動最大可減少9.57%，定向天線的數(shù)據(jù)波動最大可減少8.42%，適用于基于GNSS 的農(nóng)田平整作業(yè)。田間測量試驗結(jié)果表明，農(nóng)田環(huán)境中AHRS垂直加速度的均方根誤差為0.256 m/s2，航向角、俯仰角和翻滾角的均方根誤差分別為0.903°、0.849°和0.531°。與傳統(tǒng)單天線農(nóng)田地形測量方法對比，3 種地塊面積（35 m×50 m、35 m×100 m 和35 m×200 m）下，雙天線測量的農(nóng)田平整度準確性分別提高了14.286%、14.063%和10.084%；最大高程差的準確性分別提高了8.642%、8.333%和8.897%。高差分布列的準確性分別提高了1.536%、3.357%和2.991%。雙天線地形測量系統(tǒng)可以應(yīng)用于精細農(nóng)業(yè)中的地形測量，可為GNSS 農(nóng)田平整作業(yè)提供精度高、持續(xù)穩(wěn)定的農(nóng)田地勢信息。