地表面積估算與資源分配的DEM模型

1、地表面積估算與資源分配的DEM模型摘要本文針對三維地形的地表面積的大小以及監(jiān)測器的合理布置以求達到最有效的監(jiān)測范圍問題進行了模型分析，建立了DEM模型和三維空間最佳面積覆蓋模型，考慮該矩形區(qū)域內(nèi)地形表面的實際情況，綜合影響模型準(zhǔn)確程度的一系列因素，對其進行模型精確度的分析。問題一：根據(jù)已知測量點數(shù)據(jù)，利用二維差值法生成格網(wǎng)結(jié)構(gòu)的數(shù)字高程模型，運用DEM三維地形地表面積遞推算法求解地表面積，其精度用殘差總和的大小與回歸分析法進行綜合評估。問題二：對于矩形區(qū)域的無盲區(qū)覆蓋問題，建立了三維空間的隨機最佳覆蓋模型，運用分布式啟發(fā)算法完成整個矩形區(qū)域內(nèi)無盲區(qū)的覆蓋問題。文章最后，對該文所應(yīng)用的模型在其他

2、一些重要領(lǐng)域進行了簡單的介紹，擴展了模型從而加強了它的普遍適用性。關(guān)鍵詞：二維差值法 DEM模型 殘差總和 回歸函數(shù) 隨機最佳模型 1 問題的提出要在長4000米，寬3000米的某矩形區(qū)域內(nèi)進行一項野外作業(yè)，需要布設(shè)一些監(jiān)控器。附件data5B.txt中給出了該區(qū)域的步長為5米的網(wǎng)格節(jié)點對應(yīng)的海拔高度值及等高線圖，其概貌如下圖所示：問題一:試通過數(shù)學(xué)建模方法，估算所給矩形區(qū)域內(nèi)海拔在12米以上部分的地表面積，并分析說明計算精度；問題二：如果要在該矩形區(qū)域內(nèi)布設(shè)某型號的監(jiān)控器，已知在無障礙情況下最遠(yuǎn)可觀測距離為500米，為實現(xiàn)整個矩形區(qū)域的無盲區(qū)全覆蓋監(jiān)視，需研究在矩形區(qū)域內(nèi)設(shè)置多少個監(jiān)視器？其

3、最佳位置在何處？若在該區(qū)設(shè)置28個監(jiān)控器，這些監(jiān)控器應(yīng)該設(shè)置在什么位置？可以覆蓋矩形區(qū)域內(nèi)的多大面積。2 問題的分析2.1問題一分析：已知該矩形區(qū)域的等高線圖及相應(yīng)海拔高度值（見附錄），要對海拔在12米以上的部分地表面積進行計算。在計算實際地形的地表面積時，因其地形是隨機起伏的，很難準(zhǔn)確的運用函數(shù)來描述這種起伏狀況，而用傳統(tǒng)的數(shù)學(xué)方法又不易得到準(zhǔn)確值，故運用一種基于逼近思想的遞歸算法來求解DEM模型，從而算出三維地形的地表面積。對其模型精度的分析評估中，影響其精度因素主要有地形表面的特征、原始數(shù)據(jù)的分布及密度、DEM表面的建模方法以及其表面特性等，而原始數(shù)據(jù)的內(nèi)插法為直接影響因素，因此運用殘差

4、總和的大小進行比較，再者結(jié)合回歸函數(shù)分析與之比較，可知其精度較高，模型較為準(zhǔn)確。2.2問題二分析：要實現(xiàn)對整個矩形區(qū)域內(nèi)無盲區(qū)全覆蓋監(jiān)測，可以通過設(shè)計一個基于位置的算法來找到一條連接網(wǎng)絡(luò)中任意S源節(jié)點和Sink指令中心的路徑，并使該路徑上所有點的最小被觀察能力最大化，轉(zhuǎn)化為無線傳感器網(wǎng)絡(luò)中的“最佳覆蓋問題”，運用分布式啟發(fā)算法進行求解。3 模型假設(shè)1）問題的描述涉及到地表的規(guī)模，當(dāng)規(guī)模發(fā)生變化后，問題性質(zhì)不發(fā)生變化。2）問題的出口條件定為判斷各個小三角形面積是否小于閥值。3）監(jiān)測器在任何情況下均可以達到最大范圍的監(jiān)測效果。4符號系統(tǒng)：表示x軸坐標(biāo);：表示y軸坐標(biāo);：表示z軸坐標(biāo);：三角形1的面

5、積;：三角形2的面積;：三角形面積;:第一點的坐標(biāo);:三角形邊長和的一半;：第二點的坐標(biāo);：第三點的坐標(biāo);：第四點的坐標(biāo);:一點二點之間的距離;:一點三點之間的距離;:二點三點之間的距離;:一點四點之間的距離;:二點四點之間的距離;:三點四點之間的距離;：與兩點之間的距離。S：殘差總和y：線性函數(shù)回歸y：二次函數(shù)回歸y：負(fù)指數(shù)函數(shù)回歸5模型的建立與求解5.1問題一針對問題一，在計算地表面積時，首先圈定一個多邊形，因三維地形隨機起伏，故此多邊形內(nèi)部的三維地形表面也是隨機起伏的。在已知平面內(nèi)建立空間直角坐標(biāo)系O-xyz，將多邊形的各個頂點通過鼠標(biāo)在三維空間中進行三維拾取，然后從目標(biāo)多邊形ABCDE

