• 沒有找到結果。

應用數位影像相關係數法於現地控制點定位以改善三維量測精度之研究

N/A
N/A
Protected

Academic year: 2021

Share "應用數位影像相關係數法於現地控制點定位以改善三維量測精度之研究"

Copied!
80
0
0

加載中.... (立即查看全文)

全文

(1)

國立高雄大學土木與環境工程學系

碩士論文

應用數位影像相關係數法於現地控制點定位

以改善三維量測精度之研究

A Study on the Accuracy Improvement of

Three-Dimensional Measurement by Using Digital

Image Correlation Method to Position Control Points

研究生:王奕皓 撰

指導教授:童士恒 博士

(2)
(3)

謝誌

終於要畢業了!在這兩年間經歷許多,不論是成功或挫折,都使我成 長許多,甚至可以說是比大學生活更加充實,包含程式設計、各項室內 外實驗以及參加國內、外研討會等,都是令我成長的養分。 回首研究所這兩年間的點滴,首先要感謝我的指導教授童士恒教授, 在研究的過程中提供我許多想法與建議,並與我充分的討論解決問題的 方向與辦法,我才能順利完成論文。除了本科學業,老師亦會要求我學 習並增強其他能力,例如學習Python 語言,或是開會報告時要求使用英 文報告,增進與國外交換生的交流,加強語文能力。而在人生道路上, 老師也總會積極分享各項人生經歷,豐富了我的見聞,也學習了很多解 決問題應具備的思維邏輯,使我獲益匪淺。 其次要感謝要感謝口試委員蕭漢威副教授與林雪淳副教授,對於我的 論文內容與研究方法提供許多寶貴的建議,另外蕭老師及其研究室團隊, 提供了空拍與研究程式方面相關的技術,對於我的研究有相當大的幫 助。 再來要感謝研究室的學弟妹們,在我分身乏術時,總是熱心地提供協 助,幫助我完成實驗,也能給我許多不同的建議。 最後要感謝我的家人,在我的求學生涯中成為我最強大的後盾,使我 能夠專心在學業上,無經濟上的負擔,並給予我精神上最大的支持。 一路走來,深受大家的幫助,謝謝大家!

(4)

I

目錄

目錄 ... I 圖目錄 ... III 表目錄 ... V 中文摘要 ... 1 英文摘要 ... 3 第一章 前言 ... 5 1.1 研究動機 ... 5 1.2 研究目的 ... 7 第二章 文獻回顧 ... 8 2.1 數位影像相關係數法(DIC) ... 8 2.2 特徵點提取 ... 9 2.3 影像匹配 ... 11 2.4 影像拼接 ... 12 2.5 衛星定位 ... 13 2.6 文獻總結 ... 15 第三章 研究方法 ... 16 3.1 尺度不變轉換特徵 ... 16 3.2 數位影像相關係數法 ... 22 3.3 全球導航衛星系統 ... 25

(5)

II 3.5 相機參數 ... 29 3.6 Pix4DMapper ... 32 第四章 實驗驗證 ... 34 4.1 取得空拍影像與控制點之真實三維座標 ... 34 4.1.1 實驗儀器與配置 ... 34 4.1.2 實驗步驟 ... 39 4.1.3 實驗結果 ... 39 4.2 探討觀測時間對 GNSS 定位之影響 ... 42 4.2.1 實驗方法 ... 42 4.2.2 結果與討論 ... 47 4.3 應用 DIC 於控制點定位 ... 51 4.3.1 控制點與檢核點 ... 52 4.3.2 實驗步驟與結果 ... 52 4.3.3 結果與討論 ... 58 第五章 結論與建議 ... 65 5.1 結論 ... 65 5.2 建議 ... 66 參考文獻 ... 67 附錄1 Raspberry Pi 3 作為 IP 分享器之設置 ... 70

(6)

III

圖目錄

圖3.1 將影像依尺度倍率以 2 倍為一單位分組,並縮小影像為原來的一半, 每一組別中進行高斯平滑,再將同一組相鄰的經高斯平滑處理後的影像 兩兩相減可得其DOG 影像(Lowe, 2004) ... 18 圖3.2 中間層的檢測點和它同尺度的 8 個相鄰點以及上下相鄰尺度的 9 × 2 個點,共 26 個點進行比較,以確保在尺度空間和二維圖像空間都 檢測到極值點(Lowe, 2004) ... 19 圖3.3 特徵點描述示意圖(Lowe, 2004) ... 22 圖3.4 物體表面之結構性斑紋與次級影像(網格)示意圖(施等,2006) ... 23 圖3.5 物體表面上之次級影像變形前後相對位置示意圖(施等,2006) ... 24 圖3.6 衛星相對定位示意圖(支秉榮,2004) ... 27 圖3.7 傳統 RTK 架構(曾清涼、儲慶美,1999) ... 29 圖3.8 偏心扭曲之最大切線扭曲軸及最小切線扭曲軸(Weng,1992) ... 31 圖4.1 GNSS 接收儀 ... 35 圖4.2 Raspberry Pi 3 ... 35 圖4.3 BT-290 天線 ... 36 圖4.4 HTC RE 迷你攝錄影機 ... 36 圖4.5 天線腳架 ... 37

(7)

IV 圖4.6 航拍任務之定翼機 ... 37 圖4.7 UAV 飛行路線圖 ... 38 圖4.8 校區空拍全貌與候選控制點分布示意圖 ... 40 圖4.9 同日不同時段之水平較差比較圖 ... 48 圖4.10 同日不同時段之高程較差比較圖 ... 48 圖4.11 同時段不同日之水平較差比較圖 ... 50 圖4.12 同時段不同日之高程較差比較圖 ... 50 圖4.13 控制點與檢核點分佈示意圖 ... 52 圖4.14 x 方向檢核點之誤差標準差 ... 59 圖4.15 y 方向檢核點之誤差標準差 ... 60 圖4.16 z 方向檢核點之誤差標準差 ... 61 圖4.17 x 方向檢核點之均方根誤差 ... 62 圖4.18 y 方向檢核點之均方根誤差 ... 63 圖4.19 z 方向檢核點之均方根誤差 ... 64 附錄圖1.1 ... 70

(8)

V

表目錄

表4.1 參考站之三維座標 ... 40 表4.2 候選控制點之三維座標 ... 41 表4.3 觀測時段表 ... 42 表4.4 點 0 同日不同時段之定位結果,實驗結果編號 A ... 43 表4.5 點 16 同日不同時段之定位結果,實驗結果編號 B ... 43 表4.6 點 0 時段 1 不同日之定位結果,實驗結果編號 C1 ... 44 表4.7 點 0 時段 2 之不同日定位結果,實驗結果編號 C2 ... 44 表4.8 點 0 時段 3 之不同日定位結果,實驗結果編號 C3 ... 45 表4.9 點 16 時段 1 之不同日定位結果,實驗結果編號 D1 ... 45 表4.10 點 16 時段 2 之不同日定位結果,實驗結果編號 B2 ... 46 表4.11 點 16 時段 3 之不同日定位結果,實驗結果編號 D3 ... 46 表4.12 同日不同時段 RTK 定位結果與 3 月 9 日 RTK 定位結果之較差 ... 47 表4.13 同時段不同日 RTK 定位結果之較差 ... 49 表4.14 人工選點標記 2 張影像之檢核點精度 ... 54 表4.15 DIC 選點標記 2 張影像之檢核點精度 ... 54 表4.16 人工選點標記 3 張影像之檢核點精度 ... 55 表4.17 DIC 選點標記 3 張影像之檢核點精度 ... 55 表4.18 人工選點標記 4 張影像之檢核點精度 ... 56 表4.19 DIC 選點標記 4 張影像之檢核點精度 ... 56 表4.20 人工選點標記 5 張影像之檢核點精度 ... 57

(9)

VI

(10)

1

應用數位影像相關係數法於現地控制點定位以改善三維

量測精度之研究

指導教授:童士恒 博士 學生:王奕皓 國立高雄大學土木與環境工程學系 摘要 近年無人飛行載具(UAV)產業蓬勃發展,相關應用亦趨成熟,可搭載影像擷取設 備取得空拍影像。而取得全區地貌常需仰賴影像拼接之技術,若要使接合影像與建立 之三維數值模型有正確的三維座標,則必須要有高精度之控制點。本研究擬以數位影 像相關係數法(DIC)進行控制點定位,以檢核點檢視控制點精度改善之成效。 本研究以單頻全球導航衛星系統GNSS 接收儀,施作即時動態測量 RTK,進行 控制點三維座標量測,因此本研究先探討觀測時間不同對於RTK 之影響。實驗結果 顯示,在同一日內,即使觀測時段不同,定位之水平較差皆在2 公分以內,高程較差 皆在2.5 公分以內,而同時段不同日之定位結果,水平較差大都在 2 公分以內,高程 較差在2.5 公分以內,應用於影像接合所需控制點之三維座標測量有其可靠性。 接著本研究自測量完成的20 個候選控制點中選定 8 個控制點及 7 個檢核點,匯 入影像拼接軟體Pix4Dmapper 中,比較以人工選點方式定出控制點影像座標,與先 用人工選點定出控制點在初始影像上之影像座標後,再以DIC 比對出控制點在其他 影像上之座標兩種方法所得檢核點精度之差異;而本研究同時探討了控制點標記影像

(11)

2

數量對檢核點精度之影響。實驗結果顯示,DIC 確實能改善現地控制點定位精度,然 而在x 方向上與 y 方向上改善幅度相差甚大,推測為 UAV 飛行方向造成的差異,若 能提升影像解析度,使得DIC 在影像有較大變形時仍能精準比對,則有望改善此一 現象。

