空載光達正高化算探討：以高屏地區為例

(1)

Volume 12, No.1, March 2007, pp. 73-82

空載光達正高化算探討：以高屏地區為例

史天元

¹

何心瑜

²

陳大科

³

摘要

空載光達進行測量作業時，以 GPS 為定位工具，故通常所得點雲之參考坐標系統為 WGS84 系統，其高程為橢球高。產製正高系統之數值高程模型過程中，需利用大地起伏模式將橢球高改算為正高，提供民生使用。本研究使用內政部高屏地區空載光達地形測量之數據，並以一檢核路線視為真值，檢核空載光達所獲得之點雲橢球高誤差，DEM 部分之標準差為 0.172 公尺。再以本研究進行當時之國內五個大地起伏數值模式進行正高化算，分析化算成果，比較檢核五模式。比較結果，若以 Hw2001 為參考，以 Hw2005 成果較為接近，標準差為 0.174 公尺， Hw2003 為 0.199 公尺。此一成果顯示，其主要誤差成分為空載光達誤差，但是由於缺乏實測正高，具體差異尚難確實估算。

關鍵詞：空載光達、正高化算、大地起伏

1. 前言

數值高程模型 (Digital Elevation Model) 簡稱 DEM，其儲存有多種數據結構，最常用之格式是以格網方式呈現地表起伏。格網大小可視實際需要決定，格網越細，越可真實呈現地表情況。另一種常用之格式是不規則三角網 (TIN, Triangulated Irregular Network)，其優點為可使用原始觀測之離散點為節點，無須內插出格網點。台灣地區現有之 DEM 已無法滿足現今各界需求，需要更高精度和解析度之 DEM，而生產高解析度的 DEM 方法有許多種，其中一者為利用空載光達 (Airborne Light Detection And Ranging，LiDAR)進行產製。空載光達所產出者為點雲(Point Cloud)形態，亦即為不規則分布於空間之大量離散點，可用以構成不規則三角網式數據，或再經內插產生格網式數據。

空載光達的點雲資料是由同時觀測之 GPS 資

料與慣性導航系統(Inertial Navigation System, INS) 所得之姿態元素，配合雷射測距與掃描角解算各雷射點之三維坐標，其高程系統通常為 WGS84 (World Geodetic System 1984)系統。

正高(H)是地面上一點沿著與等位面垂直的重力方向到大地水準面的距離。橢球高(h)是地面上一點沿著法線方向到橢球面的距離。同一地面點的橢球高和正高間的差值即為大地起伏(N)，如圖 1 所示。將橢球高轉換為正高的過程稱為正高化算。全球定位系統所提供之三維坐標，其高程基本上係以橢球面起算的橢球高。如進行差分解算，其高程起算點可能受參考站之高程值性質影響，但是所得高差仍為幾何高。因為一般民生應用所需的高程資料是以大地水準面起算的正高。故需要進行正高化算的步驟。本研究以實證之方式，探討由空載光達的點雲資料經以全台灣數值大地起伏模式進行正高化算所可達到之高程精度。

1國立交通大學土木工程學系教授

2國立交通大學土木工程學系碩士

3工業技術研究院能源與資源研究所研究員

收到日期:民國 95 年 08 月 24 日修改日期:民國 96 年 04 月 11 日接受日期:民國 96 年 04 月 12 日

(2)

74

圖 1、正高、橢球高與大地起伏關係圖 (Zilkoski et al., 2005)

2. 大地起伏模式介紹

大地起伏模式之建立，需使用重力、地形等觀測量，為物理大地領域一項重要議題。本研究進行時，台灣地區大地起伏模式有黃金維 2001、黃金維 2003、陳國華 2004、呂誌強 2004 和黃金維 2005(黃金維，2003；陳國華，2004；呂誌強，2004；黃金維，2006 個人聯絡)所提供之五種不同數值大地起伏模式；其中黃金維 2001 和陳國華 2004 模式解析度為 30 秒，黃金維 2003、呂誌強 2004 有 3 秒和 30 秒兩種不同解析度，黃金維 2005 則提供 3 秒解析度。各模式之坐標系統均為 TWD97 坐標系統，

而五者間因方法、DEM 解析度和使用資料不同，

計算出的大地起伏值亦不相同。

五種模式簡述如下：

● 黃金維 2001