6、（如圖一）的某一個固定的頂點出發(fā)，依次連接各個頂點，這樣便可把多邊形分解為一個一個的三角形ABC（如圖二），為了能夠更加精確的計算出三角形固定的三維地形的表面積，再對三角形依照地形特征進行拆分，關(guān)鍵是要求出三角形內(nèi)的一個點，這個點能夠代表出這個地形的地表特征，將其與其余的點相連接，得到三個更小的三角形（如圖三），以此類推，直到三角形的面積小于一個閥值或者拆分位于三角形的平面內(nèi)，最后利用海倫公式，可得三角形面積。將逼近思想運用于本題，分割成n*n個目標(biāo)多邊形，把每個多邊形分成三個三角形，運用海倫公式求此三角形的面積，利用循環(huán)解把n*n個目標(biāo)多邊形面積相加可得面積總和。 圖一 目標(biāo)多邊形示意圖 圖

7、二 未分割的三角形示意圖 圖三 分割后的三角形示意圖 求解如下：建立空間直角坐標(biāo)系O-xyz，取四個峰頂為四個頂點，連接四個點可得到一個四邊形，任意連接一對角線，將四邊形分割成兩個三角形（如下圖四）利用海倫公式及逼近的思想，可近似的得到地表面積的大小。圖四 將四邊形分割為兩三角形若把本題表面分解成3000*4000個矩形，把每個矩形分成2個三角形，利用海倫公式和MATLAB得出循環(huán)解求得總面積為.因此地表面積為。對運用DEM模型求得的三維地表面積大小的精確度進行殘差總和和回歸函數(shù)模型進行綜合評估，公式內(nèi)容如下：殘差總和公式S=，其中殘差總和越小，模型的精確度越高。線性函數(shù)回歸y=+x+r二次函

8、數(shù)回歸y=+x+x+r負(fù)指數(shù)函數(shù)回歸y=e+r式中，r為回歸方程的函數(shù)，是判斷回歸函數(shù)擬合效果的重要指標(biāo)，運用MATLAB軟件分別求出線性、二次、負(fù)指數(shù)函數(shù)回歸三個函數(shù)各自的殘差總和，由算法結(jié)果可知：負(fù)指數(shù)函數(shù)的殘差總和最小，因此采用負(fù)指數(shù)函數(shù)回歸作為對DEM模型精確度的分析。經(jīng)MATLAB軟件求得，該模型的計算誤差約為1.02%，精確度相對較高。5.2問題二針對問題二，因為所有節(jié)點的監(jiān)測半徑是相同的，當(dāng)找到一個所有節(jié)點覆蓋半徑的最小值r滿足條件：只要經(jīng)過監(jiān)測器節(jié)點為圓心，r為半徑所確定的多個圓內(nèi)部就可以找到一條連接網(wǎng)絡(luò)中S源節(jié)點與Sink指令中心的路徑。求解過程如下：根據(jù)現(xiàn)有的3維WSN拓?fù)?/p>

9、圖生成化簡但不失地理特性的3維德洛內(nèi)三角圖D（x）（如下圖五）留映射為對應(yīng)的2維德洛內(nèi)三角圖D（x）（如下圖六），頂點集就可以盡可能多的保持原有3維WSN拓?fù)鋱D的基本地理特征。在保持3維WSN拓?fù)鋱D一定級別的地理特征條件下，3維WSN拓?fù)鋱D的覆蓋即可轉(zhuǎn)化為 D（x）的覆蓋問題，故只需保證其構(gòu)成的三角形可被傳感器節(jié)點覆蓋監(jiān)測即可。圖五 3維德洛內(nèi)三角圖圖六 2維德洛內(nèi)三角形示意圖為了實現(xiàn)整個矩形區(qū)域的無盲區(qū)全覆蓋監(jiān)視，在該矩形區(qū)域內(nèi)應(yīng)設(shè)置監(jiān)視器的數(shù)量如下圖所示，圖中所有圓的個數(shù)即為所要求得的監(jiān)測器的數(shù)量，共有30個監(jiān)測器，最有可能的位置便是在圖中所示圓的中心處設(shè)置。圖七 區(qū)域內(nèi)監(jiān)控器個數(shù)示意圖如

