LIDAR资料建构不规则三角网.docx

地籍測量稿至：93.11.29第廿三卷第4期，第1-24頁，民國93年12月修正：94.01.03接受：94.01.03利用資料蒐評法自動萃取空載光達資料中的屋頂面資訊Automatic Extraction of Building Roofs by Data Snooping from Airborne LIDAR Data邱 式 鴻*Shih-Hong Chio摘 要空載光達資料是大量分布於地表的三維點位資料，經過處理之後，可以應用於高精度的DEM生產；萃取如建物、道路或者樹木等地物，作為景觀規劃或是其他方面之應用；甚至可以用來進行災害評估。由於空載光達資料中隱含著空間中有意義的點、線、面特徵，對於萃取其特徵建構都市區或含建物區的建物模型供後續應用是目前相當重要的研究主題。而由包含大量精確平面特徵的光達資料中，進一步處理或結構化之後萃取有意義的面特徵，甚至是更有意義的地物資訊供後續應用則需要發展不同的演算法。本文即嘗試由包含建物資訊區域範圍內的光達資料中，萃取屬於建物屋頂面資訊供後續處理應用；也就是在假設地表建物屋頂面是由三維平面所建構而成的前提之下，利用資料蒐評法發展可以確實排除不屬於屋頂面上光達資料的演算法，並精確決定建物屋頂面的資訊供建立建物模型之用。關鍵詞：光達資料、資料蒐評、建物重建。AbstractLIDAR data are consisted of a large number of 3-D points on the terrain surface. After processing, LIDAR data can be used for the generation of high accuracy *國立政治大學地政學系助理教授，.tw ，02-29393091 # 51657DEM, for the extraction of terrain object, e.g. buildings, roads and trees, as well as for disaster assessments. Urban LIDAR data contains lots of meaningful 3-D planes for building reconstruction. Therefore many algorithms are developed to extract the meaningful 3-D plane features. Based on data snooping theory, this paper proposed an algorithm to exclude non-roof point data and extract building roof planes from LIDAR data for the application of building reconstruction. The test results show the feasibility of the proposed algorithm.Key Words: LIDAR Data, Data Snooping, Building Reconstruction.一.前言空載光達(LIght Detection And Ranging, LIDAR)資料是利用GPS (Global Position System)全球定位系統和INS (Inertial Navigation System) 慣性導航系統即時定出飛機飛行瞬間的位置和姿態參數，同時利用雷射測距和光學掃描原理，以每秒約數千次的觀測頻率對地面進行掃瞄，可以快速得到大量分布於地表的三維點位資料，又稱為點雲(Point Cloud)資料 Wehr and Lohr, 1999。空載光達資料可以應用於高精度的DEM生產；萃取如建物、道路或者樹木等地物，作為景觀規劃或是其他方面之應用；甚至可以用來進行災害評估Axelsson, 2000; Priestnall et al., 2000; Steinel et al., 2001。而因為空載光達資料是大量分布於地表的三維點位資料，所以空載光達資料中隱含著對建構地表地物有實質助益且有意義的點、線、面特徵，所以如何由空載光達資料中萃取有意義的特徵物，供後續建構地物資訊或其他應用，是目前相當重要的研究主題。尤其是由包含大量精確平面特徵的都市區或建物區光達資料中，進一步處理或結構化之後，以萃取有意義的平面特徵，甚至是萃取更有意義的地物資訊供後續應用則需要發展不同的演算法。由三維光達點雲資料萃取平面資訊的方法很多，如賴志恆2002以及湯凱佩和曾義星2004均以八分樹法將雷射點雲階層式三維網格結構化之後進行最小二乘平面擬合和平面特徵萃取；Lee2002則以Perceptual Organization分析雷射點雲之空間幾何結構，並將點雲資料群集成平面；Mass 1999則將LIDAR資料轉換成網格式的影像資料，然後利用影像分類輔以其他資訊萃取其面狀特徵。空載光達資料雖然是大量分布於地表的三維點資料，但是除了掃描點空間位置坐標和其反射強度之外，在萃取有意義的平面特徵過程中幾乎沒有相關的點位屬性資料可資運用，因此如何由光達點雲資料中進一步處理得到有意義的平面特徵供後續應用就變成相當重要的工作之。本文就是嘗試由包含建物資訊區域範圍內的空載光達資料中萃取平面特徵，並萃取屬於屋頂面的點位資訊供後續處理應用；也就是在假定地表建物由平面所建構而成的前提之下，嘗試由包含建物資訊的區域範圍內的空載光達資料中發展可以確實排除不在屋頂面上光達資料的演算法，並精確地決定建物各屋頂面的資訊達到精確建立建物模型的目標。空載光達資料既然是大量的三維點雲資料，因此萃取建物平面特徵的過程中可以預期除了構成建物主要屋頂面上的點雲資料外，還包括不是構成主要建物屋頂面上的點資料，比如位於女兒牆、水塔、樓梯間或其他附屬建築物上的點資料，甚至是雜訊所構成的點資料，而為了更精確決定平面特徵資訊，萃取過程中都希望可以排除這些與建物主要屋頂面不相關的點資料。因此，在本研究中利用Baarda1967所提出的資料蒐評法(data snooping) 李德仁，袁修孝，2002，期望在萃取平面特徵過程中，將這類的點位資料視為雜訊排除掉。以下各章節就針對本文的研究方法加以詳細介紹。二、研究方法概述本研究在假定某區域範圍內含屋頂面資訊的情況之下，先將光達點雲資料依據不規則三角網(Triangulated Irregular Network, TIN) 建構方式，將區域範圍內的光達資料建構成由諸多的TIN單元表示。組成相同屋頂面的TIN單元其彼此之間的法線向量應該一致，且高程應該在某一合理範圍內。因此由高程在某一範圍內且具有相似法線向量的TIN單元群組進行屋頂平面之擬合。擬合的原理是採最小二乘原理，而擬合平面過程中則利用資料蒐評法(Data Snooping)概念排除不屬於此屋頂面的TIN單元資料。而區域範圍內的屋頂面可能不只一個，因此必須重複地對可能構成屋頂面的TIN單元資料群組進行平面參數的擬合計算。由於資料擬合平面過程中是以TIN單元的法線向量和高程作為挑選的依據，並以光達資料能夠達到的高程精度為條件來建構屋頂面。因此，極有可能地表上同一屋頂面的資訊會被建構成好幾個不同的平面資訊，因此還要將初步建構完成的平面特徵以其較精確的法線向量和高程進行合併的工作，最後並求出建構完成之屋頂面和其他屋頂面的交線並輸出屬於屋頂面的光達點資訊以便供後續應用。整個研究方法的流程如圖一所示，以下兩章節將研究工作分成空載光達資料前處理以及平面擬合兩部份分別在第三節和第四節中詳細說明。圖一：研究流程圖三、空載光達資料前處理3.1 建構不規則三角網(TIN)單元不規則三角網是常用於地形資料呈現的方法之一，它是以不規則連結的三角形平面來表示數值地形資料，並以三角形節點的紀錄結構為基礎，首先將所有節點的名稱及3D坐標紀錄成一表格，依各三角形分別紀錄其節點名稱再結合所有的不規則三角形，形成三角網來表示三度空間的數值地形資料。此種表達方式可以由紀錄三角形的屬性推測其相鄰關係並進行地形資料分析，較能掌握地形資料的變化，有效的紀錄地形資料的起伏並且節省紀錄空間。其缺點是紀錄形式較複雜，且建構TIN過程中容易因為最佳三角形的條件限制，產生與真實地形資料差異甚大的數值地形模型。當然如果已知某些地形特徵線，則相當建構過程中在物空間已知其空載光達資料的地表特徵線位置。就可以用約制性Delaunay三角網的建構方式建構不規則三角網Shewchuk,1996。此方法可讓使用者指定角度和三角形面積的約束，並指定洞(hole)和凹面(concavities)，運用Rupperts Delaunay refinement algorithm產生無小角度的三角網。實際上建構Delaunay三角網或約制性的Delaunay三角網的演算法甚多 Midtb ,1993; Shewchuk, 1996; Su and Drysdale, 1995，均可以應用於建構二維的Delaunay三角網。然而並非所有地形特徵斷線均可事先獲悉，因此針對沒有特徵線之間的點資料建構三角網還是以幾何約束條件為主。上述方法是要求點集合在構成三角網時，形成不規則三角形上的三個點必須符合二維平面上的某些幾何條件。此外還有一種構網方式，要求在三維空間中建構時，點集合中必須符合每四點形成一外接球時，都不會有任一點位於此外接球的内部，此稱三維Delaunay四面體化(Delaunay tetrahedrization) Midtb, 1993。雖然這種構網方法考慮到空間的第三維度高程資料，但所形成的空間三角網仍以幾何條件為要求，所形成的邊也未必能真正反映地表特徵線。而且因為由光達資料找尋地形特徵線仍是一大困難，且建構之後的TIN仍須作後續處理，所以本研究仍以二維Delanuay三角網化方式建構不規則三角網。而二維Delaunay三角網化，由Midtb1993文章中所列之定義Delaunay, 1934如下：即Delaunay三角網是指在一群點集合中，由不重疊的三角形構成二維三角網，且此不重疊的任一三角形的三個點所形成外接圓中，都不會有任何一點位於此外接圓的内部。亦即本研究先以純幾何的觀點從空載光達資料中建構Delanuay三角網，使其形成的三角形在某些幾何條件下在XY平面上認為是最佳的三角形。其建構方式如以下步驟：A.產生外殼凸多邊形(convex hull)，使得全部地形點上任意兩點的連線均在該多邊形內B.由外向內產生Delaunay三角形C.由外殼凸多邊形的任一個邊開始向內逐一尋找第三個光達點，使通過該三點的圓內不包含其他光達點，此稱為Delaunay三角形D.重複此一步驟直到完成所有三角形二維Delaunay三角網建構完成之後，則區域範圍內的光達資料則由無數個三角形所呈現，此三角形在本文中稱作TIN單元。3.2 以TIN單元的相關資訊群組分類屋頂面資訊如前所述，一個TIN單元是由光達資料的三個點坐標所構成，其屬性資訊包括中心坐標、面積、法線向量、以及三個點位的反射強度所計算出來的平均反射強度等等，在建構過程中這些資訊可作為屋頂面的群組約制條件。在本研究中選用面積、法線方向和高程作為群組屋頂面的約制條件。若組成TIN單元(如圖二所示)三點坐標分別為(X1, Y1, Z1)、 (X2, Y2, Z2) 、(X3, Y3, Z3)，則三點坐標分量的平均值則為代表此TIN單元的位置和高程，而面積計算式如下 (3-1)另外，TIN單元法線向量即為向量與向量 的外積，即 (3-2)其中a1 = X3 - X1; a2 = Y3 - Y1; a3 = Z3 - Z1; b1 = X2 - X1; b2 = Y2 - Y1; b3 = Z2 - Z1由於法線向量表示法在分析上較不易看出其特性陳世崇，1999，故將其轉換成法線的方位角和相對於水平面的仰角來表示之，以利於法線向量的分析。如圖二所示。法線的方位角和相對於水平面的仰角計算式如(3-3)、(3-4) (3-3) (3-4)YXZ法線向量(n1,n2,n3)法線向量方位角n法線向量仰角Angn法線向量位置(X,Y,Z)(X1,Y1,Z1)(X2,Y2,Z2)(X3,Y3,Z3)圖二、 TIN單元法線向量示意圖為了避免平面參數計算受TIN單元面積太小之影響，本研究設定面積大於某一臨界值的TIN單元才納入後續處理，如此一來也可減少計算時間。將面積作為挑選TIN單元的第一個條件，其隱含的意義在於太小面積所計算出來的TIN的法線方向可能比較不精確，意即面積之考量在於挑選足夠好的法線資料作為計算之依據。接下來將滿足面積限制條件的TIN單元，依據法線方向和高程依序進行群組分類。也就是由法線方向和高程條件依序分類並萃取可能存在的屋頂面資訊。假定在空間中不考慮各個平面的位置關係，則依序分析平面法線向量方位、仰角和高程的關係，可以得知若平面具有不同法線之方位角、不同仰角和不同高程，則會構成相異的平面，但也有可能相異的平面具有相同法線方位角，但是其法線仰角卻不同，如圖三所示。而存在相同法線方位角和仰角的平面，有可能因為高度之不同而構成相異平面，如圖四所示。相同法線方位角YXZ相異法線仰角圖三：相同法線方位角，但是法線仰角不同的兩相異平面相同法線方位角YXZ相同法線仰角 圖四：相同法線方位角以及法線仰角但高程相異之兩相異平面（平行平面）實做上，本研究是依序根據TIN單元之法線方位角、仰角、和高程與以群組分類，並分析所給定區域範圍內可能存在之屋頂面資訊，亦即將區域範圍內光達資料所組成的TIN單元進行分類處理，提供後續平面擬合及資料蒐評之用。首先，將TIN單元依據法線方位可能分佈的範圍(0360)以每6為一區間進行類別區分，加上考慮法線垂直向上的特例，此時法線方位角將因X與Y方向的分量均為零而無法計算其方位角，所以最後將法線之方位角分成61種。然後將區域範圍內所有的TIN單元資料依據其法線方位角進行分類，接著以類別區間中TIN單元數量最高且超過某一界限值的法線方位角區間作為可能存在平面資訊的判斷依據。而由於未對TIN單元的位置和高程進行約制，因此對於屬於相同法線方位角的TIN單元群組中，仍可能存在對應不同的法線仰角的TIN單元；也就是此區域範圍內可能具有相同法線方位角但是不同法線仰角的相異平面，如圖三所示。因此，為了區分相同法線方位角但是不同法線仰角的相異平面資訊，需要針對相同法線方位角的TIN單元資料群再依據其法線仰角進行類別區分。本研究中將法線仰角可能分佈的範圍(090)同樣以每6為一區間進行分類，再依據法線仰角類別中數量最多且超過某一界限值的類別再進行下一階段的處理。同樣地因為未對TIN單元的位置進行約制考量，具有相同法線方位角和仰角之的平面，也有可能具有不同之高程，如圖四所示。因此，最後挑選一組最有可能對應某一平面資訊的TIN單元資料群組，進行資料蒐評擬合計算時，須先再依據該法線仰角區間內TIN單元資料群內最高和最低高程的資訊進行類別區分，本研究中以1公尺的區間進行高程類別的區分，再以高程類別區間中數量最多且超過某一界限值的類別中進行平面擬合和資料蒐評之計算。整個計算流程如以下流程圖所示。圖五：將TIN單元依序以法線方位角、法線仰角、以及高程進行類別區分和約制之後實施平面擬合計算之流程圖每一次選取可能存在之平面類別進行擬合和資料蒐評時，對類別區間中TIN單元的個數均加以限制，如相同法線方位角區間內的TIN單元個數必須大於15個，相同法線仰角區間內的TIN單元個數必須大於10個，以及相同高程區間內TIN單元個數必須大於8個。才進入該平面之參數擬合計算和資料蒐評。平面擬合和資料蒐評的程序詳見下一節，如此反覆對於可能構成平面之資料群組逐一進行平面參數之求定，直到沒有其他資料群組可以進行平面擬合計算和資料蒐評為止。四、平面擬合與資料蒐評4.1 平面擬合計算 本研究假定構成某一平面的數學函數為 (4-1) 其中(X,Y,Z)為點坐標，a,b,c為表示該平面的參數。 如此一來，區域範圍內屬於同一平面內的TIN單元都應該滿足方程式 (4-2) 其中(Xi,Yi,Zi)為TIN單元的中心坐標。但是眾所皆知點位量測難免有誤差，在本文中僅假設觀測誤差是出現在位置(Xi,Yi)的高程值Zi 上, 也就是假設位置(Xi,Yi)觀測並無誤差，只有高程Zi觀測有誤差，因此觀測方程式寫成如下所示 (4-3) 其中vi為高程改正數。 而以矩陣型式表示則為 (4-4) 各矩陣之元素如下所示 其中XT=a b c T 為待求平面參數之未知數向量。若假設觀測值為不等權，則觀測值的權矩陣P由下列矩陣形態表示。而在本研究中假定觀測值為等權。 因為上述方程組為線性方程組，因此由最小二乘法原理可知其解X=(ATPA)-1ATPL (4-5)4.2 資料蒐評空載雷射掃描所得的光達資料，除了點分布密度高和精度高的特性之外，還具有多重反射(multi return)、是離散點資料，且在都市地區具有豐富的平面資訊的特性賴志恆，2002。如前所述，本研究試著將原始光達資料經由不規則三角網(TIN)化之後，藉由依序分析TIN單元間的法線向量和高程特性群組分類並自動萃取出地表建物的屋頂面資訊。但是由於TIN單元建構的過程中僅僅考慮點位的幾何特性，所以難免有的TIN單元並不是真正對應到地物的屋頂面，也就是相對於某屋頂面而言這些TIN單元是有問題的；同樣地若是屋頂面還有其他建築結構的存在(如女兒牆、煙囪、樓梯間、水塔、堆積物等等)、或者其他因素的影響，使所產生的TIN單元資料是有問題的。無論如何這些錯誤有問題或者是不屬於相同平面的TIN單元在平面擬合過程中都應該予以濾除，因此本文嚐試應用資料搜評概念，將屋頂面上因其他建築結構所形成的非屋頂面點排除，進而決定精確可靠的平面資訊。資料蒐評是一個觀測值偵錯的方法。依據(4-4)和(4-5)式，改正數向量 (4-5)依據誤差傳播理論，由上式得改正數的協因數矩陣(權係數矩陣) (4-6)若以表示觀測值的平差值，則 所以改正數的協因數矩陣可以改寫為。因此改正數向量 ，此式描述平差系統中，改正數由平差中觀測值向量L和矩陣所決定。而改正數與觀測值真誤差l 的關係可以用 表示。由此式可知觀測值的改正數與觀測值真誤差之間的關係，詳細分析見李德仁及袁修孝，2002。其綜合分析得知(1).任一改正數受到所有觀測值真誤差的影響，(2).某一觀測值真誤差對所有改正數均有影響，(3).單一觀測值真誤差反應到對應觀測值改正數之影響量由矩陣中相對應的對角線元素值所決定，意即 。由上述可知，進行觀測值偵測時，單獨由平差之後的改正數判斷是無法真正得到正確的結果，因此進行觀測值偵錯時必須先找到反應在觀測值的真誤差量，而此誤差量可由平差之後的改正數和矩陣中相對應的對角線元素值計算而得。矩陣中第i個對角線元素稱為觀測值li的多餘觀測分量ri，且平差的多餘觀測值r等於所有相對應的對角線元素值的和，意即矩陣的跡。所以 且 所以由，可以看出觀測真誤差反應在改正數的百分比。而由改正數中誤差 ，對於實際特別重要觀測值不相關的情形下，此時Qvv和Pll均為對角線元素，可得或是。假設觀測值最多只存在一個粗差，若已知觀測值的單位權變方和權矩陣為對角矩陣，則可計算出標準化殘差 (4-7) 式中qvii為Qvv矩陣的第i個對角線元素 若觀測值li不存在粗差，wi服從標準化常態分布wi|H0N(0,1)則可通過標準化殘差的統計檢定來判斷觀測值是否存在粗差。也就是選定一信心水平0(通常為0.1%)，由常態分布表得到檢定的臨界值Ka。若wiKa 則觀測值不含粗差；若wiKa則觀測值可能含有粗差。如果單位權變方未知時，則檢定量ti如下式所列，且為符合t分布的檢定量 (4-8) 式中， 同樣地選定一信心水平0，由t分布表得到檢定的臨界值ta。若tita 則觀測值不含粗差；若tita則觀測值可能含有粗差。本研究所使用的資料蒐評方式乃由3.2節中的分類結果中，由最有可能形成一組平面的TIN單元群組資料出發，計算其觀測值的中誤差，由前述的數學模式中可知此時僅將高程值Z當做觀測值，原本LIDAR資料中搜集到位置坐標X、Y理論上應該也有某種程度的誤差量，但在本文的處理模式中是將位置資訊(X,Y)認為是正確無誤，僅有蒐集到的高程資料有所謂觀測誤差或是存在粗差。所以，若此假設正確無誤則共面的TIN單元的資料群組，其平差之後的觀測值中誤差反應的就是光達系統測量的高程精度。也就是觀測之高程精度若超過某一臨界值，就捨棄不進入本節所述共面資料蒐評的程序中。亦即為了更精確求得一組較佳之觀測值進入資料蒐評的程序中，本研究就是以LIDAR資料本身所能達到最佳高程精度的三倍中誤差作為挑選之標準，而依據史天元和彭淼祥2003，一般光達資料點位可達到之高程精度約為15cm。決定一群最佳的共面資料之後，由此資料群所決定的平面資訊，如平面法線的方位、仰角和高程作為篩選其他TIN單元進行資料蒐評之依據。而一維資料蒐評進行方式，可分為三種探測方法，即(i)向後選擇法、(ii)向前選擇法、(iii)組合向前向後選擇法李德仁，袁修孝，2002。本文採用第二種向前選擇法進行一維資料的資料蒐評。就是由一組最小可能的粗差數據組開始平差，根據平差結果來計算那些被懷疑而未參加平差的觀測值改正數。所以其程序如下：當有任一TIN單元納入資料蒐評時，以最小二乘原理重新計算平面參數之外，也計算其高程觀測誤差於常態分布下變量的檢驗量，在此是假定未知單位權變方時，計算如(4-8)式的t分布的檢驗量在給定顯著水平0下，與t分布表得到檢定的臨界值ta相比較從而判斷相對應的觀測值是否包含粗差。若tita 則觀測值不含粗差；若tita則觀測值可能含有粗差。如此反覆進行，對未參予平差的TIN單元進行資料蒐評，直到沒有TIN單元可以加入平差程序為止。4.3 共面合併與交線計算在上述的平面類別區分以及平面擬合計算過程中，極有可能將屬於同一平面之資訊擬合成數個平面。因此，此步驟主要在於針對應該合併處理的平面擬合結果依據平面高程和法線向量的方位角和仰角進行合併的動作。也就是將找到平面的法線向量的方位、仰角和高程差異在某一臨界值者重新合併計算其平面參數。本研究由精度最好的平面開始考量其他平面和此平面的法線方位和仰角相差在6度，高程相差在1m的平面重新取其所有的TIN單元重新平差計算其平面參數；然後再從未合併處理的平面中，再選取精度最佳的平面反覆進行平面參數的重新計算直到無其他平面需要合併為止。最後將光達點資料跟對應該平面高程差值與決定該平面高程精度三倍中誤差之內的光達點位輸出。五、實驗與分析本實驗使用行政院農委會所提供，於2002年4月14日掃描的的空載光達資料，涵蓋區域為新竹地區，約11.7km * 5km的面積範圍，使用的光達掃描儀器為LHALS40、掃瞄時航高約800公尺、視角35度、掃描點密度約1.7點/平方公尺。從實驗資料中選取四組區域範圍內含建物平面資訊的光達資料進行測試。第一組光達資料含一近似水平面的屋頂面資訊，第二及第三組資料內含兩斜面所構成的山形屋頂資訊，第四組資料內則含兩棟山形屋頂資訊。圖六至圖九的左圖就是原始光達資料的俯視圖，而右圖則為利用Surefer 8.0軟體所建構的二維Delaunay三角網，然後加入其點位對應的高程資訊之後所顯示之3D 不規則三角網形狀。在圖中，不同顏色表示不同的高程資訊。而為了顯示資料，圖中原始點位的顯示大小均誇大表示。 圖六、第一組原始光達資料俯視圖(左)和2D Delaunay TIN三角網(右)圖七、第二組原始光達資料俯視圖(左)和2D Delaunay TIN三角網(右)圖八、第三組原始光達資料俯視圖(左)和2D Delaunay TIN三角網(右)圖九、第四組原始光達資料俯視圖(左)和2D Delaunay TIN三角網(右)在利用資料蒐評法進行平面擬合參數求定時，先對TIN單元資料進行前處理，供後續平面參數擬合之用。由於點密度約為每平方公尺1.7點，因此每個點約佔0.58 m2，依此為依據將TIN單元在平面投影面積大於0.5m2的TIN單元才納入後續處理中。如此設定的用意在於將面積過小的TIN單元資料排除在擬合計算之外。接著依據3.2節所述，將TIN單元的法線方位角以每6度進行類別區分，共區分成61類。然後將法線方位角類別區間內TIN單元群組個數大於某一臨界值(比如50個)的所有TIN單元資料群再依據其法線仰角每6度進行類別區分；同樣地將此法線仰角類別區間內TIN單元資料群個數大於某一臨界值(比如30個)的所有TIN單元資料再依據其法線的高程以每一公尺進行平面資訊的類別區分。在本研究中，進行可能平面資訊分類時，依序規定的法線方位角類別區間內TIN單元個數必須超過30，法線仰角類別區間內的TIN單元個數必須超過15個，以及法線高程類別區間內的TIN單元個數必須超過10個才納入平面參數擬合計算的資料蒐評處理中，而第二組和第四組資料在研究過程中由於上述臨界值的設定使得求定之平面太多，因此，經由反覆試驗之後，將臨界值調整如表一所列之後得到較滿意的成果，各組資料前處理的臨界值設定如表一所列。表一：光達資料前處理進行平面分類之各臨界值設定表TIN單元面積平面分類時TIN單元個數的各項臨界值法線方位個數法線仰角個數高程個數第一組0.5301510第二組1503030第三組0.5301510第四組1301510在可能平面資訊類別區分時，TIN單元個數依序滿足法線方位角區間、法線仰角區間、以及高程區間的個數臨界值設定之後，進入平面參數擬合的平差計算中，擬合計算之後以LIDAR資料本身所能達到最佳高程精度的三倍中誤差作為是否進入資料蒐評之依據，而依據史天元和彭淼祥2003，一般光達資料點位可達到之高程精度約為15cm，而本研究以TIN單元資料群組的平差計算中誤差不大於三倍中誤差(0.45m)再放鬆至0.6m，才進入資料蒐評的處理。這一組TIN單元資料群組就是在4.2節中所述資料蒐評法中一組最佳平面擬合的資料來源。緊接著，選定一信心水平0為0.1%之後，將未參加擬合的TIN單元逐一加入資料蒐評的程序中，過程中則是將未參與平差的TIN單元資料的法線方位角和仰角與這一群最佳資料組所計算出的平面法線方位角和仰角差異在5度者才納入資料蒐評的處理中。各組區域資料經過平面擬合和資料蒐評處理之後平面個數分別是1、9、43、6，此時所擬合計算的平面尚未經合併處理。經過資料蒐評之後，由找到的不同平面之間考慮其合併的可能性。本研究由精度最好，也就是擬合之後中誤差最小的平面開始，考量此平面和其他平面的法線方位和仰角相差在6度，且高程相差在1m的平面，重新將屬於這些平面的所有TIN單元資料重新平差計算其平面參數；然後再從尚未合併處理的平面中選取精度最佳的平面，重複上述步驟，反覆進行平面參數的重新計算直到無其他平面需要合併為止。合併處理之後第二組和第三組的平面個數變成2和12。各組資料合併所找到的最後平面相關資料如表二所示。在此必須說明，所有各組的TIN單元資料在平面擬合計算過程中均已經過坐標平移的步驟，也就是各組分別選取各組中第一個TIN單元的中心坐標作為坐標的平移量，以避免擬合之平面參數存在高度相關的情形。表二：為利用本文所發展之演算法擬合之平面參數。資料組平面編號平面方位角平面仰角平面中誤差平面參數A參數A中誤差平面參數B參數B中誤差平面參數C參數C中誤差第一組1335.489.90.0260.0009 0.0007 -0.0019 0.0007 26.32 0.02 第二組129085.70.1530.0705 0.0018 -0.0257 0.0008 12.50 0.02 210984.50.041-0.0911 0.0006 0.0313 0.0003 14.64 0.02 第三組1218.781.50.0470.0930 0.0007 0.1160 0.0008 26.19 0.06 2265.789.90.0790.0013 0.0003 0.0001 0.0005 -3.49 0.03 3202.589.60.2050.0027 0.0004 0.0065 0.0006 -3.26 0.03 4105.9900.056-0.0003 0.0003 0.0001 0.0003 -3.47 0.02 5217.489.90.080.0014 0.0003 0.0019 0.0005 -3.43 0.02 622689.90.0940.0014 0.0003 0.0013 0.0005 -3.45 0.03 7264.289.90.0590.0013 0.0003 0.0001 0.0004 -3.49 0.02 8590.0013 0.0009 0.0144 0.0011 -2.88 0.07 9296.889.90.1030.0011 0.0004 -0.0005 0.0006 -3.51 0.03 10230.5900.0780.0006 0.0003 0.0005 0.0007 -3.49 0.03 113981.40.059-0.0954 0.0013 -0.1178 0.0012 12.67 0.06 12244.5900.0790.0007 0.0007 0.0003 0.0013 -3.46 0.07 第四組137.889.70.143-0.0036 0.0027 -0.0047 0.0021 17.22 0.12 2175.889.60.269-0.0005 0.0046 0.0069 0.0044 15.29 0.38 383.989.80.209-0.0036 0.0025 -0.0004 0.0033 15.80 0.24 4235.489.70.1660.0043 0.0031 0.0030 0.0047 14.44 0.45 5103.488.80.297-0.0211 0.0224 0.0050 0.0062 15.86 0.31 6570.0048 0.0051 0.0046 0.0036 17.31 0.22 由於第二組和第三組資料包含一山形屋頂面資訊，由表二看出將平面的法線仰角不大於88度的平面挑出(即第二組的平面編號1,2；和第三組的平面編號1,11)，其法線方位差約180度，求其交線，則可得到屋頂面和屋脊線的資訊，如圖十一和圖十二所示。此外，由於第四組資料所找到的平面其法線向量的仰角均大於88度，所以雖然第四組資料中有兩山形屋頂的資訊，但卻無法找到其交線。以第二組和第三組找到的屋頂面資訊和第一組和第四組的屋頂面資訊為依據，將各組資料範圍內的光達資料點位其高程和對應此平面的高程差異不大於求定此平面3倍中誤差的點位輸出並顯示在圖十到圖十三中。在圖十至圖十三所找到屋頂面點位資訊中的左圖為輸出屋頂面光達資料的俯視圖，右圖則為視覺上的三維透視圖。緊接著進一步對圖十到圖十三的成果加以分析。圖十、第一組光達資料萃取出的屋頂點位俯視圖(左)和三維透視圖(右)屋脊線屋脊線圖十一、第二組光達資料萃取出的屋頂點位俯視圖(左)和三維透視圖(右)屋脊線屋脊線圖十二、第三組光達資料萃取出的屋頂點位俯視圖(左)和透視圖(右)圖十三、第四組光達資料萃取出的屋頂點位俯視圖(左)和三維透視圖(右)首先由圖十一和圖十二中可看出求到的平面交線，雖然其精度如何仍有待驗證。但是由視覺上仍然可以清楚看出此交線即為屋脊線的位置。接著由圖十到圖十三的成果中可以清楚地看出本研究中所使用的資料蒐評法的確將不屬於屋頂平面特徵的資料排除，且由表二可以看出找到平面特徵的中誤差均小於0.3m，將合併找到的平面特徵資訊搭配適當的法線約制，可以進一步找到屋頂面資訊；但是第四組的資料中則是將原本四個斜屋頂分解成六個屋頂面，其原因是本文所發展的資料蒐評法在一開始進行可能平面特徵類別區分時，僅僅考慮TIN單元法線向量的方位角、仰角以及高程，但並未將TIN單元彼此之間的空間位置關係納入考慮，所以才會發生這樣的狀況，因此如何將TIN單元彼此之間的空間相關性也在資料蒐評過程中考量進去會是將來的工作重點。比如利用遙測的分類方法或是影像處理中的區塊成長法，先獲取可能的平面資訊類別，取代本文僅僅利用TIN單元法線向量和高程的簡易類別區分方法，甚至在類別區分過程中將光達資料的反射強度在過程中加以運用，都將是本研究的後續工作，以避免類似像第四組的成果出現。此外，本研究僅僅將輸出屋頂面光達資料或是求到屋頂平面的擬合參數，因此如何將這些資料進行後續處理或配合其他相關資料進行資料融合，如航測資料的納入，得到更完整的資訊會是將來另一個重要的工作之一。除了上述的現象之外，另一個不爭的事實是計算量的問題，還有就是本研究中臨界值設定的問題。如一開始限制挑選TIN單元的面積；以及決定一組最佳平面擬合資料時，設定TIN單元在同一法線方位區間內的個數、同一法線仰角區間的個數、以及同一高程區間個數的臨界值問題。這部分臨界值會決定一組最佳平面擬合資料的品質，而不論如何改變其值，相信都不如選用考慮空間關係性的影像區塊成長的結果來的有用。另外是將資料納入資料蒐評的臨界值。其實只要最佳的一組結果決定，這部分的臨界值可以給的寬鬆些，雖然會增加計算時間，但是不合適的TIN單元資料還是會在資料蒐評過程中加以排除。合併處理時，也有設定不同平面合併的法線方位角和仰角之臨界值以及高程臨界值的問題。擬合的平面特徵之所以會需要合併也正是因為同一平面被區分成數個平面，假定可以找到合適的區塊成長方法則這部分的問題應該可以減到最低。最後是在萃取山形屋頂資訊時，屋頂面特徵選擇的法線仰角臨界值的問題。本研究中僅以法線仰角為88度作為初步篩選的條件，雖然也得到不錯的結果，不過這部分仍有待進一步研究和發展相關的演算法。六、結論與建議本研究利用空載光達資料所組成的TIN單元法線向量方位角、仰角和高程進行可能平面之類別區分之後，再利用資料蒐評法和最小二乘擬合求定平面參數。由研究中可以清楚看出資料蒐評法的確可以將不屬於求定屋頂面的光達資料排除掉，同樣的經由上述的研究和討論可知本文所發展的方法也存在改善的空間，比如如何利用現有的TIN單元資料或是加入航遙測資料，利用影像分類或分塊的方法將TIN單元的空間關聯性也在資料蒐評的過程中考慮進去以得到更完滿的成果，甚至將處理得到的光達和平面資料進一步處理，或是跟相關資料進行資料整合得到更精確完整的資料都是將來的重要工作。七、致謝 本研究承蒙行政院農業委員會提供新竹地區空載光達資料供本研究實驗之用；以及國立成功大學測量與空間資訊學系博士生 王淼 提供光達資料顯圖程式，使得本研究得以順利完成，特此致謝。八、參考文獻1. 史天元，彭淼祥， 2003，”空載雷射掃描數據精度評估程序之研究”，第二十一屆測量學術及應用研討會論文集，第247-256頁。2. 李德仁，袁修孝，2002，”誤差處理與可靠性理論”，武漢大學出版社，第97-102頁。3. 陳世崇，1999，”以屋頂面法限為約制條件協助立體像對中建物邊緣線之篩選”，國立成功大學測量工程學系碩士論文。4. 湯凱佩，曾義星，2004，”以八分樹三維網格結構組織光達點雲資料並進行平面特徵萃取”，第二十三屆測量學術及應用研討會論文集，pp. 143-150.5. 賴志恆，2002，”雷射掃描點雲資料八分樹結構化之研究”，國立成功大學測量工程學系碩士論文。6. Axelsson, P., 2000, DEM Generation from Laser Scanner Data Using Adaptive Models”, International Archives of Photogrammetry and Remote Sensing, Amsterdam, pp.110-117.7. Baarda, W., 1967, “Statistical Concepts in Geodesy”, Ge