基于三次樣條插值的農(nóng)機(jī)作業(yè)面積計(jì)算

李亞碩，賈曉峰，伏偉浩，趙 博，王瑞雪，蘇仁忠，徐名漢，龐在溪

(1.中國農(nóng)業(yè)機(jī)械化科學(xué)研究院集團(tuán)有限公司，北京 100083； 2.土壤植物機(jī)器系統(tǒng)技術(shù)國家重點(diǎn)實(shí)驗(yàn)室，北京100083；3.中機(jī)建工有限公司，北京100083； 4.湖北省農(nóng)業(yè)機(jī)械化技術(shù)推廣總站，湖北 武漢 430017)

0 引言

隨著國內(nèi)農(nóng)機(jī)化水平的提高和農(nóng)機(jī)規(guī)模的擴(kuò)大，田間作業(yè)全過程已基本實(shí)現(xiàn)機(jī)械化。農(nóng)機(jī)作業(yè)面積是補(bǔ)貼發(fā)放、費(fèi)用收取、效率評(píng)估、工資發(fā)放等的依據(jù)，精準(zhǔn)有效的作業(yè)面積計(jì)算是農(nóng)機(jī)管理部門和農(nóng)機(jī)合作社管理人員的迫切所需。

在農(nóng)機(jī)上安裝機(jī)載定位設(shè)備可以記錄作業(yè)軌跡信息，每隔幾秒向平臺(tái)發(fā)送一次數(shù)據(jù)，利用北斗定位信息提取行駛軌跡點(diǎn)[1-2]。如何利用這些散點(diǎn)還原農(nóng)機(jī)真實(shí)行駛軌跡，并計(jì)算出真實(shí)作業(yè)面積，是值得關(guān)注和需要解決的問題。魯植雄等[3]利用圖像處理方法計(jì)算作業(yè)面積，同時(shí)標(biāo)注重耕、漏耕區(qū)域，但該方法受軌跡點(diǎn)密集度和分辨率影響較大。有學(xué)者使用柵格法通過計(jì)算軌跡覆蓋區(qū)域，減去空白區(qū)域和重疊區(qū)域統(tǒng)計(jì)作業(yè)面積[4-6]。但柵格大小不好確定，占用內(nèi)存空間較大，并且邊界會(huì)有毛刺。傳統(tǒng)距離測量法是簡單易行的作業(yè)面積計(jì)算方法，將軌跡點(diǎn)連成折線，計(jì)算折線總長度，根據(jù)農(nóng)機(jī)作業(yè)軌跡長度和幅寬計(jì)算。當(dāng)農(nóng)機(jī)沿直線行駛時(shí)，該方法計(jì)算面積誤差不大。當(dāng)農(nóng)機(jī)出現(xiàn)轉(zhuǎn)向或角度發(fā)生變化時(shí)，轉(zhuǎn)向處用折線計(jì)算會(huì)產(chǎn)生誤差。有研究采用矢量緩沖區(qū)算法，將拐角處用圓弧代替，減少了拐角誤差[7-8]。緩沖區(qū)法以幅寬作為轉(zhuǎn)角直徑不能真實(shí)代表所有轉(zhuǎn)向幅度，對(duì)行駛過程中偶然轉(zhuǎn)向有效，當(dāng)?shù)貕K是弧形時(shí)，緩沖區(qū)法仍將兩點(diǎn)間軌跡用直線代替。

本研究基于農(nóng)機(jī)運(yùn)行軌跡特點(diǎn)及作業(yè)地塊不確定性，將農(nóng)機(jī)軌跡用弧線表示，代替?zhèn)鹘y(tǒng)直線法計(jì)算，更能反映農(nóng)機(jī)行駛的真實(shí)情況。本研究設(shè)計(jì)了基于三次樣條插值的作業(yè)面積計(jì)算方法，對(duì)軌跡點(diǎn)進(jìn)行插值，計(jì)算得到的軌跡曲線長度，結(jié)合作業(yè)幅寬得到作業(yè)面積，準(zhǔn)確率較緩沖區(qū)法和距離測量法都有提高。

1 農(nóng)機(jī)作業(yè)面積計(jì)算

農(nóng)機(jī)作業(yè)軌跡分道路行駛、田間行駛和原點(diǎn)停留，有效的作業(yè)面積為田間行駛軌跡，如圖1所示。農(nóng)機(jī)在田間作業(yè)多采用往返作業(yè)形式，即到地塊邊界后掉頭反向作業(yè)。因此，計(jì)算作業(yè)面積時(shí)，除了意外停車軌跡點(diǎn)，也應(yīng)將掉頭轉(zhuǎn)向點(diǎn)去除。保留的軌跡點(diǎn)為連續(xù)單向的線段，將這些線段分別插值計(jì)算再累加。

圖1 農(nóng)機(jī)作業(yè)軌跡Fig.1 Operation track of agricultural machinery

1.1 軌跡點(diǎn)預(yù)處理