使用去除 - 回復技術 (Remove-Restore Technique)配合最小二乘配置法(Least Square Collocation, LSC)，並結合全球大地位模式 EGM96，利用美國國防製圖局編修之全球數值圖圖庫，由其中取數化等高線資料所製成台灣地區 30 秒網格數值地形資料與台灣地區五組重力觀測資料，考慮剩餘地形模型 (residual terrain model, RTM) 理論，計算出 30＂x 30＂黃金維 2001 大地起伏模式。之後

以 Hw2001 表示黃金維 2001 模式。

● 黃金維 2003

本計畫使用最小二乘配置法，以整合後之六組重力及測高資料計算大地起伏。計算之步驟為去除-回復法,此法中所需之大地起伏之長波長及短波長部分(剩餘地形效應)，將分別以 EGM96 及之農林航測所經提供之 40 公尺解析度 DTM 資料，以最小曲率法建構 3 秒網格之 DTM 求得，其精度較 DTM 採用 30 秒解析度的網格資料，且重力資料未包含一等一級重力資料之大地起伏模式為優(黃金維，2003) ，以 Hw2003 代表黃金維 2003 模式。

● 呂誌強 2004

該研究使用方法和黃金維 2003 相同，使用之 DEM 資料為農委會產製之台灣地區數值地形資料、交通大學所處理之衛星測高資料、與台灣地區六組重力觀測資料。除了計算出大地起伏模式，另外用水準點觀測值修正大地起伏面。該研究結論：較高解析度之 DEM 格網對大地起伏模式之精度有提升效果，在平坦地區改善程度有限，在地形起伏劇烈的山區則有較顯著的改善 ( 呂誌強， 2004)。之後以 Lu2004 代表此模式未修正前，

Lu2004v2 則表示修正後之模式。

● 陳國華 2004

利用 TWVD2001 水準與 GPS 網的平差結果，以 LSC 技術結合 1920 個一等水準點的橢球高、正高與重力法大地起伏等資料進行高程控制網整體平差計算，使用兩種不同的平差模式：一是帶有附加參數的「高程值」平差模式；一是為帶有附加參數的「高程差」

平差模式。結果顯示前者較重力法大地起伏改善 74%；後者則改善 69%(陳國華，2004)。

之後以 KH2004 代表此模式。

● 黃金維 2005

使用最小二乘配置法以內政部原有重力資料及內政部委託國立交通大學執行之「空

(3)

載重力測量工作」測算成果整合後之空載重力、船測重力及衛星測高資料，利用最佳化的大地位模式如 EGM96 或由 GRACE 衛星運動所求得之大地位模式 GGM01S 計算求得計算大地起伏與垂線偏差，並利用一等一級水

準點修正大地起伏值。提供台灣地區 3＂x 3＂解析度之大地起伏模式。本研究所使用之版本為已修正者，以 Hw2005 代表此模式。

圖 2 為五種模式中提供 3 秒解析度者之彩色圖與等值線圖。

(a) 黃金維 2003-3 秒解析度 (b) 呂誌強 2004-3 秒解析度

(c) 呂誌強 2004-3 秒解析度(修正後) (d) 黃金維 2005-3 秒解析度(修正後)

圖 2、台灣地區大地起伏模式(3 秒解析度)彩色圖及等值線圖

3. 空載光達資料高程精度分析

3.1 研究數據

本研究所使用之空載光達數據為工研院能資所於民國九十四年間所蒐集，使用之空載光達系統為 Leica ALS50 。測試區為高屏地區 ( 東經 120°22＇30＂~120°43＇30＂，北緯 22°23＇25＂~22°57＇55＂)，其面積約為 2027 平方

公里，平均密度每平方公尺 1.4 點~1.5 點。高屏地區基本資料如圖 3 所示。掃瞄區內覆蓋可區分為五種：裸露地、矮植被(短草、茶區、矮樹)、疏遮蔽林、森林密遮蔽區域、都會區。如圖 4 所示。此五種分類，主要為根據不同程度之地形面遮蔽情況所分。由幾乎無遮蔽之裸露地至高遮蔽率之森林密遮蔽區域，表示不同植披密度。而都會區之主要遮蔽因素為人工建物。不同之遮蔽情況，會造成不同之雷射點雲穿越度，由而形成地面點雲密度之差異。

(4)

76