(12)

3

A Study on the Accuracy Improvement of

Three-Dimensional Measurement by Using Digital Image

Correlation Method to Position Control Points

Advisor: Dr. Shih-Heng Tung Student: Yi-Hao Wang

Department of Civil and Environmental Engineering National University of Kaohsiung

ABSTRACT

In recent years, the unmanned aerial vehicle (UAV) is getting popular, and related applications are also fully developed. It can be used to take aerial photos. In order to obtain the geomorphology of the whole area, it has to rely on the technology of image stitching. High-accuracy control points are important to regulate the coordinate of stitching images and achieve high-accuracy three-dimensional digital model. This study uses the Digital Image Correlation Method (DIC) to position the control points in order to improve the accuracy of the control points.

This study uses Global navigation Satellite System (GNSS) receivers to measure the coordinates of control points. Real-time Kinematic (RTK) technique is adopted. So, this study explores the influence of time on RTK. The result shows that, in the same day, even if the observation period is different, the result of positioning moves less than 2 cm in horizontal, 2.5 cm in vertical. In the same observation period but different day, most of the result of positioning moves less than 2 cm in horizontal, 2.5 cm in vertical. This result shows that using GNSS receivers along with RTK technique to measure the coordinate of

(13)

4

control points is reliable.

Then, 8 control points and 7 check points are chosen from 20 candidate control points in this study. These points are imported into Pix4Dmapper, and the accuracy of check points by two different control point positioning methods is explored. One of the

positioning methods is done by manual. The other determines the control point in the first image by manual, and then using DIC to find the control point on the rest of the images. And the study also explores the influent of the numbers of the images including the control points on the accuracy of check points. The result shows that DIC can improve the

accuracy of control point positioning. However, the improvements in x- and y-direction are different. The flying direction of UAV could be the reason for this phenomenon. If higher resolution images could be used to raise the match precision while using DIC, this

(14)

5

第一章

前言

1.1 研究動機

台灣位於亞熱帶地區,且剛好座落於歐亞大陸板塊與菲律賓海板塊接 觸帶上,常年受到颱風與地震的威脅,當颱風侵襲時,常會伴隨大量雨 水以及強風,降低山坡地的穩定性,大量雨水夾帶砂石開始沿著坡面滑 動產生土石流,且水位急速上升的河流亦會危害鄰近的村落。如1994 年 道格颱風,豪雨造成公路多處坍方,台鐵內灣支線竹東鐵橋遭溪水沖斷; 1999 年,921 地震造成全國性的嚴重災情,更因為強烈的震動使得山坡 地區土石鬆動,震後土石流發生所需雨量僅為震前的四分之一,往後遭 遇豪大雨時更易發生土石流;2001 年,桃芝颱風造成花蓮光復鄉大興村 爆發土石流災情;2009 年,莫拉克風災造成數百人死傷,高雄市甲仙區 小林村更因為遭受到土石流侵害而滅村。 災害發生後,必須取得災區即時資訊,以利評估救災行動與防範後續 災害發生,唯災害發生之初通常伴隨著災區道路交通中斷或豪雨不斷, 往往使得災區即時資訊取得困難重重,這時通常必須仰賴衛星影像或航 拍圖取得災區資訊,但衛星高度過高,使得影像的解析度不佳,且在天 候不良的情況下,地表狀況容易被雲層所遮蔽;而有人機雖然可就近拍 攝,取得解析度較高的影像,但須考量飛行安全,亦受天候狀況影響甚 鉅。

近年無人飛行載具(Unmanned Aerial Vehicle, UAV)產業蓬勃發展,相 關應用亦趨成熟,可搭載影像擷取設備取得空拍影像,成本、解析度皆

(15)

6 優於傳統衛星影像圖,且因不須駕駛員,故在天候不良的情況下,依然 能較無顧慮地執行任務;無人飛行載具亦可飛行至較難到達的區域進行 拍攝,克服災後交通、地理條件不佳的困境。航照圖可經由一系列的影 像處理,取得相關空間資訊,進而提供即時救災、坡地監測與崩塌量計 算等資料。 由於航照圖為低空飛行拍攝,若所需觀測區域範圍較大,單一影像無 法涵蓋欲評估之區域全貌,故須拍攝多張影像,並透過影像接合之技術 完成欲知區域影像之全貌。過去進行影像匹配時,往往採用特徵點匹配 法,如常見的尺度不變特徵轉換(Scale-Invariant Feature Transform, SIFT)、 加速穩健特徵(Speeded-Up Robust Features, SURF)等,精度約略為 pixel 或sub-pixel,然而數位影像係數法(Digital Image Correlation Method, DIC) 在影像匹配方面,精度可達0.01pixel;而在進行後續影像拼接時,現地 控制點的三維座標不易量測,較室內實驗困難,因為室內實驗可透過精 密的高程測微計及水平測微計進行三維座標定位,現地控制點則需透過 大地測量或全球定位系統(Global Positioning System, GPS)定位方能完成。 在過去,僅有GPS 可供衛星定位,其定位精度不佳;近年中國北斗衛星 導航系統、俄羅斯格洛納斯系統及歐盟伽利略系統亦相繼問世,現今衛 星定位儀大多整合了多系統衛星定位,並使用即時動態測量技術

(Real-Time Kinematic, RTK),定位精度可達公分等級,適合作為控制點 真實之三維座標量測方法。

(16)

7

1.2 研究目的

現今航照圖拼接操作多採用套裝軟體,如Photoscan、Pix4Dmapper 與Acute3D 等,大多是利用特徵點擷取與特徵匹配的方式來進行影像拼 接,並利用影像控制點與外方位資訊為輔助。同一個影像控制點需標記 在不同影像上,而傳統方法為透過人工選定控制點,進而校正影像拼接 之三維座標。本研究擬透過人工選定第一張影像之控制點影像座標後, 再以數位影像相關係數法比對出該控制點在其他影像上之座標,並透過 檢核點與單純使用人工選點進行精度比較。

(17)

8

第二章

文獻回顧

2.1 數位影像相關係數法(DIC)

數位影像相關係數法(Digital Image Correlation method,DIC)是以“找 尋演算法”為基礎,使用數位影像相關係數,比對兩張影像的局部相關性, 判定變形前後對應關係。利用有限元素法的概念將變形前後之試體分割 成小網格,比較變形前後網格內灰階值之相關性,在變形後影像中尋找 變形前網格的相對位置。經由變形前後網格的相對位置,即可計算出各 網格之應變量。以下為相關研究: (1) Peters et al.(1982)提出利用計算機分析物件受力前後之數位影像, 求得表面位移量與應變的方法。 (2) Chu et al.(1985)提出以結合變形理論及數位影像相關係數法之量測 技術,及以內插理論來擴展此技術之應用範圍,並以實際案例來驗證 此技術之實用性。 (3) Lu et al.(2000)對於以影像關係法進行應變分析時,使用高次內插函 數模擬位移場,考慮高次位移項之梯度的影響進行探討,實驗結果顯 示在小變形時,忽略二次位移項之梯度影響不大,但忽略二次位移項 之梯度的誤差會隨著變形的增加而愈來愈大。 (4) Vellinga et al.(2000)所發表的研究中提到,由於數位影像相關係數 法之解析度只與像素有關,與實際之長度大小無關,所以結合掃瞄式 電子顯微鏡與數位影像相關係數法,成功的觀測到微小範圍內的應變, 換算成實際空間中的長度,其量測到的位移大約是微米等級,並認為

(18)

9 有可能達到次微米的解析度的發展空間。 (5) Lopez-Anido et al.(2004)提出以兩台高解析度CCD鏡頭試驗性的進 行三維的DIC量測,此法成功判別出複合材料中各材料行為的差異。 (6) Tung et al.(2005)利用DIC量測剛體位移,其x方向與y方向的位移誤 差值為0.004個像素。 (7) Lecompte et al.(2006)的研究,提出了不同大小的分析網格和不同尺 寸的斑點記號之間的關係對DIC分析精度的影響,並建議最佳化網格 大小與記號大小,以提高DIC分析的精度。 (8) Shih et al.(2010)的研究中,利用一台攝影機進行動態的量測,結果 顯示一懸臂樑之模態可利用DIC量測出,並可以此來判斷是否有缺 陷。 (9) Tung et al.(2011)的研究,利用單一一台相機進行三維DIC量測,可 避免因兩臺相機不同的機械與光學性質造成的誤差,經扭曲校正後影 像,物距精度可達0.0043%。

2.2 特徵點提取

圖形在電腦視覺運算中,以灰階矩陣的形式存在,而每個像素的灰階 值會因為影像尺度或是視角的改變而有所變化,因此在影像拼接之前, 必須先找到不受尺度變化及視角改變影響的點,這些具有代表性及獨特 性的點即為特徵點。特徵點主要透過數學演算法來提取,其中較具代表 性的有Harris 角點偵測演算法、SIFT 演算法與 SURF 演算法。此外,亦 有不少針對上述方法所提出的改善方法,以下回顧特徵點提取相關方法:

(19)

10 (1) Harris (1988) 提出 Harris 角點偵測演算法,將影像中的角點做為特徵 點,角點為兩邊線的交點,是一種改變視角時能保持穩定的特徵點, 與其任意方向鄰近點差異性大且具獨特性,為不錯的特徵點。角點的 提取方法主要利用搜索小網格在影像中做小距離滑動,觀察影像強度 變化,如果偵測為角點則網格於任一方向搜索,強度值皆有相當大的 變化。此方法的優點為,偵測效率高且具有較高的重複率,而缺點為, 當影像尺度改變時特徵點位也會隨之改變的情況與運算時間較長的 問題。

