數據採集是DEM的關鍵問題,研究結果表明,任何一種DEM內插方法,均不能彌補取樣不當所造成的信息損失。數據點太稀會降低DEM的精度;數據點過密,又會增大數據量、處理的工作量和不必要的存儲量。這需要在DEM數據採集之前,按照所需的精度要求確定合理的取樣密度,或者在DEM數據採集過程中根據地形複雜程度動態調整採樣點密度。 由於很多DEM數據來源於地形圖,所以DEM的精度決不會高於原始的地形圖。例如U.S.G.S.用數位化的等高線圖,通過線性插值生產的最精確的DEM的最大均方誤差(RMSE)為等高線間距的一半,最大誤差不大於兩個等高線間距。通常用某種數學擬合曲面生產的DEM,往往存在未知的精度問題,即使是正式出版的地形圖同樣存在某種誤差,所以在生產和使用DEM時應該注意到它的誤差類型。 DEM的數據質量可以參考美國U.S.G.S.的分級標準,共分為三級:第一級,最大絕對垂直誤差50米、最大相對垂直誤差21米,絕大多數7.5分幅產品屬於第一級;第二級DEM數據對誤差進行了平滑和修改處理,數位化等高線插值生產的DEM屬於第二級,最大誤差為兩個等間距,最大均方誤差為半個等間距;第三級DEM數據最大誤差為一個等間距,最大均方誤差為三分之一個等間距
5.1格網DEM套用 5.1.1地形曲面擬合 DEM最基礎的套用是求DEM範圍內任意點的高程,在此基礎上進行地形屬性分析。由於已知有限個格網點的高程,可以利用這些格網點高程擬合一個地形曲面,推求區域內任意點的高程。曲面擬合方法可以看作是一個已知規則格網點數據進行空間插值的特例,距離倒數加權平均方法,克里金插值方法,樣條函式等插值方法均可採用。 5.1.2立體
透視圖從
數字高程模型繪製透視立體圖是DEM的一個極其重要的套用。透視立體圖能更好地反映地形的立體形態,非常直觀。與採用等高線表示地形形態相比有其自身獨特的優點,更接近人們的直觀視覺。特別是隨著計算機圖形處理工作的增強以及螢幕顯示系統的發展,使立體圖形的製作具有更大的靈活性,人們可以根據不同的需要,對於同一個地形形態作各種不同的立體顯示。例如局部放大,改變高程值Z的放大倍率以誇大立體形態;改變視點的位置以便從不同的角度進行觀察,甚至可以使立體圖形轉動,使人們更好地研究地形的空間形態。 從一個空間三維的立體的
數字高程模型到一個平面的二維
透視圖,其本質就是一個透視變換。將“視點”看作為“攝影中心”,可以直接套用共線方程從物點(X,Y,Z)計算“像點”坐標(X,Y)。透視圖中的另一個問題是“消隱”的問題,即處理前景擋後景的問題。 調整視點、視角等各個參數值,就可從不同方位、不同距離繪製形態各不相同的透視圖
製作動畫。計算機速度充分高時,就可實時地產生動畫DTM透視圖。 5.1.3通視分析 通視分析有著廣泛的套用背景。典型的例子是觀察哨所的設定,顯然觀察哨的位置應該設在能監視某一感興趣的區域,視線不能被地形擋住。這就是通視分析中典型的點對區域的通視問題。與此類似的問題還有森林中火災監測點的設定,無線發射塔的設定等。有時還可能對不可見區域進行分析,如低空偵察飛機在飛行時,要儘可能躲避敵方雷達的捕捉,飛行顯然要選擇雷達盲區飛行。通視問題可以分為五類[Lee,J.(1991)]: 1)已知一個或一組觀察點,找出某一地形的可見區域。 2)欲觀察到某一區域的全部地形表面,計算最少觀察點數量。 3)在觀察點數量一定的前提下,計算能獲得的最大觀察區域。 4)以最小代價建造觀察塔,要求全部區域可見。 5)在給定建造代價的前提下,求最大可見區。 根據問題輸出維數的不同,通視可分為點的通視,線的通視和面的通視。點的通視是指計算視點與待判定點之間的可見性問題;線的通視是指已知視點,計算視點的視野問題;區域的通視是指已知視點,計算視點能可視的地形表面區域集合的問題。基於格網DEM模型與基於TIN模型的DEM計算通視的方法差異很大。 圖9-13:通視分析,圖上灰色區域為不可見區域 1)點對點通視 基於格網DEM的通視問題,為了簡化問題,可以將格網點作為計算單位。這樣點對點的通視問題簡化為離散空間直線與某一
地形剖面線的相交問題。(圖9-13) 已知視點V的坐標為(x0,y0,z0),以及P點的坐標(x1,y1,z1)。DEM為二維數組Z[M][N],則V為(m0,n0,Z[m0,n0]),P為(m1,n1,Z[m1,n1])。計算過程如下: (1.1)使用Bresenham直線算法,生成V到P的投影直線點集{x , y},K=||{x , y}||, 並得到直線點集{x , y}對應的高程數據{Z[k], ( k=1,...K-1 )},這樣形成V到P的DEM剖面曲線。 (1.2)以V到P的投影直線為X軸,V的投影點為原點,求出視線在X-Z坐標系的直線方程: (0<KH[k],則V與P不可見,否則可見。 2)點對線通視 點對線的通視,實際上就是求點的視野。應該注意的是,對於視野線之外的任何一個地形表面上的點都是不可見的,但在視野線內的點有可能可見,也可能不可見。基於格網DEM點對線的通視算法如下: (2.1)設P點為一沿著DEM數據邊緣順時針移動的點,與計算點對點的通視相仿,求出視點到P點投影直線上點集{x, y},並求出相應的地形剖面{x, y, Z(x, y)}。 (2.2)計算視點至每個 與Z軸的夾角 : (2.3)求得 。 對應的點就為視點視野線的一個點。 (2.4)移動P點,重複以上過程,直至P點回到初始位置,算法結束。 3)點對區域通視 點對區域的通視算法是點對點算法的擴展。與點到線通視問題相同,P點沿數據邊緣順時針移動。逐點檢查視點至P點的直線上的點是否通視。一個改進的算法思想是,視點到P點的視線遮擋點,最有可能是
地形剖面線上高程最大的點。因此,可以將剖面線上的點按高程值進行排序,按降序依次檢查排序後每個點是否通視,只要有一個點不滿足通視條件,其餘點不再檢查。點對區域的通視實質仍是點對點的通視,只是增加了排序過程。 5.1.4流域特徵地貌提取與地形自動分割 地形因素是影響流域地貌、水文、生物等過程的重要因子,地形屬性的空間分布特徵一直是人們用於描述這些空間過程變化的重要指標。高精度DEM數據和高解析度、高光譜、多周期的遙感影像,為人們定量描述流域空間變化過程提供了日益豐富的數據源,而且人們對流域地貌、水文和生物等過程空間變化機理理解的不斷加深,可以說人類已經進入了一個“空間模擬”的時代。基於DEM數據自動提取流域
地貌特徵和進行流域地形自動分割是進行流域空間模擬的基礎技術。 基於格網DEM自動提取流域特徵地貌和進行地形自動分割技術主要包括兩個方面:1)流域
地貌形態結構定義,定義能反映流域結構的特徵地貌,建立格網DEM對應的
微地貌特徵。2)特徵地貌自動提取和地形自動分割算法。格網DEM數據是一些離散的高程點數據,每個數據本身不能反映實際地表的複雜性。為了從格網DEM數據中得到流域地貌形態結構,必須採用一個清晰的流域地貌結構模型,然後針對該結構模型設計自動提取算法。 1)流域結構定義 可以使用一個具有根的樹狀圖來描述流域結構[Shreve],目前絕大多數算法都沿用這一描述方法。在此結構中主要包括三個部分,即結點集、界線集和匯流區集。如圖9-14所示。 圖9-14:流域結構 (a.內部溝谷段 b. 外部溝谷段 c. 內部匯流區 d. 外部匯流區 e. 溝谷結點 f. 匯流源點 g.分水線段 h. 分水線源點) 其具體內容包括幾個概念: 1)溝谷線段:一條具有兩側匯流區的線段; 2)分水線段:一條具有兩側分水區的線段; 3)溝谷結點:兩條或兩條以上溝谷線的交點; 4)分水線結點:兩條或兩條以上分水線的交點; 5)溝谷源點:溝谷的上游起點; 6)分水線源點:分水線與流域邊界的交點; 7)內部匯流區:匯流區邊界不包含流域部分邊界的匯流區; 8)外部匯流區:匯流區邊界包括部分流域邊界的匯流區。 溝谷結點和溝谷源點共同組成溝谷結點集,所有的溝谷段組成溝谷段集,形成溝谷網路;所有的分水線組成分水線段集,形成分水線網路。溝谷段集和分水線段集共同把流域分割成一個匯流區集。 溝谷段是最小的溝谷單位,溝谷段可以分為內部溝谷段和外部溝谷段。內部溝谷段連線兩個溝谷結點,外部溝谷段連線一個溝谷結點和溝谷源點。同樣,分水線段是最小的分水線單位,也分為內部分水線段和外部分水線段。內部分水線段連線兩個分水線結點,外部分水線段連線一個分水線結點和一個分水線源點。 匯流網路中每一溝谷段都有一個匯流區域,這些區域由分水線集控制。外部溝谷段有一個外部匯流區,內部溝谷段有兩個內部匯流區,分布在內部溝谷段兩側。整個流域被分割成一個個子流域,每個子流域如同樹狀圖上的一片“葉子”。 2)流域特徵地貌自動提取和地形自動分割 特徵地貌定義與提取:根據格線點高程與周圍高程值的關係,將格網點分為坡地、窪地、分水線、谷地、階地和鞍部等幾類。先計算中心點與八鄰點的高程差,然後對高程差進行排序,再根據高程差序列的特性給中心點格網賦一個特徵編碼。然後通過一系列特徵碼的組合特徵,用
模式識別的方法,將格網點劃分到已知的特徵地貌類別。 山脊線和山谷線提取:山脊線和山谷線的自動探測實際上是凹點和凸點的自動搜尋。較為簡單的運算元是2*2的局部運算元。將運算元在DEM數據中滑動,比較每個格網點與行和列上相鄰格網點的高程,標出其中高程最小(探測山谷線)或高程最大(探測山脊線)的格網點。對整個DEM數據計算一遍後,剩下的未標記格網點就是山脊線或山谷線上的格網點。 流域地形自動分割:流域地形自動分割的目標是將整個流域分割成一個個子匯流區。大多數算法是利用3*3視窗計算流向和基於“溢流跟蹤”算法確定匯流網路。算法過程如下: (2.1)格網點流向定義 採用3×3視窗按8方向搜尋計算最大坡向為各格線點的流向。分別為8方向賦不同的代碼,如右圖所示。每個格網有一個從1到9的數值,代表它流向相鄰象元的方向,如該象元為凹點,則其值為5(圖9-15)。 圖9-15:格網水流方向定義 幾種例外情況的處理: A.如果一個格線點的最大坡向格網點與之具有相同的高程值,且之前沒有其它格網點流向這個相鄰格網,則強制流向它。如果還有另外的格網點流向這個相鄰格網,則當前格網點為凹點。 B.當兩個或多個相鄰格網點的最大坡向相等時,先比較各自相鄰格網點坡向,如果仍沒解決,繼續比較相對格網點的坡向,決定賦一個流向。 C.對於具有相同高程值的區域則擴大搜尋視窗半徑,用7×7視窗,如果需要還可以使用更大視窗。 D.在DEM數據的外圍加一圈高程值為0的格網點,強制其最大坡向流向研究區之外。 當所有的格網點處理完畢後,生成一個編碼1—9的流向圖。 (2.2)凹點處理算法 由於凹點的存在,有一些流路不會流向流域出口,而是終止於凹點,所以在進行流域自動分割之前,還要對凹點進行處理。流域中凹點既可能是真實的凹點,也可能是由於插值誤差造成的,所以不能使用簡單的濾波或平滑函式,將凹點全部去除,目的是將凹點造成的斷路連線到主溝谷網路。 搜尋所有凹點的相鄰最低點(有時可能有多個高程相等的最低點),作為凹點的溢出點,以溢出點為起點繼續搜尋比它的高程低或相等的鄰點(已經搜尋的點忽略),判斷是否有比原凹點更低的格網點,如果沒有則以該凹點的溢出點為起點,重複上述搜尋過程;如果搜尋到比原凹點低的格網點,將凹點和最低鄰點的方向倒轉。如圖9-15所示:高程為48的點為一個凹點,搜尋到高程最低的鄰點為49,以它為起點繼續搜尋,找到高程點49,仍比原凹點高程高,則繼續搜尋,又找到另一個高程點49,再找到高程47的點,比原凹點高程低,結束搜尋,按搜尋方向修改流向,如圖9-16中實線箭頭方向所示。 圖9-16:凹點處理 (2.3)提取匯流網路 根據修改後的流向圖,給定一個點,所有流向它的格網點的總和就是該點的匯流區。計算方法是給定一個點,搜尋8鄰點,記錄所有流向它的格網點的位置,然後再以找到的格網點為基點繼續搜尋記錄流向它的格網點,直到沒有新的匯流點為止,所有記錄的格網點構成該點的匯流區。 通常溝谷的匯流區面積大於其它格網點的匯流區面積,可以通過設定一個閾值,將匯流區面積大於此閾值的格網點,標識為溝谷點。很明顯,不同的閾值得到的溝谷網路的複雜性是不同的,這種方法雖然為確定溝谷網路的複雜性提供了靈活性,但也使得溝谷網路的確定具有太大的隨意性。 得到溝谷網路後,可以對溝谷網路進行編碼。首先對溝谷結點編碼。從流域出口開始搜尋遍歷整個匯流網路,對每個溝谷段的上下游結點進行編碼標識,標識值是溝谷段的編碼值,並記錄下這些結點的位置。其次,把溝谷段中的每個格網點標識為溝谷段的編碼值。第三,根據溝谷段上游結點的類型判定溝谷段是內部溝谷段還是外部溝谷段。 (2.4)提取分水網路 遞歸搜尋溝谷段中的每個格網點的匯流區,將匯流區的格網點賦為該溝谷段的標識值,形成各溝谷段的子匯流區。然後進行邊界跟蹤,提取子匯流區的邊界線為分水線,得到分水線網路。最後,對溝谷網路和分水線網路及子匯流區進行拓撲編碼,以完成流域地形的自動分割。 5.1.5 DEM計算地形屬性 由DEM派生的地形屬性數據可以分為單要素屬性和複合屬性二種。前者可由高程數據直接計算得到,如坡度因子,坡向。後者是由幾個單要素屬性按一定關係組合成的複合指標,用於描述某種過程的空間變化,這種組合關係通常是經驗關係,也可以使用簡化的自然過程機理模型。 單要素地形屬性通常可以很容易地使用電腦程式計算得到,包括: 1)坡度、坡向 坡度定義為水平面與局部地表之間的正切值。它包含兩個成分:斜度——高度變化的最大值比率(常稱為坡度);坡向——變化比率最大值的方向。
地貌分析還可能用到二階差分凹率和凸率。比較通用的度量方法是:斜度用百分比度量,坡向按從正北方向起算的角度測量,凸度按單位距離內斜度的度數測量。 坡度和坡向的計算通常使用3*3視窗,視窗在DEM
高程矩陣中連續移動後,完成整幅圖的計算。坡度的計算如下: 坡向計算如下: ( 為了提高計算速度和精度,GIS通常使用二階差分計算坡度和坡向,最簡單的有限二階差分法是按下式計算點i,j在x方向上的斜度: 式中 是格網間距(沿對角線時 應乘以 )。這種方法計算八各方向的斜度,運算速度也快得多。但地面高程得局部誤差將引起嚴重得坡度計算誤差,可以用數字分析方法來得到更好得結果,用數字分析方法計算東西方向得坡度公式如下: 同理可以寫出其它方向的坡度計算公式。 2)面積、體積 (2.1)剖面積 根據工程設計的線路,可計算其與DEM各格網邊交點Pi(Xi,Yi,Zi),則線路剖面積為 其中n為交點數;Di,i+1為Pi與P i+1之距離。同理可計算任意橫斷面及其面積。 (2.2)體積 DEM體積由四稜柱(無特徵的格網)與三稜柱體積進行累加得到,四稜柱體上表面用拋物雙曲面擬合,三稜柱體上表面用斜平面擬合,下表面均為水平面或參考平面,計算公式分別為 其中S3與S4分別是三稜柱與四稜柱的底面積。 根據兩個DEM可計算工程中的挖方、填方及土壤流失量。 3)表面積 對於含有特徵的格網,將其分解成三角形,對於無特徵的格網,可由4個角點的高程取平均即中心點高程,然後將格網分成4個三角形。由每一三角形的三個角點坐標(xi,yi,zi)計算出通過該三個頂點的斜面內三角形的面積,最後累加就得到了實地的表面積。 5.2三角網DEM分析套用 5.2.1三角網內插 在建立TIN後,可以由TIN解求該區域內任意一點的高程。TIN的內插與矩形格網的內插有不同的特點,其用於內插的點的檢索比格線的檢索要複雜。一般情況下僅用線性內插,即三角形三點確定的斜平面作為地表面,因而僅能保證地面連續而不能保證光滑。進行三角網內插,一般要經過以下幾個步驟: 1)格網點的檢索 給定一點的平面坐標P(x,y),要基於TIN內插該點的高程Z,首先要確定點P落在TIN的哪個三角形中。一般的做法是通過計算距離,得到據P點最近的點,設為Q1。然後就要確定P所在的三角形。依次取出Q1為頂點的三角形,判斷P是否位於該三角形內。可利用P是否與該三角形每一頂點均在該頂點所對邊的同側(點的坐標分別代人該邊直線方程所得的值符號相同)加以判斷。若P不在以Q1為頂點的任意一個三角形中,則取離P次最近的格網點,重複上述處理,直至取出P所在的三角形,即檢索到用於內插P點高程的三個格網點。 2)高程內插 若P(x,y)所在的三角形為ΔQ1Q2Q3,三頂點坐標為(x1,y1,z1),(x2,y2,z2)與(x3,y3,z3),則由Q1,Q2與Q3確定的平面方程為 或 令 則P點高程為 5.2.2等高線追蹤 基於TIN繪製等高線直接利用原始觀測數據,避免了DTM內插的精度損失,因而等高線精度較高;對高程註記點附近的較短封閉等高線也能繪製;繪製的等高線分布在採樣區域內而並不要求採樣區域有規則四邊形邊界。而同一高程的等高線只穿過一個三角形最多一次,因而程式設計也較簡單。但是,由於TIN的存貯結構不同,等高線的具體跟蹤算法跟蹤也有所不同。 基於三角形搜尋的等高線繪製算法如下: 對於記錄了三角形表的TIN,按記錄的三角形順序搜尋。其基本過程如下: 1)對給定的等高線高程h,與所有網點高程zi(i=1,2,?,n),進行比較,若zi=h,則將zi加上(或減)一個微小正數ε> 0(如ε=10-4),以使程式設計簡單而又不影響等高線的精度。 2)設立三角形標誌數組,其初始值為零,每一元素與一個三角形對應,凡處理過的三角形將標誌置為1,以後不再處理,直至等高線高程改變。 3)按順序判斷每一個三角形的三邊中的兩條邊是否有等高線穿過。若三角形一邊的兩端點為P1(x1,y1,z1),P2(x2,y2,z2)則 (z1-h)(z2-h)<0表明該邊有等高線點; (z1-h)(z2-h)>0表明該邊無等高線點。 直至搜尋到等高線與網邊的第一個交點,稱該點為搜尋起點,也是當前三角形的等高線進入邊、線性內插該點的平面坐標(x,y): 4)搜尋該等高線在該三角形的離去邊,也就是相鄰三角形的進入邊,並內插其平面坐標。搜尋與內插方法與上面的搜尋起點相同,不同的只是僅對該三角形的另兩邊作處理。 5)進入相鄰三角形,重複第(4)步,直至離去邊沒有相鄰三角形(此時等高線為開曲線)或相鄰三角形即搜尋起點所在的三角形(此時等高線為閉曲線)時為止。 6)對於開曲線,將已搜尋到的等高線點順序倒過來,並回到搜尋起點向另一方向搜尋,直至到達邊界(即離去邊沒有相鄰三角形)。 7)當一條等高線全部跟蹤完後,將其光滑輸出,方法與前面所述矩形格網等高線的繪製相同。然後繼續三角形的搜尋,直至全部三角形處理完,再改變等高線高程,重複以上過程,直到完成全部等高線的繪製為止.