3a 位置圖 3b 衛星影像

3c 高程變化圖 3d 坡度變化圖

圖 3、高屏地區基本資料(工研院，2006)

圖 4、五種具代表性之遮蔽物分區示意圖(工研院，2006)

(5)

空載光達資料高程部分之比較，以經過人工編修，未格網化之不規則 DEM 及 DSM 比較其精度。

DEM 係用類別為地面點者作比較對象；DSM 採用第一回波(first return)和唯一回波(only return)，刪除距離地形趨勢面顯著較遠之離散點，如浮在空中之點(air point)後進行比較。又各分別比較橢球高及正高。由於 DEM 與檢核點比較部分，會受地物影響可能造成點位對應錯誤，亦即由光達點雲中匹配出對應檢核點之高程時發生錯誤，此一錯誤可歸之於採樣誤差(Sampling Error)，肇因於空載光達點雲雖然密集，但是並非連續。由經驗顯示，如檢核點位於橋樑上，而最近之雷射點在河床，往往因由點雲內插所得之對應高程亦在河床，造成兩者差值過大。此類情況當然應予刪除，DEM 檢核時不得採用位於橋樑上的檢核點，DSM 則不受此影響。此外，由於點雲過濾之過程，往往產生疏化之現象，

DEM 產生對應點之內插誤差較 DSM 之內插誤差較大。基於上述原因，本研究分別比較 DEM 與 DSM 兩者。

3.2 探討方法

因高屏地區地形西部較平坦，越往東向則進入山區，故於測試區內擇一東西向道路做為檢核路線，檢核路線地理位置如圖 5 所示。施測時，先以靜態 GPS 測量引測控制點，再利用 RTK 測得檢核路線上檢核點之三維坐標，此時所獲得之高程為橢球高。正高化算則利用 Hw2001 大地起伏模式進行正高化算求得正高(陳明華，2005、2007 個人聯絡)。

在點雲部分，以檢核路線為中心線外擴約 100 公尺的範圍，分別讀取範圍內的 DEM 和 DSM 資料 ( 仍為橢球高 ) 。以 TerraScan 軟體 (TerraSolid, 2004)，利用不同大地起伏模式化算點雲資料之高程值後，再與檢核點比較，比較方式為讀入檢核點資料(E, N, H)，點雲資料再自行以組 TIN 方式內插出和檢核點(E, N)相同之點的高程資料，輸出比較結果。檢核比對過程可依要求設定參數，本研究採用之參數為 Max triangle：5 m；Max slope：30 degrees；H tolerance：0.15 m，之後以 5-30-0.15 表示

此參數值。Max triangle 為以檢核點為中心的搜尋半徑；Max slope 為檢核點至包含三角形頂點間之最大坡度值，凡大於此值，檢核點與模型內插點間高程差值不予計算，以上兩者因為要求檢核點周圍 5 公尺內需保持平坦，故將參數值設如上述； H tolerance 為正常雷射點高程值之精度，此值使用於前項坡度之計算，以避免小面積三角形造成大坡度，而超越 Max slope。在本研究中，皆採用相同內插方式，故不同內插方法對成果之影響不予討論。但是採用兩組比對參數，成果分列於表 1、2。

圖 5、檢核路線地理分布示意圖

檢核點數總計 241 點，DEM 部分扣除橋樑上之點位 45 點(如圖 6)；DSM 則無影響。下列比較分為橢球高及正高兩類討論。

(6)

78

1. 橢球高

表 1 為參數為 5-30-0.15 之比較成果。根據比較結果得知不論是 DEM 或 DSM，差值過大的點通常位於被遮蔽處。圖 7 展示 DEM 部分差值最大 (-0.793 公尺)的點，DSM 部分差值大於 1 公尺者如圖 8 及圖 9 所示，此外，DEM 部分：其他檢核點與周圍地面點差值，值域在-0.438 至 0.350 公尺間。

造成少數點位差值較大原因可能為檢核點位置靠近道路邊緣等高程變化大之不理想處。DSM 部分：其他檢核點與周圍雷射點差值值域在-0.277 至 0.427 公尺間。

由於使用參數 5-30-0.15，DEM 部分有 2 點坡

度條件不合，調整參數為 5-41-0.3，則無此情況發生；DSM 部分，調整參數為 25-65-0.5 後，坡度條件不合點數亦由 36 點降為 21 點，因角度(Max Slope) 已相當大，故不再調整參數。比較結果見表 2。