(2) Lowe (2004) 改良Harris演算法,提出了SIFT(Scale-Invariant Feature Transform)演算法,改善Harris演算法無法抵抗影像尺度改變的問題。 此方法先對影像做高斯差分( Difference of Gussian, DOG )與影像金 字塔(Image pyramid)分層的方式,使影像特徵點俱備抵抗影像尺度改 變的能力,以及對特徵點做向量化描述來抵抗影像旋轉。此方法優點 為,當影像產生旋轉、尺度與視角改變時的有相當的穩定性。 (3) Bay (2006) 提出了SURF(Speeded Up Robust Feature)演算法,改善

SIFT演算法的特徵點提取法,具有抵抗影像旋轉與影像尺度改變能力 且運算快速的方法。此方法快速的主要原因為,利用積分影像搭配 Hessian matrix探測器找出特徵點,且藉由改變filter大小而非影像大小 來創造尺度空間,並利用X、Y方向的一階Haar wavelet response來描 述特徵點。優點為,運算快速,且特徵點重複性高,抵抗影像視角改 變、尺度變化與模糊能力強。

(20)

11

2.3 影像匹配

透過數學演算法提取特徵點之後,接著便將在不同影像上的特徵點進 行匹配,影像匹配的精度更是直接影響了影像拼接的精度。影像匹配可 以透過比對特徵點的梯度值、距離差及灰階強度等,找出在另一張影像 上相同的特徵點,作為後續拼接使用。以下回顧影像匹配的相關文獻: (1) Lowe (2004) 提出的次近鄰居法(Second-Closest Neighbor),可用於如 SIFT、SURF與BRISK等特徵點演算法之特徵描述子比對,將影像中 的特徵點置入另一張影像中尋找匹配時,利用描述子的向量,計算兩 描述子的最近距離與次近距離,利用這兩個值評估特徵點匹配性。 (2) 謝銘倫 (2002) 提出彩色碼比對法,將影像梯度方向做分類並量化產 生彩色碼做為特徵角點比對資訊。影像比對工作便可經由彩色碼的比 較來達成,若兩點的彩色碼越相似,則表示此兩個角點非常有可能是 互相對應的兩點。透過距離限制和雙向對應一致性的要求,以及與相 鄰點關係的比較,並計算出一評估對應關係的分數,最佳的對應點即 可被選擇。優點為,準確性高且效率高。缺點為,特徵點的選取上必 須避免擷取到空間中的歪斜線在影像上所形成的交點,或是不屬於角 點的特徵點,否則對應時發生錯誤的機率會較高;當場景中有若干移 動過的物體,容易因為運動的方向不一致,在參考對應組合時出現錯 誤。 (3) Rankov (2005) 提出交互關聯匹配法,為一種基於交互關係的匹配方 法,用來找出兩影像間最相似的點。主要利用搜尋視窗在欲偵測之影 像中進行搜索,並利用相關係數來比對出相互匹配的點,此篇研究中

(21)

12 並無探討非整數點的內插,所以精確度相對較低,且無利用特徵點做 匹配,運算量較大,費時較長。

2.4 影像拼接

影像匹配完後,可利用匹配點進行影像間相互關係的運算,利用此關 係將影像進行拼接,拼接後於影像結合處會有明顯的接縫產生,融合影 像演算可用來處理接合縫問題,以及處理模糊帶。以下為影像拼接相關 研究: (1) 王君如 (2014) 提出影像拼接幾何校正法,先利用除錯方法,將錯誤 的特徵點做初步的刪除,再計算均方根誤差(Root Mean Square Error; RMSE)值的方式,挑選出誤差小的特徵點作為幾何校正轉換模型的 參數,將特徵點代入幾何校正模型將單張影像進行坐標轉換,將不同 影像校正至同一坐標系來拼接影像,最後再進行影像重新採樣。從此 文獻中的幾何校正模型發現,無法處理徑向扭曲問題,且文獻中誤差 值計算方法較為繁複,花費時間較多。 (2) Rankov (2005) 提出以梯度演算法為基礎的互補影像融合法,將影像 接合處的地方做調整,利用影像間重合部分的梯度權重來做分配混合, 主要用來消除拼接影像尖端的強度變化。利用此方法混合出的影像, 容易會有雙影影像模糊情形產生,且對不同光照下影像需適當處理光 線補償問題。 (3) BURT (1983) 提出多重帶混合法,製作影像拉普拉斯影像金字塔, 並對金字塔中各層影像做相互融合,並依照權重分配,來處理影像間

(22)

13 接合帶問題。此方法能確實消除影像接合處的模糊帶,且不會有雙影 現象產生。

2.5 衛星定位

建構三維數值地表模型後,將內插出的檢核點高程與實際高程進行比 較,透過標準差來判斷模擬精度。於室內實驗時實際高程可以透過高程 測微計進行量測,室外之現地控制點與檢核點三維座標,則須仰賴大地 測量或GPS 等相關定位技術,為了提高 GPS 定位之精度,可透過 RTK 技術來測量出控制點及檢核點之三維座標,RTK 乃是基於相對定位的概 念發展而來,可達公分等級之精度。以下為衛星定位相關研究: (1) 支秉榮 (2004) 就單一參考站 RTK 定位與 VRS 技術解算之成果加以 比較與分析。由於GPS 原始觀測資料中存在之系統誤差會隨著基線 長度增加而增加,一般要求單主站RTK 基線長度不超過 10 km。在 基線長度30 公里狀況下,使用虛擬參考站(VRS)的成果確實優於 單一參考站RTK 定位,水平精度可達 0.5~2 cm,e 高程經度約 1~4 cm, 也是傳統大地測量和靜態GPS 測量所無法達到的即時性監測。 (2) 沈力洋 (2010) 藉由整天 24 小時真實單頻接收儀觀測量與模擬產生 的 GPS/Galileo 觀測量來驗證單頻 RTK 之效能。透過比較雙系統單 頻、單系統雙頻、單系統單頻之週波值解算效益,顯示雙系統單頻是 最高的;精度方面也有較傳統單系統雙頻近百分之四十五的提升,在 未來使用雙系統之單頻接收儀作為移動站,搭配e-GPS 服務,為最 有效、最經濟之定位方法。 (3) 陳鶴欽等(2005)比較嚴密網形規劃時段觀測靜態測量成果與雙主站

(23)

14 移動站快速靜態成果差異,發現後者作業搭配最小二乘配置法,所得 較差小於3 公分,符合四等控制測量作業之需求。 (4) 沈鍵偉(2004)應用 GPS 分別針對即時動態與靜態進行量測,找出精度 與時間的關係,由結果發現GPS 靜態基線測量在接收 3 小時衛星資 料之精度,適用於地滑監測,並觀察不同日期量測3 小時所得基線長 度,其結果與地表伸縮計有相同位移量。 (5) 陳鶴欽(2009)使用低成本之單頻 L1 GPS 接收儀,結合台灣地區已建 構之e-GPS 即時動態定位系統,並使用虛擬參考站技術。在零基線 之96%解算成果差值小於 1mm,與一般雙頻儀器相當;在一小時靜 態基線測試中,基線長度3 km 以內所得基線固定解之解算率達 100%, 最大較差為2.4 cm,符合一般加密控制測量之精度要求 (30mm±6ppm),7km 範圍內仍有 72%基線可滿足此要求;搭配 e-GPS 系統產生之虛擬參考站,10km 內之靜態基線測量有 75%可求得固定 解,最大較差小於3cm。而 15m 動態測量,可達 90%解算成功率, 水平定位精度約1cm。 (6) Bakula(2013)在施作快速靜態測量及超快速靜態測量時,同一點同時 使用三台多系統衛星定位儀進行控制點定位,每台接收機之間的距離 不超過0.5 公尺,可以在衛星訊號品質不佳的情況下,仍獲得可靠的 定位資訊。 (7) Bakula(2012)為了解決接收衛星訊號時,因為障礙物造成衛星觀測品 質不佳,顯著降低了定位結果的可靠性。於是在施作RTK 測量時, 單一測站同時使用三台衛星定位儀,將其視為重複測量,並將結果重

(24)

15 新初始化,以獲得公分等級之定位座標,且降低了獲得可靠定位資訊 所需要的時間。

2.6 文獻總結

數位影像量測技術是一種新興之技術,由於近年來相機的快速發展, 使用此方法進行量測之精度日益提高,而且應用範圍相當的廣泛,例如 工程界中經常需要材料的應力應變關係。本研究希望將數位影像相關係 數法應用於航拍圖拼接,改善影像拼接精度與接合度。 透過無人飛行載具搭載錄影設備取得航拍圖後,仍需外方位資訊方能 建置三維模型,雖然無人飛行載具上也有搭載GPS,但其精度不高,可 能有數公尺至數十公尺的誤差,因此需要透過高精度的控制點方能校正 三維模型之精度。回顧過去文獻,雙頻GPS 雖然具有高精度的定位能力, 然而成本較單頻GPS 高出許多;而近年俄羅斯的格洛納斯系統、歐洲的 伽利略系統及中國北斗系統亦相繼問世,隨著多系統衛星定位的發展, 可視衛星數量大幅提升,單頻亦可達到公分等級的精度需求。 本研究擬以兩組單頻衛星接收儀,施作RTK 量測控制點及檢核點, 並比較人工選定控制點與DIC 選定控制點之精度差異,探討應用數位影 像相關係數法於現地控制點定位,是否能有效提升影像拼接與三維模型 之精度。