機(jī)載設(shè)備在開機(jī)情況下，會(huì)定時(shí)發(fā)送定位信息。有一些特殊點(diǎn)會(huì)影響作業(yè)面積計(jì)算，應(yīng)提前去除。

(1)農(nóng)機(jī)停止散點(diǎn)。農(nóng)機(jī)在停止時(shí)，軌跡點(diǎn)不是集中在一個(gè)點(diǎn)，而是散落在農(nóng)機(jī)附近小范圍內(nèi)，稱為農(nóng)機(jī)停止散點(diǎn)，速度為0或者很小值，如圖2所示，將這些軌跡點(diǎn)直接去除。

圖2 停止散點(diǎn)和軌跡漂移點(diǎn)Fig.2 Dwell scatter point and trajectory drift point

(2)斷點(diǎn)。由于設(shè)備故障問題，兩點(diǎn)之間時(shí)間間隔會(huì)大于正常時(shí)間間隔，稱為數(shù)據(jù)點(diǎn)丟失。將前一斷點(diǎn)作為前段軌跡的終點(diǎn)，后一斷點(diǎn)作為后段軌跡的起點(diǎn)。

(3)漂移點(diǎn)。由于定位精度原因，個(gè)別軌跡點(diǎn)與前后兩點(diǎn)距離遠(yuǎn)大于正常距離間隔，稱該點(diǎn)為漂移點(diǎn)。去除該漂移點(diǎn)，將漂移點(diǎn)前一點(diǎn)作為前段軌跡的終點(diǎn)，后一點(diǎn)作為后段軌跡的起點(diǎn)。

(4)農(nóng)機(jī)在田間作業(yè)一般是往返行駛，到地塊邊界后會(huì)掉頭行駛，地塊邊界轉(zhuǎn)向部分不屬于作業(yè)面積。某軌跡點(diǎn)與前一軌跡點(diǎn)和后一軌跡點(diǎn)分別連成直線，兩直線夾角范圍為0°～180°。當(dāng)夾角大于一定閾值時(shí)，認(rèn)為該軌跡點(diǎn)處于地塊邊界調(diào)轉(zhuǎn)，將該點(diǎn)去除。

所有田間異常點(diǎn)處理后，農(nóng)機(jī)田間作業(yè)軌跡被劃分為多條沿地壟的單向獨(dú)立軌跡，如圖3所示。圖3a為原始軌跡點(diǎn)圖，圖3b為根據(jù)軌跡點(diǎn)速度和轉(zhuǎn)向角度，將地塊轉(zhuǎn)向軌跡點(diǎn)與田間行駛軌跡點(diǎn)標(biāo)記為不同類型。按照時(shí)間序列，計(jì)算每條作業(yè)行地頭之間作業(yè)軌跡點(diǎn)的長度，雖然由于信號(hào)問題導(dǎo)致個(gè)別軌跡點(diǎn)丟失，不影響整體軌跡長度計(jì)算。將計(jì)算得到的各段軌跡累加，得到所有軌跡總長度。

圖3 作業(yè)軌跡預(yù)處理效果Fig.3 Effect of job track preprocessing

1.2 三次樣條插值算法

兩個(gè)點(diǎn)確定一條直線，最簡單的方法就是每兩個(gè)點(diǎn)之間連一條線，最后得到一種折線圖。但是，這種方法顯然具有很明顯的缺陷：曲線不夠光滑，連接點(diǎn)處的斜率變化過大。

如圖4所示，將折線連接處用弧形代替，雖然解決了拐角處誤差，但對(duì)于其余折線P1P2、Pn-1Pn、PnPn+1等處，若地塊形狀如圖所示為弧形，用直線來代替弧形計(jì)算軌跡長度，也會(huì)產(chǎn)生誤差。

圖4 折線連接處用弧線代替Fig.4 Broken line connection is replaced by an arc

如圖5所示，用折線P1P2+P2P3+P3P4+P4P5代替弧P1P5，必然會(huì)出現(xiàn)計(jì)算誤差。因此需要盡可能還原出農(nóng)機(jī)作業(yè)的真實(shí)軌跡，才能更準(zhǔn)確計(jì)算行駛距離和作業(yè)面積。

圖5 兩點(diǎn)之間直線與弧線對(duì)比Fig.5 Comparison of straight line and arc between two points

如圖6所示，離散的坐標(biāo)點(diǎn)為軌跡點(diǎn)，穿過每個(gè)軌跡點(diǎn)的平滑曲線為基于三次樣條插值得到的弧線，未完全經(jīng)過每個(gè)軌跡點(diǎn)的弧形線為最小二乘法擬合出的曲線。計(jì)算農(nóng)機(jī)作業(yè)面積，首先要計(jì)算農(nóng)機(jī)軌跡線長度，每個(gè)采集的軌跡點(diǎn)都是真實(shí)點(diǎn)，數(shù)據(jù)點(diǎn)擬合曲線只能統(tǒng)計(jì)軌跡點(diǎn)趨勢和規(guī)律，不能保證所有點(diǎn)都在統(tǒng)計(jì)內(nèi)，因此三次樣條插值算法得到的弧線更代表農(nóng)機(jī)真實(shí)軌跡[9-11]。通過插值細(xì)分線段，將兩點(diǎn)之間的線段平滑，得到的弧線作為作業(yè)軌跡線。