2. 正高

正高部分利用前述之提供 3 秒解析度大地起伏模式進行比較，3 秒解析度模式在測試地區之差異，如圖 10 所示。在點雲資料與檢核路線比較過程中，使用參數和橢球高部份相同(5-30-0.15)，表 3 為利用不同大地起伏模式(3 秒解析度)進行正高化算後之正高值和檢核點之正高值比較成果。

表 1、參數為 5-30-0.15 刪除錯誤點後比較成果

表 2、DEM 參數為 5-41-0.3；DSM 參數為 25-65-0.5 比較成果

DEM DSM 最大值(m) 0.335 0.901 最小值(m) -0.793 -12.797 平均值(m) -0.042 0.070 標準差(m) 0.172 1.000 均方根(m) 0.177 1.000 使用點數 196 241 條件不合點數 slope：2 slope：36

DEM DSM 最大值(m) 0.335 1.514 最小值(m) -0.793 -12.797 平均值(m) -0.042 -0.045 標準差(m) 0.173 0.981 均方根(m) 0.177 0.980 使用點數 196 241 條件不合點數 0 slope：21

圖 6、位於橋上不使用之檢核點

圖 7、與周圍地面點高程內插值相差-0.793 公尺之檢核點

(7)

圖 8、DSM 部分：紅色點為檢核點，橘色點為地面點，差值為-12.797 公尺

圖 9、與周圍雷射點內插值相差 0.901 公尺之檢核點

(a)Hw2003-Hw2005 (b)Lu2004-Hw2005 (c)Lu2004v2-Hw2005 圖 10、3 秒解析度之大地起伏模式差異彩色圖(黑色粗線為海岸線)

(8)

80

表 3、利用不同大地起伏模式進行正高化算成果

DEM DSM

單位：公尺 Hw2003- 3sec

Lu2004- 3sec

Lu2004- 3sec-v2

Hw2005- 3sec

Hw2003- 3sec

Lu2004- 3sec

Lu2004- 3sec-v2

Hw2005- 3sec 最大值 0.435 0.365 0.338 0.378 0.916 0.856 0.906 0.936 最小值 -0.773 -0.823 -0.783 -0.753 -0.457 -0.507 -0.477 -0.437 平均值 0.079 0.009 -0.027 0.014 0.144 0.074 0.037 0.079 標準差 0.199 0.194 0.175 0.174 0.206 0.200 0.179 0.177 均方根 0.214 0.193 0.176 0.174 0.251 0.213 0.182 0.194 坡度條件不合點數 2 2 2 2 14 14 14 14

表 3 中可看出四個模式值域略有差異。以均方根而言，除 Hw2003 約為 21 公分，其餘皆小於 20 公分，尤以 Hw2005 值最小。就表 1 和表 3 中 DEM 部分標準差和均方根而言，以經過水準點修正之模式差異較小，尤以 Hw2005 與化算前之值無顯著差異，且與 Lu2004-3sec-v2 相近，可能因為檢核路線未曾經過山區，無法完整顯示兩者相異。但是因為缺乏實測正高，表三僅顯示此四個模型與 Hw2001 之差異。

4. 結語

如同光達點雲之三向座標具有觀測誤差一樣，大地起伏模式亦具有誤差。Hwang et al. (2005) 及 You (2006)便針對本研究所使用數種大地起伏模式之誤差有所討論。隨同種種重力觀測資料之增加，與演算方法之精進，大地起伏模式精度之逐年提高，實為預期。除以嚴謹物理大地之分析方式產出大地起伏模式外，較小面積之作業區，亦可以面擬合之方式進行區域性大地起伏值之求定，如高書屏等(2006)，即以台中市轄區之大地起伏進行討論。

謹就本研究之發現，整理如下述。

(1) 檢核路線之檢核點橢球高為實測，DEM 部分，表 1、2 顯現平均誤差小於 5 公分，均方根差約 17 公分。此與一般認為之空載光達精