(25)

16

第三章

研究方法

本研究主要目的為利用數位影像相關係數法進行控制點定位,以提高 影像拼接的精度。研究中將分別先以人工選點方式定出第一張影像上控 制點影像座標後,再以數位影像相關係數法比對出在其他影像上的控制 點影像座標,與單純使用人工選點方式定出控制點影像座標,再比較兩 方法所得結果之精度差異。影像拼接採用商用軟體Pix4Dmapper,此軟體 影像匹配法為尺度不變轉換特徵(SIFT),且可以自由匯入連結點供影 像拼接。而控制點在空間中的定位則採用單頻全球導航衛星系統,透過 即時動態測量定出控制點及檢核點之三維座標。

3.1 尺度不變轉換特徵

SIFT(Scale-Invariant Feature Transform) 演 算 法 (Lowe,2004) 改 良 自 Harris演算法,改善Harris演算法無法抵抗影像尺度改變的問題,並具備 抵抗影像旋轉的能力。先對影像做不同尺度的縮小,得到一系列縮小後 的影像,將影像尺寸由大到小排列構築成影像金字塔(Image pyramid),接 著對金字塔中每一層的影像做高斯平滑,減少影像圖像雜訊以及降低細 節層次,增強影像在不同比例大小下的影像效果,如邊緣效果,再將相 鄰的高斯影像相減,做高斯差分(Difference of Gaussian, DOG ),利用DOG 來近似高斯拉普拉斯(Laplace of Gaussian, LOG),LOG可用來檢測影像邊 緣,但運算量較大,所以利用DOG來近似,以減少運算時間。

(26)

17 利用影像金字塔搭配高斯差分,使特徵點具備尺度不變的能力,接 著對高斯差分影像做極值檢測,找出區域極值為特徵點,然後過濾低對 比度特徵點,最後進行特徵點向量化描述,描述特徵點的梯度方位及強 度值。其方法步驟如下: (1) 假設L( x, y, σ )表示不同高斯尺度下的圖像,G(x, y, σ )為不同尺度下的 高斯函數,I( x, y )為輸入的影像灰階值,σ為高斯影像尺度改變參數, Lowe在其研究中進行不同尺度下特徵點的重複性,得出最佳值為1.6, L( x, y, σ )與G(x, y, σ )表示如下: 2 2 2 ( )/2 2

1

( , , )

2

x y

G x y

e

(3.1)

( , , )

( , , )

( , )

L x y

G x y

I x y

(3.2) (2) 將金字塔每一層中相鄰的高斯影像做高斯差分(DOG)(圖 3.1),DOG 可 用來找出影像的邊緣與斑點位置,DOG 的算法如式(3.3),其中 k 為常 數,設定k 21/ss 為影像尺度。 (3.3)

( , , ) ( ( , ,

)

(x, y, ))

,

( , ,

)

( , , )

D x y

G x y k

G

I x y

L x y k

L x y

(27)

18 圖3.1 將影像依尺度倍率以 2 倍為一單位分組,並縮小影像為原來的一半, 每一組別中進行高斯平滑,再將同一組相鄰的經高斯平滑處理後的影像 兩兩相減可得其DOG 影像(Lowe, 2004) (3) 利用同一組內 DOG 金字塔中相鄰兩層圖像間進行比較可以獲得初步 的特徵點,若某像素的 DOG 值與相鄰的 26 個像素相比為極大值或極 小值時,則將此像素視為區域極值所在的位置,也就是初估的特徵點, 如圖 3.2。為了在每組中能檢測 S 個尺度的極值點,則 DOG 金字塔每 組需 S+2 層圖像,而 DOG 金字塔是由高斯金字塔相鄰兩層相減得到, 所以高斯金字塔每組需 S+3 層圖像。

(28)

19 圖3.2 中間層的檢測點和它同尺度的 8 個相鄰點以及上下相鄰尺度的 9 × 2 個點,共 26 個點進行比較,以確保在尺度空間和二維圖像空間都 檢測到極值點(Lowe, 2004) (4) 為了提高特徵點的精度到次像素,並將低對比度或不夠顯著的特徵 點位過濾掉,利用DOG 函數於尺度空間中以初估特徵點為原點的泰 勒展開式(曲線擬合)D x( )直到二次項,計算特徵點的子像素偏移量, 其中x 為到此初估特徵點的偏移量: 2 2

1

( )

2

T T

D

D

D x

D

x

x

x

x

x

(3.4) 將D x( )對x 做一次微分並令它等於零,得到極值ˆx (3.5) 當ˆx 在任一維度的絕對值大於 0.5,表示極值已偏移到鄰近點上,故 會將原本極值點移往該點位,並重新計算子像素位置與差值,反覆迭 代直到收斂。 1 2 2

ˆ

D

D

x

x

x

(29)

20 接著計算D x( )位於ˆx 的值D x( )ˆ ,D x( )ˆ 過小則特徵點容易受雜訊影響, 所以將其剔除,D x( )ˆ 的式子如下: ˆ

1

ˆ

( )

2

T x

D

D x

D

x

(3.6) 文獻中設定D x( )ˆ 門檻值為 0.03,若值小於 0.03 則表示此特徵點的穩 定性不足,將予以刪除。 (5) 因為 DOG 函數對於偵測邊緣上的點相當敏感,但邊緣上的點穩定性 較差,為了刪除位於邊緣上的特徵點,可以利用特徵點的主曲率進 行評斷,對於邊上的特徵點,因為穿過邊緣方向的主曲率會遠大於 沿著邊緣方向上的主曲率,因此穿過邊緣方向與沿著邊緣方向的主 曲率比值遠大於位於角落的比值。為了計算主曲率比值,可以利用 Hessian matrix 的特徵值比例與主曲率比例相同的特性來去除位於邊 緣上的點。Hessian matrix如下: xx xy xy yy

D

D

H

D

D

(3.7) 矩陣的特徵值α 與 β 代表 x 和 y 方向的梯度 2 xx yy xx yy xy

Tr H

D

D

Det H

D D

D

(3.8) Tr H 代表矩陣對角元素之和,Det H 代表矩陣行列式值,令 ,

(30)

21 2 2 2 2 2

1

Tr H

Det H

(3.9) 為了剔除邊緣特徵點,定一門檻值 th ,通常為 10,需滿足(3.10)才 保留 (3.10) (6) 為了使特徵點俱備旋轉不變性,於特徵點描述前,先找出特徵點主 梯度方向,利用圖像局部特徵給每一個特徵點分配一個基準方向做 為特徵點向量,以相鄰像素的梯度方位及強度值找出特徵點的主梯 度方向。m x y, 梯度強度及 x y, 梯度方向分別如式(3.11)及(3.12) 所示。 2 2

( , )

( (

1, )

(

1, ))

( ( ,

1)

( ,

1))

m x y

L x

y

L x

y

L x y

L x y

(3.11) 1

( ,

1)

( ,

1)

( , ) tan (

)

(

1, )

(

1, )

L x y

L x y

x y

L x

y

L x

y

(3.12) 完成梯度的強度與方向計算後,使用直方圖統計鄰域內像素的梯度 和方向。梯度直方圖將 0~360 度的方向依每柱 10 度,分為 36 個柱, 為了減低離特徵點較遠點位的影響程度,每個相鄰像素依據其高斯 圓形權重(強度值大小與梯度方向)加入特徵點直方圖中,最後直方圖 中峰值的方向即為此特徵點的方向,若直方圖統計中,有其他值大 於或等於峰值的 80%時,代表此點可能屬於多方向的情形,此時新 增一個與此點同樣位置同樣尺度,但主要方向為不同方向的特徵點, 通常僅 15%的特徵點被賦予多個方向。 2 2

1

th th

Tr H

Det H

(31)

22 (7) 為特徵點建立一個描述子向量,使特徵點在不同光線與視角下能保 持其不變性,以及使特徵點間可以更容易區別,做為影像間特徵點 匹配依據。以特徵點為中心取16×16 窗口,再將其分為 4×4 子窗口, 計算每個子窗口內 8 個維度(即 0、45、90、135、180、225、270、 315 度)的梯度方向直方圖,最終以 4×4×8 共 128 維的向量來表示。 特徵向量形成後,為了去除光照變化的影響,需要對它們進行正規 化處理,假設得到的描述子向量為H h h1, ,...,2 h128 ,正規化後的特 徵向量為L l l1 2, ,...l128 則 128 1 , 1, 2,3,... i i j j h l i h (3.13) 為了減少非線性亮度的影響,把大於 0.2 的向量值設為 0.2,再進行 一次正規化處理,提高特徵的鑑別性。 圖3.3 特徵點描述子示意圖(Lowe, 2004)

3.2 數位影像相關係數法

數位影像相關係數法(Digital Image Correlation method,DIC)的原理是 利用分析物體表面的灰階分佈特徵,對變形前後的影像進行特徵的比對 來獲得變形前後影像的相對位置,並藉此相對位置計算出影像中任一位

(32)

23 置的位移向量,其它物理量如應變場、剪應變都可以利用此位移向量求 得。 數位影像相關係數法常應用於影像識別技術等相關領域,以物體表面 上的結構性斑紋為影像特徵,結構性斑紋如圖3.4 中散佈的黑色斑點所示, 此特徵則被用來做為分析比對變形前後兩張影像局部相關性的依據。我 們以變形前後兩張影像的局部相關性為目標函數,此相關性影響變形前 後兩張影像的局部對應關係,故以對應函數的參數為最佳化的變數,以 數值方法進行最佳化求解,茲以圖3.5 為例說明其對應關係。圖中變形前 中心點位置在點P ,變形後點 P 位置改變為P*,兩者之間的函數關係如 式(3.14)、(3.15)所示: *