圖6 擬合和插值方法效果對(duì)比Fig.6 Comparison of fitting and interpolation methods

對(duì)于區(qū)間[a，b]，每個(gè)點(diǎn)坐標(biāo)為[xi，yi]，i=0，1，…，n。將區(qū)間分成x0=a

S0(x)=y0+b0(x-x0)+c0(x-x0)2+d0(x-x0)3
S1(x)=y1+b1(x-x1)+c1(x-x1)2+d1(x-x1)3

(1)

……

Sn(x)=yn+bn(x-xn)+cn(x-xn)2+dn(x-xn)3

其中，S0(x)在區(qū)間[x0，x1]上，S1(x)在區(qū)間[x1，x2]上，……，Sn(x)在區(qū)間[xn-1，xn]上。函數(shù)S(x)若為區(qū)間[a，b]上的三次樣條插值函數(shù)，需滿足以下條件。

(1)在每個(gè)子區(qū)間[xi-1，xi]上S(x)是三次多項(xiàng)式。

(2)S(xi)=yi，i=0，1，…，n。

(3)在區(qū)間[a，b]上S(xi)的二階導(dǎo)數(shù)S″(xi)連續(xù)。

Si(x)=aix3+bix2+cix+di，i=0，1，…，n

(2)

S(xi)=yi，S(xi-0)=S(xi+0)，i=0，1，…，n

(3)

S′(xi-0)=S′(xi+0)

(4)

S″(xi-0)=S″(xi+0)，i=0，1，…，n

(5)

式(2)中的ai、bi、ci、di待定，需滿足式(3)、式(4)和式(5)。

式(3)、式(4)和式(5)共給出4n-6個(gè)條件，需要待定4(n-1)個(gè)系數(shù)，因此要唯一確定三次插值函數(shù)，還要附加2個(gè)邊界條件。

S′(x0)=y1′，S′(xn)=yn′

(6)

S″(x0)=y1″，S″(xn)=yn″

(7)

2 試驗(yàn)及結(jié)果分析

試驗(yàn)采用弧形地塊農(nóng)機(jī)作業(yè)數(shù)據(jù)作為數(shù)據(jù)集，地塊內(nèi)作業(yè)面積為122 544 m2。在農(nóng)機(jī)上安裝采集設(shè)備，每5 s采集一次作業(yè)點(diǎn)數(shù)據(jù)，通過無線發(fā)送至后端服務(wù)器。每個(gè)軌跡點(diǎn)包含農(nóng)機(jī)作業(yè)點(diǎn)經(jīng)緯度、車輛行駛速度等信息。通過三次樣條插值方法計(jì)算農(nóng)機(jī)作業(yè)軌跡長度，乘以農(nóng)機(jī)幅寬得到農(nóng)機(jī)作業(yè)面積。對(duì)比試驗(yàn)為采用折線累加、矢量緩沖區(qū)算法和最小二乘法擬合算法3種方法計(jì)算作業(yè)軌跡長度。

試驗(yàn)步驟如下。

(1)軌跡點(diǎn)預(yù)處理。

(2)利用三次樣條插值法獲取農(nóng)機(jī)作業(yè)軌跡線。

(3)計(jì)算單條軌跡線長度。

(4)累計(jì)軌跡線總長度，乘以幅寬，得到總作業(yè)面積。

準(zhǔn)確率計(jì)算方式為

(8)

式中S——真實(shí)測量面積，m2

S0——計(jì)算面積，m2

采用不同算法計(jì)算的作業(yè)面積結(jié)果如表1所示。本文提出的三次樣條插值法在該地塊的作業(yè)面積計(jì)算中準(zhǔn)確率高于其他幾種方法。作業(yè)緩沖區(qū)算法在處理有轉(zhuǎn)向弧度的作業(yè)軌跡時(shí)，準(zhǔn)確率要優(yōu)于折線累加法，而曲線擬合方法準(zhǔn)確率低于折線累加法。

表1 不同算法計(jì)算作業(yè)面積結(jié)果Tab.1 Calculation results of working area by different algorithms

3 結(jié)論

針對(duì)農(nóng)機(jī)在弧形地塊作業(yè)時(shí)，作業(yè)面積依靠軌跡點(diǎn)連成折線方式計(jì)算誤差較大的問題，提出了利用三次樣條插值法對(duì)軌跡點(diǎn)進(jìn)行插值，將軌跡線用弧線表示。該方法既能保證計(jì)算涵蓋所有軌跡點(diǎn)，也根據(jù)已有軌跡點(diǎn)信息還原真實(shí)農(nóng)機(jī)運(yùn)行軌跡。試驗(yàn)結(jié)果表明，本文方法面積計(jì)算準(zhǔn)確率為97.83%，高于其他面積計(jì)算方法。