度相符。唯誤差值域為(0.335, -0.793)，其發生原因尚須進一步探討，初步認為由檢核點匹配之確定性造成。DSM 表 1、2 顯現均方根差約一公尺，推測因為檢核點於空載光達點雲中產生對應點之匹配錯誤造成，因為 DSM 包含第一回波和唯一回波，僅刪除距離地形趨勢面顯著較遠之離散點，所包含之點數較多，符合設定條件匹配成功者亦較多，而檢核點多在地形面上，DSM 有非地形面之點，故匹配錯誤之機率高。

(2) 正高比對部分，若以 Hw2001 為參考，以 Hw2005 成果較為接近，標準差為 0.174 公尺，

Hw2003 為 0.199 公尺。此一成果顯示，其主要較差成分為空載光達與檢核點間之誤差。由於缺乏實測正高，具體正高誤差尚難確實估算。

(3) 以本研究之探討經驗，顯現以檢核點檢核光達之作業，就目前所有之檢核作業而言，尚有改善必要。建立標準作業程序，以確實檢驗成果，有其迫切性。根據現有數值大地起伏模型如 Hw2005 之成果報告，其精度應已滿足多數地形圖測製之應用需求，而未來仍有精進空間。如保存橢球高之各項數據，未來便有再次改算之機會。對於未來比對，深具意義。

(9)

致謝

本研究進行過程中承蒙工研院劉進金先生、徐偉城先生提供意見，交通大學黃金維教授提供協助，以及兩位匿名評審提供指正，謹此深致謝意。

參考文獻

工研院，2006。辦理LiDAR測區之高精度及高解析度數值地形測繪、資料庫建置與應用推廣工作案總報告初稿，內政部主辦，工業技術研究院能源與環境研究所執行。

呂誌強，2004。DEM解析度對大地起伏模式之影響，國立交通大學土木工程學系碩士論文。

高書屏、甯方璽、王文安，2006。應用不同方法推求台中市區域性大地起伏值之研究，測量工程， 48(4):3-22。

陳國華，2004。整合TWVD2001水準及GPS資料改進台灣區域性大地水準面模式以應用於GPS 高程測量，國立成功大學測量與空間資訊學系博士論文。

黃金維，2003。台灣地區大地起伏模式精度評估期末報告書，內政部。

Hwang, C., Y.S. Hsiao, and H.C. Shih, 2005.

Geodetic and geophysical results from a Taiwan airborne gravity survey: Data reduction and accuracy assessment, submitted to Journal of Geophysical Research (in press).

Terrasolid, 2004. TerraScan user guide (18.11.2004), Terrasolid.

You, R.J., 2006. Local geoid improvement using GPS and leveling data: Case study, Journal of Surveying Engineering-ASCE, 132 (3):

101-107.

Zilkoski, David B., Edward E. Carlson, Curtis L.

Smith, 2005. Guidelines for establishing GPS-Derived orthometric heights(Version 1.4), National Geodetic Survey.

(10)

82

Orthometric height reduction in airborne lidar operation, A Study with Kao-Hsiung and Ping-Tung area

T.Y. Shih

¹

H.Y. Ho

²

T.K. Chen

³

ABSTRACT

In the general process, GPS is used for the positioning in the airborne lidar operation. Therefore, most frequently the height obtained is referenced to WGS84 ellipsoid and the measurements are in ellipsoidal height. In order to produce the digital elevation model in the orthometric height, the geoid undulation correction should be performed. This paper investigates the differences of five currently available digital geoid models of Taiwan area as applied to the airborne lidar point cloud obtained for the Kao-Hsiung and Ping-Tung area. Based on the differences between the result from the elevation obtained in the field and the elevation interpolated from the Lidar point clouds, the ellipsoid height are examined. The standard deviation for the ellipsoid heights is 0.172 m. Then, the orthometric height differences among those obtained from different models are compared. The standard deviation for Hwang-2005 is 0.174 m, and 0.199 m for Hwang-2001. Due to the lack of directly observed orthometric heights, the real error in terms of orthometric height is still difficult to confirm.

Key Words: Airborne lidar, Geoid undulation, Orthometric height

Received Date: Aug. 24, 2006 Revised Date: Apr. 11, 2007 Accepted Date: Apr. 12, 2007 1 Professor, Department of Civil Engineering, National Chiao-Tung University

2 Master, Department of Civil Engineering, National Chiao-Tung University

3 Manager and Researcher, Geo-Information Technology Lab, Resources Technology Division, Industrial Techonolgy Research Institute

空載光達正高化算探討：以高屏地區為例