( , )

x

x u x y

(3.14) *

( , )

y

y v x y

(3.15) 其中( , )x y 為變形前影像座標,( *, *)x y 為變形後影像座標,u x y( , )、v x y( , ) 為變形量。 圖3.4 物體表面之結構性斑紋與次級影像(網格)示意圖(Tung et al., 2005)

(33)

24 圖3.5 物體表面上之次級影像變形前後相對位置示意圖(Chu et al., 1985) 對於未變形的影像,可以利用有限元素法的觀念將其分割成數個次 級影像,如圖3.4 中一個方形的格子即代表一個次級影像。假設變形前影 像為影像 A,變形後影像為影像 B,二者之間存在著前述的對應關係, 則我們可以用數位影像相關係數來判定二者間的相關程度,數位影像相 關係數的定義(施等,2006)如下式: 2 2 ij ij ij ij g g COF g g (3.16) 其中gijgij 分別為影像 A 在( , )i j 座標上及影像 B 在( , )i j 座標上的灰階 值,而( , )i j 座標為影像 A 上( , )i j 座標依函數計算後所得在影像 B 上之 對應座標 以最佳化程序求解出每一次級影像最大相關係數來決定對應函數的 最佳參數,即可獲得次級影像變形前後的相對應座標,再運用這些相對 應座標即可計算出個別的位移量,進而獲得位移場與應變場。

(34)

25

3.3 全球導航衛星系統

由全球衛星導航系統國際委員會(International Committee on Global Navigation Satellite Systems, ICG)所認定全球性衛星定位系統共有四支, 分別為美國的全球定位系統GPS、俄羅斯的格洛納斯系統 GLONASS、 中國的北斗導航系統BDS 及歐洲的伽利略系統 Galileo,整合多系統衛星 定位,提供人們更多可視衛星,獲得更可靠的觀測量,可大幅改善定位 精度,而這新一代的衛星導航系統稱為全球導航衛星系統(Global

Navigation Satellite Systems, GNSS),以下介紹四種定位系統: (1) 全球定位系統 GPS GPS 為美國國防部研製之中距離圓形軌道衛星導航系統,1994 年建置 完成初期,在訊號中仍有刻意加入的選擇性誤差,即SA 政策(Selective Availability),此乃為了防範軍事攻擊而降低其精度至 100 米左右,至 2000 年才取消對民用訊號的干擾,精度可達 10 米。GPS 總共由 24 顆 衛星組成,平均分布於6 個軌道面上,在任何時刻至少可以觀測到 4 顆衛星(沈鍵偉,2004)。 (2) 格洛納斯系統 GLNOSS GLNOSS 由蘇聯於 1982 年開始研發,然而一度因為經費不足,造成 衛星設備老化,曾只剩 6 顆衛星運行,目前由俄羅斯維護與運作。另 外 GLONASS 所使用衛星編碼調制技術與 GPS 不同,增加 GPS 硬體 製造廠成本,無法打開民用市場,而使 GLONASS 較 GPS 處於劣勢(陳 鶴欽,2009)。至 2018 年 3 月為止,總計有 24 顆衛星運作中。 (3) 北斗衛星導航系統 BDS

(35)

26 BDS 為中華人民共和國獨立架設之衛星導航系統,主要可分成兩個部 分,第一代為2000 年起運作的實驗系統,第二代為建設中的全球導航 系統。其中第一代官方名稱為北斗衛星導航試驗系統,又稱作北斗一 號,主要於中國境內服務。第二部分官方名稱為北斗衛星導航系統, 又稱作北斗二號。北斗二後於2012 年始提供亞太地區用戶定位服務, 未來將有35 顆衛星運行,並預計於 2020 年完成對全球的覆蓋(楊開元, 2014)。 (4) 伽利略定位系統 Galileo Galileo 為興建中的衛星定位系統,由歐盟通過歐洲太空總署和歐洲導 航衛星系統管理局建造。有別於GPS 的軍事用途,Galileo 起初便以 民用為目的進行開發(陳鶴欽,2009),其為歐盟國家自主的高精度定 位系統,Galileo 於 2005 年發射第一顆衛星,並預計於 2020 年發射完 成30 顆衛星,其中 24 顆為工作衛星,6 顆備用衛星。 GNSS 衛星定位是由位在地面的接收儀接收來自衛星所發射的無線 電磁波訊號,以時間差求得待測點位與衛星之距離,並定出點位之三維 座標。GNSS 衛星提供了兩種觀測量,電碼虛擬距離及載波相位觀測量, 而不論哪種觀測量皆伴隨各種系統誤差。誤差包含衛星軌道誤差、衛星 時錶誤差、對流層誤差、電離層誤差、多路徑效應、接收儀時錶誤差及 天線相位中心誤差(沈力洋, 2010),透過衛星相對定位,可以有效削弱 上述系統誤差,以下將介紹衛星相對定位之相關理論。

(36)

27 圖3.6 衛星相對定位示意圖(支秉榮,2004) 如圖3.6 所示,參考站 A 之位置向量為已知,以𝑅𝐴表示,兩站間基 線向量為∆𝑅𝐴𝐵,則待測站B 之位置向量可表示成: B A AB R R R (3.17) 其中: i i i A A A R R e r (3.18) i i i B B B R R e r (3.19) (3.20) 𝑅𝑖為各衛星之位置向量,𝑒𝐴𝑖、𝑒𝐵𝑖為各衛星之單位向量,𝑟𝐴𝑖為參考站A 與各衛星之距離,𝑟𝐵𝑖為待測站 B 與各衛星之距離,i 為衛星編號。 B A AB i i i i AB B A B A AB A A B B B A AB X X X R R R Y Y Y e r e r Z Z Z

(37)

28

相對定位要求參考站及待測站需有相同時段之觀測時間,藉由空中 及地面的一次差組成二次差分計算,以差分過程來抵銷大部分的系統誤 差,以提升基線求解之精度(曾清涼、儲慶美,1999)。

3.4 即時動態定位(Real-Time kinematic Position, RTK)

傳統的單一參考站 RTK 系統是由參考站、移動站、通訊設備所組 成,如圖3.7 所示。將參考站設在一個已知且固定的點位上和移動站做同 步觀測,透過無線電通訊設備將接收到的載波相位觀測資料傳送到移動 站,讓位於移動站的計算機能夠進行二次差分計算。為了可以達到即時 性的求解需求,因此移動站會在移動的狀態下求解週波未定值,此求解 法稱為 OTF(On-The-Fly),顧名思義就是接收儀在移動的狀態下僅需短 時間的觀測量資料即時求解未定值,一旦週波未定值解算出來後,就可 以很快的求出待測點的坐標位置(支秉榮,2004)。

(38)

29 圖3.7 傳統 RTK 架構(曾清涼、儲慶美,1999) 雖然大部分之系統誤差能在相對定位的概念中,能藉由二次差分消 除,其中電離層與對流層造成之系統誤差仍存在且影響頗大,隨著基線 長度增加而增加,這將導致定位精度降低,甚至無法正確求解週波值, 因此單一參考站 RTK 之施測範圍有相當的限制,一般基線長度不會超 過 10 km。(楊名,1997)。

3.5 相機參數

相機參數包含了內方位與外方位,內方位參數是相機與鏡頭內部的 機組零件所產生,例如焦距、像主點偏移量與扭曲參數,其中扭曲參數 是用來將因鏡頭上的透鏡造成的影像扭曲變形座標轉換成未扭曲座標的

(39)

30 參數;外方位參數則是當相機進行拍攝時所在的空間座標以及其對於 XYZ軸的旋轉量。 1. 內方位參數:指相機內部之參數,包含了像主距 f 、透視中心到像主 點之偏移量X Y0, 0 與鏡頭所產生的扭曲參數。 在扭曲校正的實驗中,通常需要求得扭曲係數以還原影像,而扭曲係 數會因扭曲方程式定義不同而異,在此引用Weng等1992年提出之扭曲 方程式,表示如下: ' ( , ) x x x x y (3.21) ' ( , y) y y y x (3.22) 其中( , )x y 為理想影像座標(無扭曲),( , )x y' ' 為拍攝影像座標(有扭曲), ( , )x y 為 x 、 y 方向之鏡頭總扭曲變形量。 (1) 徑向扭曲(Radial distortion) 鏡頭扭曲絕大部分是因為徑向扭曲的關係,其所造成之扭曲會使 得影像沿著由中心點往外發散的徑向直線變形,令座標點拉近中 心或是遠離中心。 徑向扭曲是徑向距離的函數,在任何角度都是一樣,所以徑向扭 曲也是一個對稱的扭曲,也叫做輻射扭曲。扭曲變形的程度為座 標點到中心距離的函數,其方程式如下: p x x x (3.23) p y y y (3.24) 2 2 2 r x y (3.25) 2 4 6 1 2 3 ... xr k r k r k r x (3.26)

(40)

31 2 4 6 1 2 3 ... yr k r k r k r y (3.27) xryr為徑向扭曲變形量,

r

為到扭曲中心的距離,

x

