基于三次樣條法的林地面積量算方法比較研究

基于三次樣條法的林地面積量算方法比較研究

面積測量是森林資源研究的基礎工作。區(qū)域測量的方法、方法和影響直接影響森林資源的研究和監(jiān)測方法、方法和影響。同樣，區(qū)域測量的精度直接關系到森林資源調查監(jiān)測結果的精度和可靠性。因此，根據不同的原則、不同的方法和不同的用途，選擇合適的測量方法，確保必要的測量精度，是測量工作者長期研究的主題。縱觀國內外對于土地面積的量測,從作業(yè)過程上講有兩大類:一類是野外直接測量法,如采用DGPS、羅盤/電子羅盤配合皮尺、全站儀、經緯儀/電子經緯儀配合光電測距儀/地形尺等;另一類則是從地圖上量測.從圖形形態(tài)上分兩大類,一類面向規(guī)則圖形的解析算法和圖解算法,一類是面向任意不規(guī)則圖形的量算方法.實際工作中,選擇哪一類方法,主要取決于圖形形態(tài)、現有設備和工具、要求的必要精度等.就森林資源調查、監(jiān)測中量測面積的現狀而言,常采用以下3種方法:①儀器實測法——主要應用林界勘測、林區(qū)地籍測量、伐前后調查、火災后調查等.②求積儀法——主要應用于林業(yè)基本圖、林相圖、資源圖、規(guī)劃圖上林班、小班面積的量算.③數字電子地圖法——可用直接野外數字化測圖或原圖的掃描或手扶跟蹤數字化圖.前者,直接采用儀器原始測量坐標求算面積;后者則依據用圖者圈定的地界,自動解析折點坐標進而求算面積.上述3種方式解算面積在精度上的主要問題是:①儀器實測法對于曲線邊界的不規(guī)則圖形面積計算失真.②數字電子地圖上對于曲線邊界的不規(guī)則圖形面積計算失真.③求積儀法確定面積的極限精度問題.④儀器實測面積精度與圖上量測面積精度的關系問題.1量測面積誤差2儀器實測面積誤差σs1的一般公式為:σ2s12s1=σ2儀2儀+σ2系12系1(1)式中,σ儀為測量儀器量測圖形拐點坐標、邊長、角度誤差引起的面積誤差,多屬偶然性質的隨機誤差;σ系1為以拐點之間的折線代替弧線引起的面積偏差,屬系統(tǒng)誤差,對折線連接的圖形,一般可認為σ系1=0.圖上量測面積誤差σs2的一般公式為:σ2s22s2=σ2s12s1+σ2量2量+σ2系22系2+σ2系32系3(2)式中,σ量為量測圖上面積所引起的誤差,如求積儀安置誤差、描跡誤差、計算誤差、計算機數字化地圖量算時的捕點移位誤差等,均屬偶然性隨機誤差;σ系2為描述偏離圖斑界線所引起的誤差,計算機數字化地圖上以折代曲量算誤差,一般屬于系統(tǒng)性質的誤差;σ系3為圖板不平、圖紙變形引起的誤差,對于數字電子地圖,σ系3=0,σ系3屬于系統(tǒng)性質的誤差.基于以上分析,不難得出結論:當野外測量圖斑的方法確定后,必然有σs1<σs2,即實測面積精度要高于任何手段獲取的圖上圖斑面積精度.因此,實際應用中,最可信的是野外數字測量直接獲得的面積精度.但這種方法需要到現地進行量測,這對于大面積林地面積的量測將會使工作量過大.2提高曲線圍成圖斑面積相對精度的有效手段以折線代替曲線簡化曲線圍成的圖斑,利用圖解解析、實測數據(如坐標、邊長、夾角等)計算面積顯然會帶來一定的系統(tǒng)誤差.現用坐標解析公式為例加以說明.為研究問題方便,假設正n邊形圖斑,其面積為s0,其外接圓半徑為R,點位誤差為σ位,得s0=n2R2sin2πnσs0=√n2σ位Rsin2πn(3)s0=n2R2sin2πnσs0=n2?√σ位Rsin2πn(3)式中,σs0為測量誤差.即正n邊形面積由于測量誤差引起的面積相對誤差為:σs0s0=√2n(σ位R)(4)σs0s0=2n?√(σ位R)(4)即n越大,σ位越小,σs0so越小,相對精度越高.若圖斑為一半徑R的圓,用接近正n邊形觀測時,則有:(σss)2=(σs1πR2)2+(σs0s0)2=[√n2sin2πnπ(σ位R)]2+[√2n(σ位R)]2(5)即n越大,σ位越小,σs0s越小,相對精度越高.當n較大時,sin2πn≈2πn,(σss)=2√n(σ位R).由此得出結論,提高曲線圍成圖斑面積相對精度的有效手段是:①增加點個數或提高測點密度;②使σ位R減小,即提高點位精度和加大圖形面積;③實際工作中,為使σ位減小,必須采用高精度的儀器測量坐標.當采用以C/A碼為觀測值RTDGPS測量曲線圍成的面積時,除適當加大測點的密度外,可用皮尺加測相鄰兩側點的邊長,按秩虧自由網平差求得各點坐標,進而使σ位減少.對于全站儀定位、RTKGPS定位則不需加測邊長和實行秩虧自由網平差,因為其σ位一般為cm級水平.顯然,實際工作中,無限增大測點數會導致龐大的工作量,因而不可取.采用三次樣條(CubicSpline)法計算曲線圍成面積,將得到良好的效果.3次樣條及其同區(qū)間的配合關系利用解析坐標、圖解坐標和實測坐標計算由曲線圍成的多邊形面積多采用梯形公式法,也可采用辛普生法或等距三次樣條法.這些方法的主要缺點是:①在梯形法中,假定不規(guī)則邊界每段為一直線,在辛普生法中,每段被假設為二次拋物曲線;而在等距三次樣條法法中,每段被假設為一條三次多項式曲線.等距三次樣條法的數學模型較接近于實際地面不規(guī)則的邊界.盡管如此,這個模型仍然有一個缺點:每個都假設是三次多項式的任意兩段的連接處將形成尖峰點,不可能很好地與地面曲線匹配.②在梯形法、辛普生法和等距三次樣條法中,必須在x(或y)方向等向距RTKGPS定位.所以,我們不能夠根據地面曲線的變化選擇合適的點來測量這些點,而在這種條件下建立的模型不能夠很好地匹配地面曲線,一種遵循這種不完善模型的方法是不可能精確的.此外,在野外按基線等區(qū)間地測量垂距是非常不方便的.為了避免這些問題,這里介紹一種新的方法,采用此法我們不僅能夠根據地面曲線的變化選擇測量垂距的點,而且還應用了比梯形法、辛普生法或等距三次樣條法能更好地擬合地面曲線的數學模型(見圖1).一條不規(guī)則的曲線ab被分成為n個相等或不等的區(qū)段[X0,X1],[X1,X2],…,[Xn-1,Xn].在每個區(qū)間上應用一個三次樣條多項式來配合區(qū)間[X0,Xn]上的這條不規(guī)則的曲線.三次樣條多項式定義如下:已知定義在[a,b]上的一個函數f和一組結點a=X0<X1<…<Xn=b,對于f一個三次樣條內插式P是滿足下列條件的函數:(a)P是一個三次多項式,對每個i=0,1,…,n-1在區(qū)段[Xi,Xi+1]上表示為Pi;(b)對每個i=0,1,…,n有P(Xi)=f(Xi);(c)對每個i=0,1,…,n-2有P(Xi+1)=Pi(Xi+1);(d)對每個i=0,1,…,n-2有P′i+1(Xi+1)=P′i(Xi+1);(e)對每個i=0,1,…,n-2有P″i+1(Xi+1)=P″i(Xi+1);(f)P″(X0)=P″(Xn)=0(自由邊界).為了構作函數f的三次樣條內插式P,定義中的條件可以應用到每個三次多項式Pi(x)上,對每個i=0,1,…,n-1,寫為下面的形式:Pi(x)=ai+bi(X-Xi)+ci(X-Xi)2+d(X-Xi)3(6)對每個i=0,1,…,n-2,根據條件(b)有Pi(xi)=ai=f(xi)(7)而根據條件(c),對每個i=0,1,…,n-1設hi=Xi+1-X,并定義an=f(Xn).對每個i=0,1,…,n-1,式(7)可以重寫為ai+1=ai+bihi+cih2i+dih3i(8)按類似的方式定義bn=P′(Xn).且P′(Xn)=bi+2cihi+2dih2i(9)則對每個i=0,1,…,n-1有P′i(Xi)=bi.應用條件(d)對每個i=0,1,…,n-1得出bi+1=bi+2cihi+2dih2i(10)Pi各系數之間的另外一種關系可以求得如下,設P″i(X)=2ci+6di(X-Xi)則P″i(X)=2ciP″i+1(Xi+1)=2ci+1P″i+1(Xi+1)=2ci+6di(X-Xi)對每個i=0,1,…,n-1應用條件(e)有ci+1=ci+3dihi(11)我們還規(guī)定cn=P″(Xn)/2.在式(11)中解算di,得di=ci+1-ci3hi(12)將式(12)代入式(8)中,并解算bibi=(ai+1-ai)/hi-hi(2ci+ci+1)/3(13)然后在下標減小情況下對bi-1有bi-1=(ai-ai-1)/hi-1-hi-1(2ci-1+ci)/3(14)將式(12)代入式(10)中并把下標減1得bi=bi-1+hi-1(ci-1+ci)(15)將式(13)和式(14)代入式(15)中,對于i=0,1,…,n-1給出了線性方程組hi-1ci-1+2(hi-1-hi)ci+hici+1=3(ai+1-ai)/hi-3(ai-ai-1)/hi-1(16)在式(16)中,由于{hi}n-1i=0和{ai}ni=0點組{Xi}ni=0的間隔和f在這些點上的值已給定,所以{ci}ni=0的值能夠計算出來.由式(13)求出常數{bi}n-1i=0的余項,由式(12)并構作三次樣條多項式{Pi}n-1i=0也就簡單了.利用式(14)和自由邊界條件P″(a)=P″(b)=0,我們能夠求得(n+1)個未知數{ci}ni=0.用通常表示法a=X1<X2<…<Xn=b,該自由邊界條件意味著cn=P″(Xn)/2=0和P″(X0)=2c0+6d0(X0-X0)=0,所以c0=0.用矩陣來表示式(16)和自由邊界條件,我們能夠寫出H×C=A(17)式中,Η={100???0h02(h0+h1)h1????0hi2(h1+h2)h2????????????????00??0hn-22(hn-2+hn-1)hn-10???001}是一個(n+1)×(n+1)的矩陣;C={c0c1?cn}是一個(n+1)×1的矩陣;A=[03(a2-a1)/h0-3(a1-a0)/h1?3(an-an-1)/hn-1-3(an-1-an-2)/hn-2]是一個(n+1)×1的矩陣,并且ai=f(Xi).解式(17),我們得到C=H-1A(18)所以Ci=H-1iA(19)式中,H-1i是一個1×(n+1)的行矩陣,由H-1的第i行元素組成,而H-1是矩陣H的逆陣.將式(19)代入式(13)中,我們得到bi=(ai+1-ai)/hi-2hiHi-1+H-1i+1A/3(20)再把式(19)代入式(12)中,我們得到di=(H-1i+1-H-1i)A/3(21)所以,利用式(7),(18),(21)和(20)確定三次樣條多項式的全部系數就能夠得出該三次樣條多項式的模型.如圖1所示,在X軸和曲線ab之間劃斜線的區(qū)域,區(qū)間[Xi,Xi+1]上的三次樣條多項式可以寫成式(22)Pi(X)=ai+bi(X-Xi)+ci(X-Xi)2+di(X-Xi)3(22)而面積Ai可用積分計算如下:A=∫xi+1xiΡi(X)dX=∫xi+1xi[ai+bi(X-Xi)+ci(X-Xi)2+di(X-Xi)3]dX=(ai-biXi+ciX2i-diX3i)(Xi+1-Xi)+(bi/2-ciXi+3diX2i/2)(X2i+1-X2i)+(ci/3-diXi)(X3i+1-X3i)+(di/4)(X4i+1-X4i)(23)利用式(23)可計算區(qū)間[Xi,Xi+1]內曲線ab和X軸之間的面積,即劃斜線的區(qū)域面積Ai.按類似的方式,應用下面的等式可計算區(qū)間[X0,Xn]內曲線ab和X軸之間的面積:A=∫xnx0Ρ(X)dX=∫x1x0Ρ(X)dX+∫x2x1Ρ(X)dX+\:+∫xnxn-1Ρ(X)dX=A1+A2+\:+An-1=n-1∑i=1Ai(24)下面以兩個例子來比較新方法與常規(guī)梯形法、辛普生法和等距三次樣條法的精度.例1我們利用中國科協(xié)北京綠化基地的一小地塊,以4種方法計算其面積來做比較,結果見表1.例2對中國科協(xié)北京綠化基地某地塊,首先用全站儀按小間距(≤3m)精確測量其面積做為真值,而后用DGPS測量,第1次測量12個點,第2次測量24個點,測點選擇在地形變化處,并用皮尺檢測距離,參與秩虧自由網平差,各種算法計算結果見表2、3.由表1～3可知,不管測點的數量是否變化,抑或是不同林地面積量測,三次樣條法量算林地面積的精度都要大大高于其他3種方法.4次樣條法的適用性1)由于圖上量測法無法避免系統(tǒng)誤差帶來的精度損失,因此,實際應用中,