10、果在該區(qū)域內(nèi)設(shè)置28個監(jiān)測器，最有可能的設(shè)置位置便是在圖六2維德洛內(nèi)三角形的各個頂點處，那樣便可以得到三維地形的最佳覆蓋面積。6 模型的評價和推廣 6.1優(yōu)點： 1）DEM模型求解面積具有客觀性、實際性，在現(xiàn)實生活中易被采納，具有可操作性實用價值。 2）遞歸算法能夠靈活調(diào)整閥值，更好的控制模型的計算精度。 3）考慮較全面，運用較廣泛，實際性強。6.2改進之處： 1）選擇恰當(dāng)?shù)母叱厅c來確定拆分三角形較為困難。 2）模型是在理想的情況下建立的，并沒有考慮野外的動態(tài)因素，因此模型仍存在著不足，但總的來說，還是具有很好的借鑒意義和指導(dǎo)價值的。6.3模型的推廣應(yīng)用DEM模型已經(jīng)在各領(lǐng)域均有廣泛的應(yīng)用，比

11、如線路勘測設(shè)計、水利建設(shè)工程、軍事工程等等方面，對于自然災(zāi)難而形成的洪水、泥石流淹沒區(qū)的估算，山區(qū)日照分析和GIS中的三維可視化，以及本文涉及的平面等高線圖、立體等高線圖，只要是與地理位置相聯(lián)系的空間定位，均可以用到數(shù)字高程模型。7模型改進本文所建模型是從地勢表面宏觀方向考慮的，從而忽略了較為微小的一些真實地勢概貌，若將插值單元的長度與寬度盡可能的縮小以至趨于無窮小，那么所建模型變更能反映出山勢地貌的真實情況，計算出的地表面積以及監(jiān)測器的合理布置的精確性也會有大大的提高。參考文獻李志林，朱慶.數(shù)字高程模型.武漢：武漢測繪科技大學(xué)出版社，2004孫家廣.計算機圖形學(xué)（第三版）.北京：清華大學(xué)出版

12、社，2003任彥，張思東.線傳感器網(wǎng)絡(luò)三維空間最佳覆蓋路由協(xié)議.電子學(xué)報，2006，2朱衡君.MATLAB語言及實踐教程.北京：清華大學(xué)出版社，北京交通出版社，2005姜啟源.數(shù)學(xué)模型.北京：高等教育出版社，1987 楊曉云，顧利亞. 字高程模型的誤差處理與精度評估.鐵道勘查，2004，5侯飛，韓奎峰.多分辨率DEM的插值殘差計算及精度評價模型.工程勘察2010， 附錄一:( 步長為5m的三維立體圖形的程序)圖為該區(qū)域等高線示意圖clc, cleara=load('data1.txt');amin=min(min(a), amax=max(max(a)x0=0:5:4000;

13、y0=0:5:3000;xx0,yy0=meshgrid(x0,y0);mesh(xx0,yy0,a)v=-49,1,12,44,75,107,171,184:80:357,357:100:614,614;figure, c=contour(x0,y0,a,v); clabel(c)圖為步長為5m的三維立體圖形附錄二：（地行地表面積的算法程序）clc, clear %該程序計算最小面積a=load('data1.txt');a=a'x0=0:5:4000; y0=0:5:3000;pp=csape(x0,y0,a)x=0:4000; y=0:3000;z=fnval(p

14、p,x,y);m,n=size(z);s=0;for i=1:m-1 for j=1:n-1 p1=x(i),y(j),z(i,j); p2=x(i+1),y(j),z(i+1,j); p3=x(i+1),y(j+1),z(i+1,j+1); p4=x(i),y(j+1),z(i,j+1); p12=norm(p1-p2); p23=norm(p3-p2); p13=norm(p3-p1); p14=norm(p4-p1); p34=norm(p4-p3); if all(z(i,j),z(i+1,j),z(i+1,j+1)>=12) z1=(p12+p23+p13)/2; s1=sqr

15、t(z1*(z1-p12)*(z1-p23)*(z1-p13); s=s+s1; end if all(z(i,j),z(i+1,j+1),z(i,j+1)>=12) z2=(p13+p14+p34)/2; s2=sqrt(z2*(z2-p13)*(z2-p14)*(z2-p34); s=s+s2; end endends附錄三：（實現(xiàn)無盲區(qū)全覆蓋監(jiān)視區(qū)域內(nèi)設(shè)置的監(jiān)視器個數(shù)的程序） MATLAB程序如下：clc, clearload pdataa=load(' data1.txt');x0=0:5:4000; y0=0:5:3000;v=-49,1,12,44,75,10

16、7,171,184:80:357,357:100:614,614;c=contour(x0,y0,a,v); clabel(c)hold onplot(P1(1),P1(2),'S')plot(P2(1),P2(2),'D'), plot(P3(1),P3(2),'D'), plot(P4(1),P4(2),'D')b=a' for i=1:801 for j=1:601 txy=(i-1)*5;(j-1)*5;b(i,j); td=norm(txy-Pz1); if td>=995 & td<=1000 plot(i-1)*5,(j-1)*5,'o') end endendfor i=1:801 for j=1:601 txy=(i-1)*5;(j-1)*5;b(i,j); td=norm(txy-Pz2); if td>=495 & td<=500 plot(i-1)*5,(j-1)*5,'o') end endendfor i=1:801 for j=1:601 txy=(i-1)*5;(j-1)*5;b(i,j); td=norm(tx