pyp為 像主點座標k 、1 k2k 為徑向扭曲參數。 3 (2) 偏心扭曲(Decentering distortion) 由於相機中的透鏡並未完全的共線,使得實際的光學系統受到不 同程度的偏心。造成像點位移的方向垂直於徑向線,造成所有的 通過中心點的徑向線彎曲變形所以此扭曲包含徑向扭曲與切向 扭曲。唯有某一個方向上的徑向線沒有變形的情況,這條線即被 稱為最小切線向扭曲軸,而與其垂直的則為最大切線向扭曲軸, 如圖3.8 所示,偏心扭曲之變形量如下: 2 2 1 2 2 2 xp p r x p xy (3.28) 2 2 2( 2 ) 2 1 yp p r y p xy (3.29) xpyp為偏心扭曲量, p1、

p

2為偏心扭曲參數。 圖3.8 偏心扭曲之最大切線扭曲軸及最小切線扭曲軸(Weng,1992)

(41)

32

(3) 薄菱鏡扭曲(Thin prism distortion)

因為非完美的鏡片設計或者製造時使得鏡片或感光元件些微傾 斜。薄菱鏡扭曲會造成類似偏心扭曲的影像變形,也包含了徑向 扭曲以及切線向扭曲。 其偏心扭曲之變形量如下: 2 2 1 xp s x y (3.30) 2 2 2 vp s x y (3.31) xpyp為偏心扭曲量,

s

1、

s

2為薄菱鏡扭曲參數。 2. 外方位參數:是指相機在地面座標內之位置與方位,相機相對於 XYZ 方向之旋轉角 、 、k與真實座標XLYLZL。可利用 GPS 定位、 測量或是共線式等方式去求得。

3.6 Pix4DMapper

Pix4Dmapper 為一套專業航測影像處裡軟體,能快速、自動化處理任 何相機拍攝之航空影像、傾斜影像、地面影像、支援各種波段,能達公 分級的三維精度,另外能迅速產製正射影像、數值地形模型、外方位參 數等成果。 軟體便利且自由度高,使用者可以依需求匯入不同資料,其中一項功 為可以匯入手動連結點輔助影像拼接操作,這是其他軟體沒有的功能, 利用這項功能可以將數位影像相關係數法與特徵匹配法所演算出的點匯

(42)

33

入軟體中,觀察匯入前後拼接影像與三維模型,比較拼接準確度與精度 做為評估的依據,此特點為本研究選用這套軟體最主要的原因。

(43)

34

第四章 實驗驗證

本章將分為三部分進行論述,首先將詳述取得空拍影像與量測控制 點三維座標之實驗方法;其次將探討以單頻GNSS 施作 RTK,量測得出 之控制點三維座標,受衛星觀測時間不同之影響,其定位結果應用於影 像接合之可靠性;最後探究將DIC 應用於現地控制點定位時,是否能有 效提升影像拼接與三維數值模型之精度。

4.1 取得空拍影像與控制點之真實三維座標

控制點之精度乃是影像拼接及三維數值模型校正之重要指標,本研 究之實驗場址為國立高雄大學校區,並以定翼機搭載攝影機飛行,取得 共計1619 張空拍影像,匯入 Pix4Dmapper 進行影像拼接。接著選定特徵 較明顯的地物標記做為控制點,並採用RTK 技術量測控制點三維座標, 共量得20 個點。以下將詳細介紹本節實驗器材、實驗步驟與實驗結果。 4.1.1 實驗儀器與配置 1. GNSS 接收儀:本研究使用美國 Trimble BD910 單頻 GNSS 模塊,在 基線長度不超過5 公里的情況下,其水平精度可達 8mm,垂直精度 可達15mm,搭配開發板自行組建為 GNSS 接收儀,如圖 4.1 所示。 本研究一共使用兩台GNSS 接收儀,一台作為參考站,一台作為移動 站。

(44)

35

圖4.1 GNSS 接收儀

2. Raspberry Pi 3:由於 BD910 模塊並無支援儲存功能,故透過 RS232 轉USB 將 NMEA 格式之衛星定位資料傳入 Raspberry Pi3 儲存,儲存 程式以python2.7 撰寫,將 NMEA 資料儲存為 txt 文字檔;另外可以 Raspberry Pi 3 連接手機網路,並以有線網路方式分享網路予 BD910, 使BD910 具備網路傳輸功能,以施作 RTK。如圖 4.2 所示,相關設 定方法請參照附錄1。 圖4.2 Raspberry Pi 3 3. 天線:本研究使用中國北天通訊之 BT-290,可以接收多系統衛星訊 號,並且可支援雙頻訊號,如圖4.3 所示。

(45)

36 圖 4.3 BT-290 天線 4. 攝影機:採用 HTC RE 迷你攝錄影機,如圖 4.4 所示。錄影後擷取 影像取得航照圖,影像畫素為1920×1080。 圖 4.4 HTC RE 迷你攝錄影機 5. 腳架:天線用之雙腳支架,固定於最低高度 137 cm,地面點位至天 線相位中心總計142.1 cm,如圖 4.5 所示。

(46)

37 圖 4.5 天線腳架 6. UAV:採用國立高雄大學資管系蕭漢威副教授及其研究團隊自行組裝 之定翼機,翼展1400 mm,機身長 925 mm,6 吋螺旋槳,使用 ArduPlane V3.8.4 飛行控制系統,搭載 RE 迷你攝錄影機,總重約 1300g,如圖 4.6 所示。 圖4.6 航拍任務之定翼機

(47)

38 7. 飛行路線:以單 S 型東西向飛行,將校區分為南、北兩部分,分為兩 次航拍任務, 如圖 4.7 所示。 圖4.7 UAV 飛行路線圖 第一次任務如黃線所示,自黑色標點處飛行至紅色標點處,以單 S 型航線東西向飛行,最後回到黑色標點處;第二次任務如橘線所示,達 飛行高度後即開始以單S 型東西向飛行錄影,至黃色標點處後回到黑色 標點處。

(48)

39 4.1.2 實驗步驟 1. 取得國立高雄大學校區之空拍影像:以定翼機搭載 RE 攝錄影機飛行 於150 公尺進行錄影,篩選錄影截圖後共計取得 1619 張航拍影像。 2. 將影像匯入 Pix4Dmapper 進行影像拼接與建立三維模型。 3. 從初步完成的拼接圖上尋找特徵分布較明顯的地物標記作為候選控 制點。 4. 選定工學院六樓頂空曠處作為參考站點為,施作連續 48 小時之單點 定位,並取得平均三維座標。 5. 依序以 RTK 量得 20 個候選控制點之三維座標。 4.1.3 實驗結果 校區影像拼接完成後,為了將控制點匯入Pix4Dmapper 中以校正整 體模型之真實座標,選定工學院樓頂作為參考站,以單點定位方式,自 2018 年 3 月 2 日起觀測至 2018 年 3 月 5 日,並為了去除架設初期精度較 差之結果,只擷取中間48 小時,共計 5760 筆定位結果,並將結果加以 平均,結果如表4.1 所示。決定參考站座標後,於 2018 年 3 月 9 日以 RTK 方式測量每個候選控制點的三維座標,每個點位觀測5 分鐘,並將結果 加以平均,如表4.2 所示,並以圖 4.8 表示控制點在拼接後影像上之分布 圖,其中黃標為參考站,紅標為控制點。

(49)

40 圖4.8 校區空拍全貌與候選控制點分布示意圖 表4.1 參考站之三維座標 工學院頂樓參考站 緯度(°N) 經度(°E) 高程(m) 3 月 2 日 ~3 月 5 日 平均 22.73317481 120.27635677 31.685 標準差(m) 0.519 0.522 1.451

(50)

41 表 4.2 候選控制點之三維座標 編號 緯度(°N) 經度(°E) 高程(m) 0 22.73339415 120.27708238 4.767 1 22.73424787 120.27751516 5.213 2 22.73515157 120.27745700 5.298 3 22.73600868 120.27995472 6.140 4 22.73644376 120.28108071 5.490 5 22.73567636 120.28377579 5.842 6 22.73526771 120.28521259 5.823 7 22.73595919 120.28653052 5.910 8 22.73713092 120.28620037 5.557 9 22.734596931 120.28934672 5.926 10 22.73332655 120.28930004 5.631 11 22.73276265 120.29022532 4.990 12 22.73331680 120.28698335 5.833 13 22.73444851 120.28599863 11.376 14 22.73302181 120.28457260 5.994 15 22.73472518 120.28332873 5.898 16 22.73197843 120.27986124 4.894 17 22.73142519 120.27816572 5.262 18 22.73182958 120.27757101 5.322 19 22.73206007 120.27739995 5.319

(51)

42

4.2 探討觀測時間對 GNSS 定位之影響

控制點之精度乃是影像拼接及三維數值模型校正之重要指標,本研 究之實驗場址為國立高雄大學校區,並採用RTK 技術量測控制點三維座 標,共量得20 個點。由於量測所有的點並非於同一時刻量得,故需驗證 時間對RTK 測量之影響。本節將分為兩個部分探討,一為同點同日不同 時段對定位精度之影響,一為同點同時段不同日對定位精度之影響。以 下將詳細介紹本實驗之方法與結果討論。 4.2.1 實驗方法 本研究選定了點0 及點 16 做為測試點位,分別於 2018 年 3 月 29 日、 3 月 30 日、4 月 2 日進行定位,於 RTK 產生固定解後,每個點持續定位 20 分鐘,並將定位結果加以平均,與 3 月 9 日全校測量時之結果比較。 每天共分為3 個觀測時段,如表 4.3 所示。定位結果將與 2018 年 3 月 9 日之全校測量結果做比較,以驗證當日量測結果之可靠性。本次實驗分 兩部分,首先比較同一日內,不同時段定位結果之差異;其次,比較同 一時段,不同日所定位結果之差異。表4.4、表 4.5 分別表示了點 0 與點 16 在同一日內,不同時段之定位結果。表 4.6、表 4.7、表 4.8 分別表示 了點0 同時段不同日之定位結果;表 4.9、表 4.10、表 4.11 分別表示了點 16 同時段不同日之定位結果。 表4.3 觀測時段表 時段 觀測時間 0 16 1 10:40-11:00 11:10-11:30 2 14:10-14:30 14:40-15:00 3 17:40-18:00 18:10-18:30

(52)

43 表4.4 點 0 同日不同時段之定位結果,實驗結果編號 A A 3 月 30 日 緯度(°N) 經度(°E) 高程(m) 時段 1 A1 平均 22.73339424 120.27708226 4.762 標準差(m) 0.002 0.002 0.006 時段 2 A2 平均 22.73339424 120.27708226 4.768 標準差(m) 0.002 0.002 0.007 時段 3 A3 平均 22.73339425 120.27708225 4.767 標準差(m) 0.002 0.002 0.005 表4.5 點 16 同日不同時段之定位結果,實驗結果編號 B B 3 月 30 日 緯度(°N) 經度(°E) 高程(m) 時段 1 B1 平均 22.73197853 120.27986115 4.914 標準差(m) 0.002 0.003 0.009 時段 2 B2 平均 22.73197851 120.27986117 4.910 標準差(m) 0.003 0.007 0.019 時段 3 B3 平均 22.73197856 120.27986114 4.919 標準差(m) 0.004 0.004 0.010

(53)

44 表4.6 點 0 時段 1 不同日之定位結果,實驗結果編號 C1 C1 時段1(10:40-11:00) 緯度(°N) 經度(°E) 高程(m) C11 3 月 30 日 平均 22.73339424 120.27708226 4.762 標準差(m) 0.002 0.002 0.006 C12 4 月 2 日 平均 22.73339423 120.27708225 4.769 標準差(m) 0.002 0.002 0.006 3 月 9 日 平均 22.73339415 120.27708238 4.767 標準差(m) 0.002 0.002 0.006 表4.7 點 0 時段 2 之不同日定位結果,實驗結果編號 C2 C2 時段2(14:10-14:30) 緯度(°N) 經度(°E) 高程(m) C21 3 月 30 日 平均 22.73339424 120.27708226 4.768 標準差(m) 0.002 0.002 0.007 C22 4 月 2 日 平均 22.73339423 120.27708226 4.768 標準差(m) 0.002 0.002 0.009 3 月 9 日 平均 22.73339415 120.27708238 4.767 標準差(m) 0.002 0.002 0.006

(54)

45 表4.8 點 0 時段 3 之不同日定位結果,實驗結果編號 C3 C3 時段 3(17:40-18:00) 緯度(°N) 經度(°E) 高程(m) C31 3 月 29 日 平均 22.73339432 120.27708224 4.767 標準差(m) 0.003 0.003 0.007 C32 3 月 30 日 平均 22.73339425 120.27708225 4.767 標準差(m) 0.002 0.002 0.005 3 月 9 日 平均 22.73339415 120.27708238 4.767 標準差(m) 0.002 0.002 0.006 表4.9 點 16 時段 1 之不同日定位結果,實驗結果編號 D1 D1 時段1(11:10-11:30) 緯度(°N) 經度(°E) 高程(m) D11 3 月 30 日 平均 22.73197853 120.27986115 4.914 標準差(m) 0.002 0.003 0.009 D12 4 月 2 日 平均 22.73197848 120.27986120 4.903 標準差(m) 0.003 0.003 0.013 3 月 9 日 平均 22.73197843 120.27986124 4.894 標準差(m) 0.002 0.003 0.010

(55)

46 表4.10 點 16 時段 2 之不同日定位結果,實驗結果編號 B2 D2 時段 2(14:40-14:60) 緯度(°N) 經度(°E) 高程(m) D21 3 月 30 日 平均 22.73197851 120.27986117 4.910 標準差(m) 0.003 0.007 0.019 D22 4 月 2 日 平均 22.73197848 120.27986119 4.913 標準差(m) 0.003 0.004 0.009 3 月 9 日 平均 22.73197843 120.27986124 4.894 標準差(m) 0.002 0.003 0.010 表4.11 點 16 時段 3 之不同日定位結果,實驗結果編號 D3 B3 時段 3(17:40-18:00) 緯度(°N) 經度(°E) 高程(m) D31 3 月 29 日 平均 22.73197846 120.27986167 4.911 標準差(m) 0.004 0.003 0.008 D32 3 月 30 日 平均 22.73197856 120.27986114 4.919 標準差(m) 0.004 0.004 0.010 3 月 9 日 平均 22.73197843 120.27986124 4.894 標準差(m) 0.002 0.003 0.010

(56)

47 4.2.2 結果與討論 本研究採用Vincenty 公式(Vincenty, 1975)計算兩點距離,將與 3 月 9 日之RTK 定位結果比較,其間的差異則定義為「較差」,並將計算結果 繪製如表4.12、圖 4.9、圖 4.10、表 4.13、圖 4.11、圖 4.12。 從表4.12 可以發現,在同一日內,即使觀測時段不同,彼此間與 3 月9 日之定位結果差異不大,水平較差皆在 2 公分以內,高程較差皆在 2.5 公分以內,顯見同日內不同時段的測量結果差異非常小,而從圖 4.10 可以明顯發現,點16 定位之高程較差普遍比點 0 來大,可能是因為點 16 附近有較多的樹,影響了衛星觀測的品質,造成定位的穩定性及可靠性 不如點0。從表 4.13 來看,僅有兩筆定位結果水平較差大於 2 公分,其 中D31 達到 0.044 公尺,很可能是人為架設儀器時產生較大的誤差。 表4.12 同日不同時段 RTK 定位結果與 3 月 9 日 RTK 定位結果之較差 同日不同時段 RTK 定位結果與 3 月 9 日 RTK 定位結果之較差 編號 水平較差(m) 高程較差(m) A1 0.016 0.005 A2 0.016 0.001 A3 0.017 0.000 B1 0.014 0.020 B2 0.011 0.016 B3 0.018 0.025

(57)

48

圖4.9 同日不同時段之水平較差比較圖

(58)

49 表 4.13 同時段不同日 RTK 定位結果之較差 同時段不同日 RTK 定位結果之較差 編號 水平較差(m) 高程較差(m) C11 0.016 0.005 C12 0.016 0.002 C21 0.016 0.001 C22 0.015 0.001 C31 0.024 0.000 C32 0.017 0.000 D11 0.014 0.020 D12 0.007 0.009 D21 0.011 0.016 D22 0.008 0.019 D31 0.044 0.017 D32 0.018 0.025

(59)

50

圖4.11 同時段不同日之水平較差比較圖

(60)

51 總結本節實驗之結果,在同一日內施作 RTK 定位時,不同時段的影 響不大,水平較差皆在2 公分以內,高程在 2.5 公分以內;而不同日同時 段之測量結果顯示,大部分實驗結果之水平較差在2 公分以內,高程較 差亦不超過2.5 公分,而本研究之校區空拍圖每個像素之實際大小約略為 12.8 公分,2 公分之較差僅 0.156 個像素,因此以單頻 GNSS 衛星接收儀 施作RTK,應用於影像接合所需控制點之三維座標測量有其可靠性。

4.3 應用 DIC 於控制點定位

Pix4Dmapper 拼接完成之影像與三維數值模型,其三維座標乃是 基於搭載於UAV 上,低精度 GPS 所提供之座標,透過空中三角測量與 內、外方位參數經最佳化程序所計算得出,往往與真實三維座標相差甚 大,為了校正接合影像與數值模型之座標,必須加入具備高精度三維座 標之控制點。在過去,往往透過人工選點的方式,定出同一個控制點在 不同影像上之影像座標,此種方法較容易產生人為造成的隨機誤差。回 顧過去相關研究,控制點對於接合影像與三維數值模型之真實三維座標 校正有顯著的效果,且應用數位影像相關係數法於具有高解析度之影像 上,的確能改善單純使用人工選點定出控制點影像座標之精度(周祐諒, 2016)。然而現地實驗中,UAV 飛行高度甚高,地面解析度不佳,因此應 用DIC 進行控制點比對之可行性猶未可知。本研究擬以人工選點方式定 出控制點在第一張影像上之影像座標後,再以DIC 比對出控制點在其他 影像上之座標,並與單純使用人工選點進行精度比較。

(61)

52 4.3.1 控制點與檢核點 本研究從20 個候選控制點中,選定其中 15 個點,其分佈如圖 4.13 所示,其中8 個做為控制點,以紅色星號表示,7 個做為檢核點,以黃色 星號表示。 圖4.13 控制點與檢核點分佈示意圖 4.3.2 實驗步驟與結果 本實驗同時探討了控制點標記影像數量對影像拼接精度影響之變化, 而控制點標記影像數量至少得為兩張以上,因此將從兩張影像逐步增加 到五張影像,並觀察檢核點精度之變化,而進行影像標記時,將以影像 中心之三維座標最接近該控制點三維座標之影像,作為初始影像。步驟 如下: 1. 在 Pix4Dmapper 中將每個控制點與檢核點,以人工選點方式標記在兩 張影像上,並將成果輸出,製成表 4.14。

(62)

53 2. 將第一步驟中,在初始影像上人工標記之控制點影像座標輸出,匯入 DIC 中,以初始影像比對出第二張影像上,控制點之影像座標並匯出, 再將比對出之影像座標與初始影像之座標,匯入 Pix4Dmapper 中進行 影像拼接,並將成果輸出,製成表 4.15。 3. 逐步增加影像標記數量,步驟如同第一步驟與第二步驟,製成表 4.16、 表 4.17、表 4.18、表 4.19、表 4.20、表 4.21 與表 4.22。

(63)

54

表 4.14 人工選點標記 2 張影像之檢核點精度

Check Point Error X [m] Error Y [m] Error Z [m] 3 0.0947 0.4000 4.4479 5 -0.1931 -0.2988 3.8927 6 -0.1229 -0.0737 3.5555 7 -0.1246 0.2054 3.6143 10 -0.0659 -0.6014 0.3695 13 -0.0421 -0.0394 1.6933 14 -0.0210 -0.1924 2.1489 Mean [m] -0.067842 -0.085739 2.817443 Sigma [m] 0.085547 0.303731 1.346800 RMS [m] 0.109183 0.315600 3.122796 表4.15 DIC 選點標記 2 張影像之檢核點精度

Check Point Error X [m] Error Y [m] Error Z [m] 3 0.1065 0.3401 3.9111 5 -0.1224 -0.5047 2.5855 6 -0.1181 -0.2876 2.2541 7 -0.0859 0.0714 2.0606 10 -0.0512 -0.5728 -0.1227 13 -0.0013 -0.0241 1.5053 14 0.0508 -0.2404 2.0213 Mean [m] -0.031658 -0.174030 2.030745 Sigma [m] 0.080930 0.300646 1.121730 RMS [m] 0.086901 0.347382 2.319958

(64)

55

表 4.16 人工選點標記 3 張影像之檢核點精度

Check Point Error X [m] Error Y [m] Error Z [m] 3 0.0500 0.4093 1.5831 5 -0.1805 -0.4466 1.0301 6 -0.0588 -0.2025 1.2033 7 -0.1182 0.1062 1.2679 10 -0.0552 -0.5958 0.4776 13 -0.0102 -0.0445 0.7859 14 0.0028 -0.1895 1.3025 Mean [m] -0.052873 -0.137631 1.092912 Sigma [m] 0.071825 0.311595 0.338870 RMS [m] 0.089188 0.340637 1.144241 表4.17 DIC 選點標記 3 張影像之檢核點精度

Check Point Error X [m] Error Y [m] Error Z [m] 3 0.0704 0.3540 1.0151 5 -0.0647 -0.4854 0.4924 6 -0.0648 -0.2755 0.8578 7 -0.0207 0.0153 0.5676 10 -0.0522 -0.5572 0.3268 13 0.0321 -0.0504 0.6839 14 0.0768 -0.2440 1.2231 Mean [m] -0.003307 -0.177598 0.738105 Sigma [m] 0.057739 0.290029 0.289791 RMS [m] 0.057833 0.340085 0.792955

(65)

56

表 4.18 人工選點標記 4 張影像之檢核點精度

Check Point Error X [m] Error Y [m] Error Z [m] 3 0.0546 0.4222 1.2087 5 -0.1812 -0.4540 0.5024 6 -0.0603 -0.1949 0.7502 7 -0.1192 0.0593 0.5332 10 -0.0892 -0.5897 0.2354 13 -0.0113 -0.0529 0.6779 14 0.0032 -0.1809 0.8789 Mean [m] -0.057632 -0.141545 0.683823 Sigma [m] 0.074174 0.309111 0.286528 RMS [m] 0.093932 0.3399978 0.741426 表4.19 DIC 選點標記 4 張影像之檢核點精度

Check Point Error X [m] Error Y [m] Error Z [m] 3 0.0872 0.3661 0.8184 5 -0.0901 -0.4612 0.3202 6 -0.0816 -0.2854 0.4692 7 -0.0533 -0.0055 0.2365 10 -0.0686 -0.5429 0.2929 13 0.0223 -0.0421 0.6575 14 0.0591 -0.2195 0.9843 Mean [m] -0.017850 -0.170071 0.539849 Sigma [m] 0.067270 0.285466 0.265965 RMS [m] 0.069598 0.3322287 0.601809

(66)

57

表 4.20 人工選點標記 5 張影像之檢核點精度

Check Point Error X [m] Error Y [m] Error Z [m] 3 0.1001 0.4617 0.1963 5 -0.1072 -0.3733 -0.0884 6 -0.0255 -0.2314 -0.0145 7 -0.0488 -0.0041 0.0259 10 -0.1085 -0.5900 0.2023 13 0.0323 -0.0813 0.4470 14 0.0450 -0.1747 0.5946 Mean [m] -0.016092 -0.141863 0.194740 Sigma [m] 0.073238 0.304822 0.231781 RMS [m] 0.074985 0.336216 0.302731 表4.21 DIC 選點標記 5 張影像之檢核點精度

Check Point Error X [m] Error Y [m] Error Z [m] 3 0.0618 0.3748 0.2495 5 -0.0718 -0.4274 -0.2079 6 -0.0595 -0.2700 0.0502 7 -0.0227 -0.0111 -0.0679 10 -0.0448 -0.5314 0.2091 13 0.0399 -0.0665 0.3404 14 0.0774 -0.2255 0.5666 Mean [m] -0.002820 -0.165312 0.162841 Sigma [m] 0.056748 0.278415 0.241304 RMS [m] 0.056819 0.323795 0.291110

(67)

58 4.3.3 結果與討論 從表4.14、表 4.15 可以發現,不論單純人工選點,或是 DIC 選點, 若是僅將控制點標記於兩張影像上,雖然相較於人工選點,DIC 選點於 高程誤差標準差精度改善約為17%,但其高程誤差標準差卻達 1.12 公尺, 並不適合應用於三維數值模型之建立。接著從表4.16、表 4.17 可以發現, 當控制點標記於3 張影像時,相較於僅標記 2 張影像,DIC 選點在 x 方 向之誤差標準差精度改善24.40%,y方向之誤差標準差精度僅改善約 7%, 高程誤差標準差之精度改善14.48%;而人工選點則在 y 方向上有變差的 情況發生,這很可能是因為標記於第3 張影像上時,產生較大的人為誤 差。表4.18、表 4.19 的結果顯示,x 方向之誤差標準差改善約為 10.26%, y 方向之誤差標準差改善約 7.65%,高程之誤差標準差改善 23.76%,其 中x 方向誤差平均值較標記於 3 張影像時大,可能是因為第 4 張影像與 初始影像差異較大,造成DIC 比對的結果於某方向累積較大的誤差;表 4.20、表 4.21 顯示,x 方向上的誤差平均值又降低了,甚至低於標記 3 張影像時之結果,可能是DIC 比對第 5 張影像與初始影像上產生的誤差, 與第4 張影像之誤差抵銷的結果所致。為了方便比較人工選點與 DIC 選 點精度差異改善趨勢,將前一節的表格重新繪製成折線圖,標如圖4.14、 圖4.15、圖 4.16、圖 4.17、圖 4.18 及圖 4.19 所示。

(68)

59 圖4.14 x 方向檢核點之誤差標準差 在x 方向上,不論標記於幾張影像上,相較於單純採用人工選點, DIC 選點皆能提升控制點定位之精度,惟當欲標記之影像與初始影像相 差較大,如影像有旋轉或是真實座標距離較遠時,會增加DIC 比對之難 度,因此選定控制點之初始影像後,建議選擇初始影像拍攝序號前後, 飛行方向相同之影像,選定包含初始影像在內奇數張影像進行標記,可 能可以降低在某方向上,DIC 所累積的誤差。

(69)

60 圖4.15 y 方向檢核點之誤差標準差 圖4.15 顯示,在 y 方向上 DIC 選點亦能改善現地控制點定位之精 度,但相比於x 方向上之誤差標準差,y 方向之誤差顯然大得多,可能是 因為UAV 飛行方向為 x 方向,而標記之控制點影像座標,其在 x 方向上 的位移較飛行方向一致,因而在x 方向上有著較大幅度的改善,而 y 方 向上的改善則相當有限,若能提高影像之解析度,使得在影像有較大變 形時仍能以DIC 比對出控制點影像座標,則有望改善 y 方向上之誤差。

數據

圖 4.1 GNSS 接收儀
圖 4.9 同日不同時段之水平較差比較圖
圖 4.11 同時段不同日之水平較差比較圖
表 4.14 人工選點標記 2 張影像之檢核點精度
+4

參考文獻

相關文件

運用 Zuvio IRS 與台日比較文化觀點於日本文化相關課程之教學研究 Applying Zuvio IRS and Perspective on Cultural comparison between Taiwan and Japan to Teaching

由於醫療業導入 ISO 9000 品保系統的「資歷」相當資淺,僅有 三年多的年資 11 ,因此,對於 ISO 9000 品保系統應用於醫療業之相關 研究實在少之又少,本研究嘗試以通過

Mutual information is a good method widely used in image registration, so that we use the mutual information to register images.. Single-threaded program would cost

譚志忠 (1999)利用 DEA 模式研究投資組合效率指數-應用

Jones, "Rapid Object Detection Using a Boosted Cascade of Simple Features," IEEE Computer Society Conference on Computer Vision and Pattern Recognition,

本研究是以景觀指數進行對 1993 年、2008 年與擴大土地使用三個時期之評 估,其評估結果做比較討論。而目前研究提供研究方法的應用-GIS 與 FRAGSTATS 之使用方法。從 1993 年至

“ Ghosts and Ancestors in Medieval Chinese Religion: The Yu-lan-p'en Festival as Mortuary Ritual.” History of Religions 26, No. Tillman,

第三節 負數概念之 負數概念之 負數概念之迷思概念 負數概念之 迷思概念 迷思概念相關研究 迷思概念 相關研究 相關研究 相關研究