運用曲面擬合提升幾何法大地起伏值精度之研究 - 政大學術集成

全文

(1)國立政治大學地政學系. 碩士論文. 私立中國地政研究所. 運用曲面擬合提升幾何法大地起伏值政治大立精度之研究. ‧. ‧ 國. 學. The Study of Applying Surface Fitting to Improve Geometric Geoidal Undulation n. er. io. sit. y. Nat. al. i n U. Ch. v. engchi 研究生：蔡名曜指導教授：甯方璽. 中. 華. 民. 國. 一. 零. 三. 年. 六. 月.

(2) 立. 政治大. ‧. ‧ 國. 學. n. er. io. sit. y. Nat. al. Ch. engchi. i n U. v.

(3) 摘要大地起伏值為正高與橢球高的差異量，如果取得高精度的大地起伏值，可以利用衛星定位測量施測橢球高並計算得到高精度的正高，其成本低廉，可望取代傳統的水準測量。而大地起伏值可以分為幾何法或重力法的大地起伏值，其中幾何法的大地起伏值計算方法簡易且精度高，可以利用曲面擬合方法取得之。但是幾何法的大地起伏值會受到地形起伏的影響，大範圍的曲面擬合會降低其精度。台灣的地形起伏大，難以進行大範圍曲面擬合。於是本研究利用環域方法搜尋待測點位鄰近的水準點參與曲面方程式擬合大地起伏，試圖找到最適合的大地起伏擬合範圍。成果顯示：環域的範圍從 10 公里至 30 公里，利用二次曲面方程式擬合大地起伏在台灣平地區域能夠達到預測精度與內部精度同時低於 5 公分。另外由於衛星定位測量橢球高的誤差較高，需進行資料品質評估並進行粗差偵測。針對粗差偵測提出新的方法，利用最佳化演算法中的量子行為粒子群演算法計算最小二乘平差法中的權矩陣，期望能夠將粗差觀測量的權重降低，達到粗差偵測的效果。成果顯示最佳化權矩陣演算法，能夠將粗差對平差系統的影響量降到最低。. 立. 政治大. ‧. ‧ 國. 學. io. sit. y. Nat. 本研究建立一套台灣地區的大地起伏擬合作業程序：利用環域搜尋鄰近水準點、曲面方程式及環域範圍選擇與資料的粗差偵測，可獲得高品質的大地起伏。. n. al. er. 關鍵字：大地起伏、曲面擬合、量子行為粒子群算法、粗差偵測. Ch. engchi. I. i n U. v.

(4) Abstract The geoidal undulation is the difference of ellipsoid height and orthometric height. We can obtain high accuracy of orthometric height by existing high accuracy of geoidal undulation and the ellipsoidal height measuring by GPS. It expected to replace the traditional leveling survey due to the less cost. This study uses buffer method to search the leveling benchmarks around the object point, attempts to find the proper range of fitting geoidal undulation to curve surface. Experimental results shows that it can archive 5cm level on both prediction error and internal precision by fitting geoidal undulation on 2nd curve surface model where the buffer range is from 10 km to 30 km. In this study, also uses the quantum-behaved particle swarm optimization to calculate the weight matrix of least square adjustment, the purpose is to down-weighting the suspicious outlier, and detect the outlier. Experimental results shows that the optimal weight matrix algorithm can reduce the influence of outlier.. 立. 政治大. ‧ 國. 學. This study establish a procedure of fitting geoidal undulation: using buffer analysis to search the adjacent leveling benchmark, selecting the proper buffer range and surface equation and detecting outlier in data.. ‧. Keyword: Geoid Height, Curve Fitting, Quantum-behaved Particle Swarm Optimization, Outlier Detection. n. er. io. sit. y. Nat. al. Ch. engchi. II. i n U. v.

(5) 目錄第一章緒論 .......................................................................................1 第一節研究背景 ................................................................................................... 1 第二節研究動機及目的 ....................................................................................... 4 第三節論文架構 ................................................................................................... 5. 第二章文獻回顧...............................................................................7 第一節大地起伏 ................................................................................................... 7. 政治大. 第二節粗差偵測理論 ........................................................................................... 8. 立. 第三節粒子群演算法 ........................................................................................... 9. ‧ 國. 學. 第三章理論基礎.............................................................................11. ‧. 第一節大地起伏概念 ......................................................................................... 11. sit. y. Nat. 第二節大地起伏值之解算方法 ......................................................................... 12. io. er. 一、以全球大地位模式解算大地起伏 ......................................................... 13 二、 Stokes 積分公式法 ................................................................................. 14. n. al. Ch. i n U. v. 三、衛星測高法 ............................................................................................. 14. engchi. 四、幾何擬合法 ............................................................................................. 14 第三節粗差偵測原理 ......................................................................................... 15 一、可靠度理論 ............................................................................................. 15 二、資料探測法 ............................................................................................. 18 三、選擇權迭代法 ......................................................................................... 20 第四節粒子群演算法原理 ................................................................................. 21 一、基本粒子群演算法理論 ......................................................................... 21 二、基本粒子群演算法算法 ......................................................................... 22 三、量子行為粒子群演算法理論 ................................................................. 24. III.

(6) 四、量子行為粒子群演算法算法 ................................................................. 24 第五節最佳化權矩陣 ......................................................................................... 26 一、基本構想 ................................................................................................. 26 二、演算方法 ................................................................................................. 27. 第四章研究方法 ............................................................................ 30 第一節研究流程 ................................................................................................. 30 第二節研究資料 ................................................................................................. 31 一、內政部公告之一等水準點 ..................................................................... 31. 政治大第三節曲面方程式分析方法 ............................................................................. 34 立二、台中市水準點資料 ................................................................................. 33. 一、過度擬合 ................................................................................................. 34. ‧ 國. 學. 二、交叉驗證法 ............................................................................................. 38. ‧. 三、研究方法 ................................................................................................. 39 第四節曲面擬合範圍分析方法 ......................................................................... 40. y. Nat. io. sit. 一、擬合的環域範圍選擇方法 ..................................................................... 40. n. al. er. 二、環域中心點選擇方法 ............................................................................. 41. i n U. v. 三、曲面擬合的限制範圍 ............................................................................. 43. Ch. engchi. 四、研究方法 ................................................................................................. 44 第五節粗差偵測分析方法 ................................................................................. 45 一、研究對象 ................................................................................................. 45 二、研究方法 ................................................................................................. 49 第六節台灣地區大地起伏擬合作業程序 ......................................................... 53 一、擬合作業程序 ......................................................................................... 53 二、研究方法 ................................................................................................. 53. 第五章實驗成果 ............................................................................ 55 第一節曲面方程式初步實驗成果與討論 ......................................................... 55 IV.

(7) 第二節曲面擬合範圍實驗成果與討論 ............................................................. 56 一、二次曲面的擬合成果 ............................................................................. 56 二、三次曲面的擬合成果 ............................................................................. 69 三、討論 ......................................................................................................... 81 第三節粗差偵測分析實驗成果 ......................................................................... 85 一、模擬資料實驗成果 ................................................................................. 85 二、台中市資料實驗成果 ............................................................................. 98 第四節臺灣地區大地起伏擬合作業程序實驗成果 ....................................... 106 一、僅搜尋內政部公告之一等水準點 ....................................................... 107. 政治大. 二、同時搜尋內政部及台中市水準點 ....................................................... 108. 立. 三、討論 ....................................................................................................... 109. ‧ 國. 學. 第六章結論與建議.......................................................................111. ‧. 第一節結論 ....................................................................................................... 111 第二節建議 ....................................................................................................... 111. y. Nat. al. v i n Ch 外文參考文獻 ....................................................................................... 114 engchi U. 中文參考文獻 ....................................................................................... 113. n. 二、. er. io. 一、. sit. 參考文獻 ........................................................................................113. V.

(8) 圖目錄圖 1 正高、橢球高及大地起伏的關係圖.............................................. 12 圖 2 正高、橢球高及大地起伏的實際關係.......................................... 12 圖 3 機率密度函數因為觀測量中的粗差而產生位移(Lehmann, 2013) ............................................................................................................ 18 圖 4 粒子群演算法粒子搜尋極端值的過程........................................... 22 圖 5 基本粒子群演算法演算的流程....................................................... 23 圖 6 量子行為粒子群演算法演算的流程............................................... 26. 政治大. 圖 7 研究方法流程圖.............................................................................. 30. 立. 圖 8 本研究之一等一級及一等二級水準點分布圖............................... 31. ‧ 國. 學. 圖 9 本研究之一等水準點橢球高標準偏差圖（單位：cm） .............. 33 圖 10 台中市政府公共管線資料庫系統建置案所測設之水準點........ 34. ‧. 圖 11 一次曲面擬合模擬點位的成果，內部精度為 0.2723 單位 ....... 36. y. Nat. io. sit. 圖 12 三次曲面擬合模擬點位的成果，內部精度為 0.2703 單位........ 36. n. al. er. 圖 13 五次曲面擬合模擬點位的成果，內部精度為 0.2637 單位........ 37. Ch. i n U. v. 圖 14 在不同模型複雜度下的訓練樣本誤差及測試樣本誤差（Hastie T., et al., 2009） ................................................................................. 37. engchi. 圖 15 交叉驗證法流程範例.................................................................... 38 圖 16 曲面方程式分析方法流程圖........................................................ 40 圖 17 以環域方法選取待擬合點位之示意圖........................................ 41 圖 18 以環域 5 公里擬合台灣大地起伏所涵蓋的區域........................ 42 圖 19 以環域 10 公里擬合台灣大地起伏所涵蓋的區域....................... 43 圖 20 以環域 15 公里擬合台灣大地起伏所涵蓋的區域....................... 43 圖 21 台灣無水準點位之山區示意圖.................................................... 44 圖 22 曲面擬合範圍分析流程圖............................................................ 45 VI.

(9) 圖 23 規則網格分布模擬點位資料（4X4）單位：m .......................... 46 圖 24 規則網格分布模擬點位資料（6X6）單位：m .......................... 46 圖 25 本研究台中市水準點位編號圖 .................................................... 48 圖 26 以模擬點位偵測粗差的實驗分類圖 ............................................ 49 圖 27 16 點規則網格分布模擬點位資料含有兩個粗差之真誤差圖 .... 50 圖 28 16 點規則網格分布模擬點位資料含有三個粗差之真誤差圖 .... 50 圖 29 36 點規則網格分布模擬點位資料含有兩個粗差之真誤差圖 .... 51 圖 30 36 點規則網格分布模擬點位資料含有三個粗差之真誤差圖 .... 51. 政治大圖 32 大地起伏擬合作業程序實驗流程 ................................................ 54 立圖 31 36 點規則網格分布模擬點位資料含有四個粗差之真誤差圖 .... 52. ‧ 國. 學. 圖 33 預測誤差與內部精度比較圖（台中實驗區） ............................ 55. ‧. 圖 34 二次曲面方程式擬合大地起伏 5 公里環域成果（內部精度，單位：cm） ........................................................................................... 57. sit. y. Nat. 圖 35 二次曲面方程式擬合大地起伏 5 公里環域成果（預測誤差，單位：cm） ........................................................................................... 58. er. io. 圖 36 二次曲面方程式擬合大地起伏 10 公里環域成果（內部精度，單位：cm） ........................................................................................... 59. al. n. v i n Ch 圖 37 二次曲面方程式擬合大地起伏公里環域成果（預測誤差，單 i U e n g c h 10 位：cm） ........................................................................................... 60. 圖 38 二次曲面方程式擬合大地起伏 15 公里環域成果（內部精度，單位：cm） ........................................................................................... 61 圖 39 二次曲面方程式擬合大地起伏 15 公里環域成果（預測誤差，單位：cm） ........................................................................................... 62 圖 40 二次曲面方程式擬合大地起伏 20 公里環域成果（內部精度，單位：cm） ........................................................................................... 63 圖 41 二次曲面方程式擬合大地起伏 20 公里環域成果（預測誤差，單位：cm） ........................................................................................... 64 圖 42 二次曲面方程式擬合大地起伏 25 公里環域成果（內部精度，單 VII.

(10) 位：cm） ........................................................................................... 65 圖 43 二次曲面方程式擬合大地起伏 25 公里環域成果（預測誤差，單位：cm） ........................................................................................... 66 圖 44 二次曲面方程式擬合大地起伏 30 公里環域成果（內部精度，單位：cm） ........................................................................................... 67 圖 45 二次曲面方程式擬合大地起伏 30 公里環域成果（預測誤差，單位：cm） ........................................................................................... 68 圖 46 三次曲面方程式擬合大地起伏 5 公里環域成果（內部精度，單位：cm） ........................................................................................... 69 圖 47 三次曲面方程式擬合大地起伏 5 公里環域成果（預測誤差，單位：cm） ........................................................................................... 70. 政治大圖 48 三次曲面方程式擬合大地起伏 10 公里環域成果（內部精度，單立位：cm） ........................................................................................... 71. ‧ 國. 學. 圖 49 三次曲面方程式擬合大地起伏 10 公里環域成果（預測誤差，單位：cm） ........................................................................................... 72. ‧. sit. y. Nat. 圖 50 三次曲面方程式擬合大地起伏 15 公里環域成果（內部精度，單位：cm） ........................................................................................... 73. er. io. 圖 51 三次曲面方程式擬合大地起伏 15 公里環域成果（預測誤差，單位：cm） ........................................................................................... 74. al. n. v i n Ch 圖 52 三次曲面方程式擬合大地起伏 20 公里環域成果（內部精度，單 U i e h n c g 位：cm） ........................................................................................... 75 圖 53 三次曲面方程式擬合大地起伏 20 公里環域成果（預測誤差，單位：cm） ........................................................................................... 76 圖 54 三次曲面方程式擬合大地起伏 25 公里環域成果（內部精度，單位：cm） ........................................................................................... 77 圖 55 三次曲面方程式擬合大地起伏 25 公里環域成果（預測誤差，單位：cm） ........................................................................................... 78 圖 56 三次曲面方程式擬合大地起伏 30 公里環域成果（內部精度，單位：cm） ........................................................................................... 79 圖 57 三次曲面方程式擬合大地起伏 30 公里環域成果（預測誤差，單. VIII.

(11) 位：cm） ........................................................................................... 80 圖 58 環域大小與曲面方程式關係圖(內部精度) ................................. 82 圖 59 環域大小與曲面方程式關係圖(預測誤差) .................................. 82 圖 60 二次曲面方程式擬合大地起伏最適的環域距離展示圖 ............. 83 圖 61 二次曲面方程式擬合大地起伏無法合理擬合之點位展示圖 ..... 85 圖 62 臺灣地區大地起伏擬合作業程序實驗：待測點位鄰近的水準點位圖（僅含內政部公告之一等水準點） ...................................... 107 圖 63 臺灣地區大地起伏擬合作業程序實驗：待測點位鄰近的水準點位圖 .................................................................................................. 108. 立. 政治大. ‧. ‧ 國. 學. n. er. io. sit. y. Nat. al. Ch. engchi. IX. i n U. v.

(12) 表目錄表 1 一等水準點統計資料（You, 2006） .............................................. 32 表 2 曲面擬合模擬資料點位表 (無單位) .............................................. 36 表 3 模擬點位資料的二次曲面方程式參數.......................................... 46 表 4 台中市水準點位資料品質數據（殘差單位：m） ....................... 47 表 5 台中市水準點位資料加入一般粗差量統計表（單位：m） ...... 53 表 6 台中市水準點位資料加入較大粗差量統計表（單位：m） ....... 53. 政治大表 8 二次曲面方程式擬合大地起伏 5 公里環域成果表....................... 58 立表 7 不同環域距離能夠找到的水準點位數量統計表.......................... 56. 表 9 二次曲面方程式擬合大地起伏 10 公里環域成果表..................... 60. ‧ 國. 學. 表 10 二次曲面方程式擬合大地起伏 15 公里環域成果表................... 62. ‧. 表 11 二次曲面方程式擬合大地起伏 20 公里環域成果表 ................... 64. sit. y. Nat. 表 12 二次曲面方程式擬合大地起伏 25 公里環域成果表................... 66. io. er. 表 13 二次曲面方程式擬合大地起伏 30 公里環域成果表................... 68. al. 表 14 三次曲面方程式擬合大地起伏 5 公里環域成果表.................... 70. n. v i n Ch 表 15 三次曲面方程式擬合大地起伏 ................... 72 e n g10c 公里環域成果表 hi U 表 16 三次曲面方程式擬合大地起伏 15 公里環域成果表................... 74 表 17 三次曲面方程式擬合大地起伏 20 公里環域成果表................... 76 表 18 三次曲面方程式擬合大地起伏 25 公里環域成果表................... 78 表 19 三次曲面方程式擬合大地起伏 30 公里環域成果表................... 80 表 20 規則網格分布 16 點的模擬點位資料含有兩個粗差成果表（粗差 =3σ） ............................................................................................... 87 表 21 規則網格分布 16 點的模擬點位資料含有兩個粗差成果表（粗差 =4σ） ............................................................................................... 87 表 22 規則網格分布 16 點的模擬點位資料含有三個粗差成果表（粗差 X.

(13) =3σ） ................................................................................................ 89 表 23 規則網格分布 16 點的模擬點位資料含有三個粗差成果表（粗差 =4σ） ................................................................................................ 89 表 24 規則網格分布 36 點的模擬點位資料含有 2 個粗差成果表（粗差 =3σ） ................................................................................................ 91 表 25 規則網格分布 36 點的模擬點位資料含有 2 個粗差成果表（粗差 =4σ） ................................................................................................ 92 表 26 規則網格分布 36 點的模擬點位資料含有 3 個粗差成果表（粗差 =3σ） ................................................................................................ 94 表 27 規則網格分布 36 點的模擬點位資料含有 3 個粗差成果表（粗差 =4σ） ................................................................................................ 95. 政治大表 28 規則網格分布立 36 點的模擬點位資料含有 4 個粗差成果表（粗差 =3σ） ................................................................................................ 96. ‧ 國. 學. 表 29 規則網格分布 36 點的模擬點位資料含有 4 個粗差成果表（粗差 =4σ） ................................................................................................ 97. ‧. sit. y. Nat. 表 30 規則網格分布 36 點模擬資料粗差個數在等權及最佳化權矩陣下的後驗中誤差表（單位=cm，先驗中誤差=3cm，粗差=3σ） 98. er. io. 表 31 台中市水準點位進行粗差偵測實驗一（*號代表含有 1m 大粗差點位，單位：m） ............................................................................. 99. al. n. v i n Ch 表 32 台中市水準點位進行粗差偵測實驗二（*號代表含 0.5m 大粗差 U i e h n c g 點位，單位：m） ........................................................................... 101 表 33 台中市水準點位進行粗差偵測實驗二（*號代表含 0.25m 大粗差點位，單位：m） ....................................................................... 104 表 34 臺灣地區大地起伏擬合作業程序實驗成果表（僅含內政部公告之一等水準點） .............................................................................. 108 表 35 臺灣地區大地起伏擬合作業程序實驗成果表 ........................... 109 表 36 台中市與內政部共同點位的差異表(單位：m) ........................ 109. XI.


(15) 第一章. 緒論. 第一節研究背景全球導航衛星系統（Global Navigation Satellite System, GNSS）測量為近年來快速發展且發展成熟的一項測量方法，GNSS 在大地測量、攝影測量及海洋測量等空間資訊領域中受到廣泛應用。以大地測量為例，以往的大地控制網使用三角測量或是天文測量等傳統方法來建立，但目前已經被 GNSS 所取代，因為 GNSS 測量可以在短時間之內以人力成本最低的狀況下取得高精度的平面坐標。儘管 GNSS 測量可以取得高精度的平面坐標，但是在高程部分受限於幾何因素，精度無法與平面坐標達到相同的等級（Meyer et al, 2006）；另外，因為高程系統本身就相對複雜，GNSS 所用的高程系統是橢球高（ ellipsoid height ）系統（ WGS84 ）並非水準測量所使用的正高（orthometric height）系統。於是 GNSS 所測得的高程必須再經過轉換後才能變成所需的正高。因為上述這些因素，GNSS 在高程測量上並不如平面測量方便及精準。. 立. 政治大. ‧. ‧ 國. 學. sit. y. Nat. 雖然 GNSS 在測量高程的精度等級劣於平面測量的精度，加上測得的高程資料並不能直接使用，但 GNSS 仍然具有低人力及時間成本的優點，使得許多研究開始發展利用 GNSS 進行高精度水準測量的方法。一般而. n. al. er. io. 言，傳統的水準測量需要大量時間成本及人力成本，更不用談及成本更高的精密水準測量，而在地形險惡地區進行水準測量更是艱鉅的挑戰。於是如果能夠使用 GNSS 測量高精度的水準資料，基於測量人員的安全及成本考量，可以避免許多過於崎嶇及危險的水準路線。於是發展出了 GPS 水準測量技術（GPS/leveling），利用 GPS 所測量得到的橢球高來推算正高。. Ch. engchi. i n U. v. GPS 水準測量技術是利用 GNSS 測量得到某一點的橢球高，加上已知的大地起伏值（geoid height 或是 geoidal undulation），利用幾何方法計算某一點的正高。由此可知，如果降低 GNSS 測量的誤差並搭配高精度的大地起伏值，可以利用 GNSS 得到高精度的正高，可望取代傳統的水準測量。臺灣高山地區佔據全島土地面積的百分之七十，因此要建立全臺灣的水準網並不容易。過去幾十年來，內政部國土測繪中心（前土地測量局）及陸軍聯勤測量隊利用傳統水準測量已建立全臺水準網，提供高精度的正高資料。儘管水準網已經建立，如果需要獲取某一點的正高，仍然要從已知的水準點引測計算高程差。如果精度的要求並不高，利用 GPS 水準測量技術可以快速獲得正高資料，省去許多時間及人力。 1.

(16) GPS 水準測量的前提是需要一個已知的大地起伏模型，而該模型可以利用物理或是幾何方法建立。以幾何概念來建立大地起伏並不考慮物理因素，只透過一個數學模式擬合已知的大地起伏值，並建立模型。利用幾何方式建立大地起伏模型的方法又被稱為擬合法，所計算出的大地起伏被稱為幾何法的大地起伏。通常使用二次曲面方程式作為擬合法的數學模式，並利用最小二乘法計算曲面方程式的係數。在進行曲面方程式擬合時，必須先進行粗差偵測將含有大誤差的觀測量剔除，以確保平差估計的參數不受到粗差影響。但是最小二乘平差的配賦誤差能力相當良好，使得少數的含有粗差觀測量會影響到所有觀測量。在討論粗差時，一般是使用觀測量經過平差後的殘差來判斷。但分析平差後的殘差時，粗差觀測量的殘差在平差過程中往往被稀釋，粗差反應在殘差中僅僅只是很小一部份。於是透過殘差分析來偵測粗差並不是相當理想的。. 政治大在粗差偵測理論中， Baarda （ Baarda, 1968 ）的數據探測法（ data 立 snooping）是透過計算觀測量的標準化殘差（normalized residuals）來判斷 ‧. ‧ 國. 學. 粗差。如果觀測量中不含粗差觀測，觀測量的標準化殘差應該符合使用者所定義的統計分布（常態分布、t 分布）；使用者決定顯著性水平（如 0.1% 或 0.5%）可以得到對應的臨界值（critical value），超過臨界值的觀測量則判斷為粗差。但數據探測法理論上對於單一粗差的偵測效果較佳。. Nat. sit. y. 如果平差系統存在多個粗差，粗差偵錯會變得相當困難，必須透過特. n. al. er. io. 殊的檢驗方法來定位粗差。例如進行多次的數據探測法，也是重複進行數據探測法，每次剔除一個可能的粗差，直到觀測量中不存在粗差。但多個粗差的狀況相當複雜，在不知道粗差數量及分布狀況下，使用數據探測法可能會誤判粗差觀測量。所以必須找到另一種可以提升偵測能力並降低誤判機率的粗差偵測方法。. Ch. engchi. i n U. v. 另一種粗差偵測的可能發展，為透過最佳化演算法來判斷粗差觀測量。目前基於計算機處理能力及運算效能的突破，許多過去需要耗時大量時間的數學問題，如今都能以合理的計算時間解算完成。許多最佳化演算法因此被發展出來解決最佳化問題（optimization problem）。最佳化問題簡單的說就是在某個問題之所有的決策方案中，找到其中一個最好的解決方案。目前有多種常用的最佳化演算法，例如：模擬退火法（ Simulated Annealing）、基因演算法（Genetic Algorithm）及粒子群演算法（Particle Swarm Optimization）。其中粒子群演算法是最佳化演算法中，利用集體智慧（Collective Intelligence）演算法概念來設計的一種方法。集體智慧演算法即透過許多個體互相競爭及合作，共同分享資訊來尋求一個最好的共同利益。粒子群演算法的優點在於計算快速、設定的參數較少、計算過程透. 2.

(17) 明及程式邏輯易懂等；但粒子群演算法的缺點在於只能保證收斂到區域最佳解。而量子行為粒子群演算法（ Quantum-behaved Particle Swarm Optimization）是粒子群演算法的一個改良算法，其需要控制的參數更少，且經過數學方法證明，該演算法可以找到全域的最佳解。基於最佳化演算法這個強力的計算工具，可以嘗試發展出有關粗差偵測的新方法。. 立. 政治大. ‧. ‧ 國. 學. n. er. io. sit. y. Nat. al. Ch. engchi. 3. i n U. v.

(18) 第二節研究動機及目的傳統上，正高資料需要透過實施高人力及時間成本的水準測量才能獲取。隨著全球定位系統的發展，發展出 GPS 水準測量技術，可以利用 GPS 直接測量獲得正高資料。但在實務上，GPS 水準測量技術需要取得待測點的大地起伏值，加上實測的橢球高計算待測點的正高值。因此理論上如果建立了一個高精度的大地起伏模型，可以取得高精度的正高資料。因此如何建立一個高精度的大地起伏模型是本文的研究重點。有鑑於臺灣的一等水準點數量充足，而且大部分點位含有橢球高資料，有利於使用擬合法建立大地起伏模型。相較於以重力方式建立的大地起伏模型，以擬合法建立的大地起伏模型的優勢在於其精度相當高，但劣勢在於需要充足的點位資料來解算擬合面。臺灣平地地區含有橢球高資料的一等水準點分布相當均勻，理論上可以建立相當高精度的大地起伏模型；. 立. 政治大. ‧. ‧ 國. 學. 本研究使用最小二乘法解算曲面方程式，接著利用環域分析（buffer analysis）來找到最適合的擬合範圍。相關研究指出，以幾何法擬合大地起伏模式，擬合精度與範圍相關：擬合範圍越大，精度越低。台灣的地形崎嶇複雜，合理的擬合範圍可能會因地而異。因此如何在台灣任何觀測點上找到最適合的擬合範圍來擬合大地起伏值，是需要透過一個完善的擬合策略來取得。. sit. y. Nat. n. al. er. io. 於是本研究的研究目標在於建立一個完善的擬合策略，包含如何選擇擬合範圍、擬合水準點位並分析擬合點位的品質。在理論上，擴大擬合範圍會使得擬合精度相對降低；但過小的擬合範圍內可能沒有足夠的已知水準點位。於是擬合範圍的大小必須透過適當取捨，而取捨的準則應該透過一定的評估基準。在決定擬合範圍大小的同時，應該同時分析擬合水準點位的資料品質，透過粗差偵測能夠以統計方法分析水準點位是否合乎精度標準。在粗差偵測方法上，本研究希望從最佳化演算法技術來調整觀測量權重的方向出發，偵測並弱化粗差對最小二乘法參數估值的的影響。. Ch. engchi. 4. i n U. v.

(19) 第三節論文架構本論文架構包含五個章節，各章節內容如下：第一章. 前言：. 說明本研究的研究背景、研究動機並回顧國內外有關的研究文章及論文。第二章. 文獻回顧：. 回顧過去有關大地起伏、粗差偵測及粒子群演算法的相關國內外研究。第三章. 理論基礎：. 政治大. 闡述大地起伏的理論基礎、粗差偵測及量子行為粒子群演算法的運算原理。. 立. ‧ 國. 研究方法：. 學. 第四章. ‧. 介紹研究區域及資料，說明本研究曲面方程式、曲面擬合範圍、粗差偵測分析及臺灣地區大地起伏擬合作業程序之實驗方法與步驟。. y. 研究成果與分析：. sit. Nat. 第五章. n. al. er. io. 檢視曲面方程式、曲面擬合範圍、粗差偵測分析及臺灣地區大地起伏擬合作業之實驗成果，並進行比較與分析。. Ch. 第六章. i n U. v. 結論與建議：. engchi. 對研究的成果進行討論，並提出結論及未來研究方向的建議。. 5.


(21) 文獻回顧. 第二章. 第一節大地起伏 GPS 水準測量技術的相關準則可以詳見美國大地測量局（National Geodetic Survey, NGS）中，Zilkoski et al.（2008）之文章內容雖然僅適用於美國地區，但內容可以做為研究臺灣 GPS 水準測量之參考。文中說明 GPS 水準測量技術在應用上是可以取代傳統水準測量；且在不超過十公里的範圍內利用 GPS 載波相位觀測，可以取得 2cm 精度的正高。NGS 的對 GPS 水準測量有三個基本規範（basic rules）、四個基本控制建議（basic control suggestions）及五個 GPS 水準測量獲取正高值的必要步驟。三個基本規範是分別是：1. 利用 GPS 測量取得的橢球高精度須依照 NGS 的規範；2. GPS 水準測量需使用 NGS 最新的混合大地水準面模型（hybrid geoid model）；3. 利用最新的國家垂直系統。四個基本控制建議分別是：1. 使用含有驗證過之 NAVD88 正高的點位作為 GPS 測站；2. 如果待測區域的一邊小於 20 公里，則該邊需要布設含有驗證過之 NAVD88 正高的水準點位； 3. 如果待測區域的一邊大於 20 公里，則該邊之含有驗證過之 NAVD88 正高的水準點位間的距離不可超過 20 公里；4. 如果待求區域在山區，則需要有區域內最高點及最低點水準點的高程資料。五個必要步驟分別是：1. GPS 測量使用 3-D、最小約制及最小二乘平差法；2. 利用第一個步驟的平差程序來偵測並移除所有粗差；3. 計算從最小約制平差計算後. 立. 政治大. ‧. ‧ 國. 學. sit. y. Nat. n. al. er. io. 的 GPS 水準測量正高值之差異量；4. 利用從步驟三以前的成果來決定哪些 NAVD88 控制點含有最小的殘差量可以視為有效點位，這個步驟決定了 GPS 水準測量獲取的正高之精度；5. 利用從步驟四以前的成果，藉由固定一個 NAD83 的控制點的經緯度坐標及所有有效的 NAVD88 點位高程來進行約制正高的平差計算。臺灣並沒有對於 GPS 水準測量的規範，本文中， NGS 對 GPS 水準測量的規範及計算程序等將有助於後續建立臺灣的 GPS 水準精度及施測規範。在實施 GPS 水準測量需要有一個大地起伏模型來提供大地起伏值資訊而大地起伏模型可以以重力方式或幾何方式來計算。. Ch. engchi. i n U. v. 國內利用重力法建立大地起伏模型的有黃金維等人（2013）採用以快速傅立葉變換為基礎的 stokes 積分理論計算的重力法大地起伏模型，並使用已實測 GPS 橢球高的一等水準點資料作為評估重力法大地起伏模型精度之用；若實測資料精度佳，則與重力大地起伏合併計算成一個混合型（hybrid）大地起伏模型。國內利用最佳化演算法擬合大地起伏的研究有：李尚訓（2010）、沈昱廷（2011）、簡子淩（2012）及陳佳菱（2012）。這些文章的研究區域皆在未改制前台中市，地勢平坦且含有橢球高之一等二級水準點點位分布密 7.

(22) 集，於是擬合之均方根誤差皆可達到±5cm 以下，惟這些文章中的點位之橢球高資料是以 GPS-RTK 方式獲取。一般來說，在 GPS 水準測量的誤差來源中，以 GPS 測量所得之橢球高含有最多誤差量，因此 GPS 橢球高施測的方式決定了 GPS 水準測量精度的品質。因此要得到更少誤差的擬合法大地起伏模型，GPS 橢球高應以靜態測量等精度較高的方法施測。. 第二節粗差偵測理論測量領域的可靠度理論（theory of reliability）主要是 Baarda（1968）從單個備選假設出發（僅考慮粗差而不考慮系統誤差），研究測量平差中的內部及外部可靠度。所謂內部可靠度代表平差系統本身能夠找到粗差的能力，而外部可靠度代表粗差對平差系統中參數估計的影響。而 Baarda 從可靠度理論所推導出檢驗粗差的資料探測法（Data Snooping），是利用單個一維的備選假設，將粗差歸入函數模型來判斷粗差。單個一維的備選假設代表觀測量中只有一個粗差的情形。. 立. 政治大. ‧. ‧ 國. 學. 資料探測法的基本假設是觀測量中僅含有一個粗差，若是觀測量中含有多個粗差，理論上可以透過多次一維的資料探測法逐步剔除粗差，這種方法又稱為迭代的資料探測法（iterated data snooping）。Schwarz et al. （1993）利用迭代的資料探測法偵測誤差，會較穩健估值法（ Robust estimation）來得有效率。並且反對過去一些對於粗差的判定標準，例如剔除殘差量最大的觀測量，或是剔除殘差量超過三倍中誤差的觀測量。. sit. y. Nat. n. al. er. io. Lehmann（2013）同樣對測量界常用的三倍中誤差準則提出質疑。作者將粗差視為隨機性粗差或系統性粗差來分別討論，這也等同於李德仁等人（2005）將粗差視為隨機模型的一部分或是函數模型的一部份，而這樣的差別使得粗差對機率分布函數的影響會不同。而不論將粗差視為隨機性或是系統性，三倍中誤差準則基本上忽略了背後所隱含的風險。三倍的中誤差，代表臨界值（critical value）為 3，臨界值的給定會連帶決定第 I 類誤差（type I decision error）及第 II 類誤差（type II decision error）的機率。而在統計假說上，每個觀測量是否為粗差都會建立一個備選假設，使得第 I 類誤差會隨著觀測量增加而倍增。於是作者評論說，在大量觀測量的情況下，僅僅使用三倍中誤差準則可能會錯判許多正確的資料。. Ch. engchi. i n U. v. 於是在粗差偵測上，一般所使用的三倍中誤差準則實際上並不是相當嚴謹，背後隱含的決策誤差（decision error）是使用者所忽略的。於是許多種的粗差偵測方法被開發出來，前面所提的迭代的資料探測法是其中之一，但其缺點在於：第一次計算資料探測法中，各個觀測量的標準化殘差中，標準化殘差最大的觀測量，可能並不是粗差，將其剃除會產生錯誤的判斷（李德仁，2005）。 8.

(23) 所以另外一種粗差偵測的方法不再去剔除粗差觀測量，而是減弱粗差對平差系統的影響。李德仁（1984）則是對粗差採用不同的處理方法，該方法將每次平差後的殘差或其他統計量，利用事先選擇的權函數計算下一次迭代中，各個觀測量的權。如果權函數適當，則含粗差觀測量的權將會在迭代計算中越來越小，達成定位粗差的效果。另外一種較廣為人知的方法是穩健估值法（Robust Estimation），Kubik et al.（1985）統整當代各種測量領域中常用的穩健估值法，最簡單的是最小範數法，還有 Hubers 方法及 Krarup 及 Kubik 等人提出的丹麥法（Danish Method）。. 第三節粒子群演算法粒子群演算法最早是由 Kennedy et al.(1995)提出了粒子群演算法最初步的概念：粒子群演算法介於基因演算法及進化演算法（evolutionary algorithm）中間。粒子群演算法的隨機性及函數適應值（fitness）的概念與進化演算法相似，而每次迭代更新個體最佳值及群體最佳值則類似基因演算法中個體交配（crossover）產生新個體的概念。. 立. 政治大. ‧ 國. 學. ‧. Van den Bergh (2001)對基本粒子群演算法的收斂性進行數學證明，研究證明基本粒子群演算法只能收斂到區域最小值（local minima），並不保證能夠收斂到全域最小值（global minima）。文中改良基本粒子群演算法中的使其能夠收斂至區域及全域最小值，而改良的方式有兩種：一種方法是在一般的粒子群依照正常搜尋過程行動外，透過一組持續不斷進行隨機搜. sit. y. Nat. n. al. er. io. 尋的獨立粒子群，來找到全域的最小值；另一種方法是粒子群找到區域最小值後，重新賦予粒子群新的隨機位置，再搜尋一次最小值。這兩種搜尋方式透過數學證明可以找到全域的最小值。. Ch. engchi. i n U. v. Sun et al.（2004）發表了基於量子力學的新形態粒子群演算法。該種演算法假設粒子是以量子力學的基礎在搜尋空間中搜尋最佳位置，實驗顯示這種改良算法相較於基本的粒子群演算法有更好的成果，並且有全域的搜尋的能力。孫俊等人（2011）證明量子行為粒子群演算法可以在一定的機率下收斂到全域的最小值。文中表示：在量子行為粒子群演算法中，固定維度的搜尋空間下，粒子群越大，全域最小值的收斂率越大，並且可以避免太早收斂到非全域最小值，但是運算速度將會降低。於是在粒子群體數量的控制上，應該要針對目標函數的複雜性來考慮數量的多或少：低維的函數應該用較少的粒子個數來增加運算速度；高維的函數應該增加粒子個數避免算法過早收斂到非全域最小值。該書除了討論量子行為粒子群演算法的收斂性外，也帶入了平均最好 9.

(24) 位置的概念（mean best position），強化了演算法的搜尋能力。該研究詳細的內容同時發表在 Xi et al.（2008）中。平均最好位置是每個粒子其最佳位置的平均數，而該位置的最大功用是避免粒子的太快收斂到目前的全域最佳位置，增加其他未位於全域最佳位置附近的粒子搜尋其他空間的機會。粒子群演算法在應用面上，Yetkin et al.（2011）在實地測量前先利用粒子群演算法選取最佳的 GPS 基線組合，能夠有效減少 GPS 測量所花費的成本；Mohammadzadeh et al.（2009）利用粒子群演算法對 fuzzy cost 函數求解最佳的函數值來自動萃取衛星影像中的道路；Yang et al.（2012）利用粒子群演算法選擇高光譜影像的波段，並利用最小二乘向量機對都市土地覆蓋物進行分類，分類精度高於使用所有波段資料或是利用 PCA （principal component analysis）及 LDA（linear discriminant analysis）方法篩選波段之成果。. 立. 政治大. ‧. ‧ 國. 學. n. er. io. sit. y. Nat. al. Ch. engchi. 10. i n U. v.

(25) 理論基礎. 第三章. 第一節大地起伏概念在高程系統中，主要有三個高程基準：大地水準面（geoid）、似大地水準面（quasi-geoid）及參考橢球面（ellipsoid）。大地水準面是正高（orthometric height）的起算基準，是一個最近似於平均海水面的重力等位面；似大地水準面是正常高（normal height）的起算基準，是以平均正常重力代替大地水準面所使用的實際重力值所定義出之曲面；參考橢球面則是橢球高（ellipsoid height）的起算基準，是以參考橢球體為依據之曲面。. 政治大. 大地起伏值（geoidal undulation）為地面上某一點其正高與橢球高之. 立. 學. H=h−N. (1). ‧. 其中：. ‧ 國. 差異量，如圖 1，三者的關係可以簡單的用式(1)表示。. y. Nat. er. io. sit. H=正高；N=大地起伏；h=橢球高；. 正高是沿地面點垂線方向到大地水準面之距離，而橢球高是地面點依法線. n. al. Ch. i n U. v. 方向到橢球面之長度。因為正高的垂線方向非直線且不垂直於橢球面，如. engchi. 圖 2，所以式(1)並不是等式；而應該修正為式(2)。 H≈h−N (2) 雖然式(1)並不夠嚴謹，但其誤差的影響 GPS 水準測量等應用來說可以接受。（Meyer et al, 2006）. 11.

(26) 政治大. 圖 1 正高、橢球高及大地起伏的關係圖. 立. ‧. ‧ 國. 學. n. er. io. sit. y. Nat. al. Ch. engchi. i n U. v. 圖 2 正高、橢球高及大地起伏的實際關係. 第二節大地起伏值之解算方法有幾種方法可以解算大地起伏值：(1) 以全球大地位模式解算大地起伏、 (2)Stokes 積分公式法、(3)衛星測高法及(4) 擬合法。（孔祥元等，2011）前面兩種方法屬於物理大地測量方法，後面兩者屬於幾何大地測量方法。. 12.

(27) 一、以全球大地位模式解算大地起伏以全球大地位模式解算大地起伏是利用實測的真實引力位與參考橢球體之正常引力位資料，代入 Bruns 公式計算大地起伏（Ardalan et al, 2001），如式(3)： N=. 𝑇 𝛾 (3). 其中： N 為大地起伏值；T 為擾動位；γ為正常重力值；. 政治大. 擾動位是指大地水準面上某一點之實測引力位與正常引力位之差值，如式. 立. 式(4)：. ‧ 國. 學 (4). ‧. 其中：. T=V−U. y. Nat. io. sit. V 為實測引力位；U 為正常引力位；. n. al. er. 經過推導，可以求得擾動位 T 的球諧函數泰勒級數展開式，如式(5)，並. i n U. v. 帶入式(3)，可以得到利用全球大地位模式解算大地起伏值之公式，如式 (6)：（陳國華，2004）. Ch. engchi. ∞. 𝑛. 𝑛=2. 𝑚=0. 𝐺𝑀 𝑎 𝑛 ̅ cos 𝑚𝜆 + 𝑆𝑛,𝑚 ̅ sin 𝑚𝜆)𝑃̅𝑛,𝑚 (sin 𝜙) T(r, θ, λ) = ∑ ( ) ∑ (𝐶𝑛,𝑚 𝑟 𝑟 (5) ∞. 𝑛. 𝑛=2. 𝑚=0. 𝐺𝑀 𝑎 𝑛 ̅ cos 𝑚𝜆 + 𝑆𝑛,𝑚 ̅ sin 𝑚𝜆)𝑃̅𝑛,𝑚 (sin 𝜙) N= ∑ ( ) ∑ (𝐶𝑛,𝑚 𝛾𝑟 𝑟 (6) 其中：. 13.

(28) G 為萬有引力常數；M 為地球質量；a 為橢球的長半徑；𝜙為地心緯度； ̅ 及𝑆𝑛,𝑚 ̅ 為完全正規化之引力位球諧函數； λ為經度；r 為地心距離；𝐶𝑛,𝑚 𝑃̅𝑛,𝑚 為完全正規化之 Legendre 函數；n,m 分為階數及次數。現今較常被使用的全球大地位模式是 EGM2008 模型，二、Stokes 積分公式法 Stokes 積分公式法是利用 Stokes 積分式計算大地水準面上一點擾動位，再利用 Bruns 公式計算大地起伏值，推導後如式(7)：（Heiskanen, et al, 1967）. N=. 立. 政治大 𝜎. (7). ‧ 國. 學. 其中：. 𝑅 ∫ ∫ Δ𝑔𝑆(𝜓)𝑑𝜎 4𝜋𝐺̅. ‧. 𝐺̅ 為全球重力平均值； σ 為單位球體；dσ 為單位球體上的面積元素； ψ 為球面弧距；Δg 為大地水準面上的重力異常；S(ψ)為 Stokes 函數。. er. io. sit. y. Nat. 三、衛星測高法. al. v i n Ch 射雷射波，經由接收地球表面反射的信號後，以雷達波之波速及訊號來回 engchi U 時間差計算地球表面到衛星的高度，本法適用於決定海水面之重力場與大 n. 衛星測高法原理是利用雷達測高儀（Radar Altimeter）向地球表面發. 地起伏。（趙嘉展等人，2010）之後可以利用式(8)計算大地起伏： h = N + H + ΔH + a + d (8) 其中： h 為橢球高；N 為大地起伏；H 為海洋地形；ΔH為海象影響量（海潮、固體潮、海流等）；a 為測高儀觀測量；d 為徑向軌道誤差。四、幾何擬合法. 14.

(29) 以幾何擬合法解算大地起伏值是利用一個數學模式來模擬大地起伏面，並以內插計算待求點之大地起伏值。數學模式可以為平面或曲面，一般常見的是採用二次曲面方程式，如式(9)： ξ = 𝑎0 + 𝑎1 x + 𝑎2 𝑦 + 𝑎3 𝑥𝑦 + 𝑎4 𝑥 2 + 𝑎5 𝑦 2 (9) 其中ξ為待求點之大地起伏值；𝑎0 ~𝑎5 為待定係數；x,y 則為待求點的平面坐標。利用曲面擬合法擬合的範圍是有所限制的，範圍過大會使得擬合精度變差。大地起伏是取決於地球內部密度變化及地形起伏等因素，因此可以. 治政誤差會越大，所以將大區域劃分為大小適中的區域有助於消除模型誤差的大立影響。（匡志威等人，2010）用一組曲面來精確確定大地起伏的分布；而區域面積越大所導致的是模型. ‧ 國. 學. 因此最理想的狀況是在合理的範圍內使用合理的曲面方程式擬合大地起伏值。而合理的擬合範圍會受到局部地形起伏變化所影響，地形起伏小的地. ‧. 區（如：平原），可以在大區域得到良好的成果；若地形含有丘嶺、山區. n. al. 一、可靠度理論. er. io. sit. y. Nat. 等，受到局部地形變化的影響，不能得到良好的擬合精度。. 第三節粗差偵測原理 iv. Ch. n U engchi. 傳統測量平差使用的最小二乘平差法，抗粗差的能力較不足夠，原因在於最小二乘平差法配賦偶然誤差的能力過於良好。在最小二乘平差中，含有較大誤差的觀測量因為誤差配賦的關係，反應在殘差量上的總是一小部分，使得使用者會漏判錯誤資料。於是 Baarda 所發展出的可靠度理論，是要在最小二乘平差的架構下，以統計假說檢定為基礎，探討平差系統尋找粗差的能力及粗差對平差系統的影響。平差系統尋找粗差的能力稱為內部可靠度，而粗差對平差系統的影響則是外部可靠度。要討論 Baarda 的可靠度理論必須從平差系統的幾何條件開始推導。先假設一個間接觀測平差方程式如式(10)： V=A∙X−l 15.

(30) (10) 其中，V 為觀測量改正數；A 為係數矩陣；l 為觀測量而 X 為估計參數。此為觀測方程式的函數模式，而觀測方程式的隨機模式如式(11)： 𝐷𝑙𝑙 = 𝜎0 2 𝑄𝑙𝑙 = 𝜎0 2 𝑃−1 (11) 其中𝐷𝑙𝑙 為觀測量的變方－協變方矩陣；𝜎0 為單位權中誤差；𝑄𝑙𝑙 為觀測量的協因數矩陣；P 為觀測量的權矩陣。而本觀測方程式的法方程式為式 (12)： 𝐴𝑇 𝑃𝐴𝑋 = 𝐴𝑇 𝑃𝑙 (12). 政治大. 由式(12)可以推導出估計參數矩陣，如式(13)：. 立. 𝑋 = 𝑁 −1 𝐴𝑇 𝑃𝑙. ‧ 國. 學. (13). ‧. 其中 N 代表法系數矩陣也就是式(12)的𝐴𝑇 𝑃𝐴。將式(13)代入式(10)可以得到式(14)：. al. y (14). er. io. sit. Nat. V = A𝑁 −1 𝐴𝑇 𝑃𝑙 − 𝑙 = (A𝑁 −1 𝐴𝑇 𝑃 − I)𝑙. v. n. 式中 I 代表單位矩陣。從式(14)可以得到改正數的協因數矩陣，如式(15)：. Ch. engchi. 𝑄𝑣𝑣 = 𝑃−1 − 𝐴𝑁 −1 𝐴𝑇. i n U. (15) 將式(15)的結果代回式(14)，可以得到式(16)： V = −𝑄𝑣𝑣 𝑃 ∙ 𝑙 (16) 式(16)說明了改正數會根據𝑄𝑣𝑣 𝑃及觀測量𝑙所決定，但是𝑄𝑣𝑣 𝑃是在未觀測前就可以計算的，理論上可以根據事前所設計的 A 矩陣及觀測量權矩陣 P 來推估觀測量誤差與改正數的關係，也因此𝑄𝑣𝑣 𝑃也被稱為平差的幾何條件（李德仁等人，2005）。在平差的幾何條件𝑄𝑣𝑣 𝑃矩陣中的對角線元素，代表各個觀測量的多餘觀測分量。當多餘觀測分量為 0 時，該觀測量為必要觀測量；而當多餘觀 16.

(31) 測分量為 1 時，代表該觀測量完全多餘，不參與平差計算。推導出平差的幾何條件後，接著可以建立一個統計檢定包含零假設（𝐻0 ）及一個備選假設（ 𝐻𝐴 ）及基於標準常態分佈（ standard normal distribution）下的檢定統計量（𝜔𝑖 ），也就是 Baarda 所提出的標準化殘差（standard residual），如式(17)： 𝜔𝑖 =. 𝑣𝑖 𝜎0 √𝑞𝑖𝑖 (17). 其中𝑣𝑖 為第 i 個觀測量的殘差；𝜎0 為先驗的單位權中誤差；𝑞𝑖𝑖 為𝑄𝑣𝑣 𝑃矩陣中第 i,i 個元素。在零假設下，標準化殘差是符合標準常態分佈的，如式 (18)：. 立. 治政𝜔 |𝐻 ~𝑁(0,1) 大 𝑖. 0. (18). ‧ 國. 學. ‧. 在可靠度理論中，內部可靠度的概念就是平差系統能夠找到粗差的下限值，下限值越小則內部可靠度越好。於是假設觀測量𝑙𝑖 含有粗差量∇𝑙𝑖 會使得檢驗統計量𝜔𝑖 的機率密度函數產生相應的位移量𝛿𝑖 （在統計學上稱之為非中心參數, non-centrality parameter），如圖 3，觀測量中存在粗差𝑒𝑔 ，使得原本的機率密度函數的平均值位移，可以接著計算平移量𝛿𝑖 ：. n. al. er. io. sit. y. Nat. ∇𝑣𝑖 = −𝑟𝑖 ∇𝑙𝑖. Ch. engchi U. v ni. (19). 其中∇𝑣𝑖 代表改正數因為粗差量∇𝑙𝑖 所產生的差異量，並藉由標準化殘差公式(17)計算位移量，如式(20)： 𝛿𝑖 =. −∇𝑣𝑖 𝜎𝑙𝑖 √𝑟𝑖. =. ∇𝑙𝑖 √𝑟 𝜎𝑙𝑖 𝑖 (20). 其中，𝑟𝑖 為第 i 個觀測量的多餘觀測分量。. 17.

(32) 圖 3 機率密度函數因為觀測量中的粗差而產生位移(Lehmann, 2013). 從式(20)可以推導出可發現粗差量的最下限值∇0 𝑙𝑖 ，也就是觀測量的內部可靠度，如式(21)： ∇0 𝑙𝑖 = 𝜎𝑙𝑖. 𝛿0 √𝑟𝑖. (21) 政治大 ∇ 𝑙 代表第 i 個觀測量可以發現粗差的下限值；而𝛿 代表零假設與備選假設立 0 𝑖. 0. ‧ 國. 學. 在給定的信心水準α（也是第 I 類誤差機率）及檢驗功效β（也是第 II 類誤差）下，最短的距離。Baarda 的建議是α=0.1%、β=80%則𝛿0 =4.13。在式(21)中，由於各個觀測量的中誤差不同，不能在同一個基準下比較每個. n. er. io. al. ∇0 𝑙𝑖 𝛿0 = 𝜎𝑙𝑖 √𝑟𝑖. sit. Nat. ′ 𝛿0,𝑖 =. y. ‧. 觀測量的內部可靠度，於是將各個觀測量的中誤差從式中除去，得到如式 (22)：. Ch. engchi. i n U. v. (22). ′ 上式𝛿0,𝑖 稱為可控性數值，可控性數值是一個倍數值，也就是觀測值𝑙𝑖 上的粗差至少為其中誤差的幾倍才可以在給定的機率（信心水準及檢驗功效）下被發現（李德仁等人，2005）。. 二、資料探測法 Baarda 的資料探測法，主要是在已知的單位權中誤差的前提之下，利用標準化殘差可以偵測平差系統中可能的粗差。資料探測法的基本架構是以常態分佈的統計檢定（Gauss Test），並在事前給定的信心水準α計算臨界值 c，來判斷個別觀測量的標準化殘差是否超過臨界值並依此判斷粗差。信心水準α與臨界值 c 的關係如式(23)（Lehmann, 2013）： 𝛼 c = Φ−1 (1 − ) 2 (23) 18.

(33) 其中，Φ為常態分佈 N(0,1)的累積分布函數。當臨界值為 3 時α=0.0027，是一般測量實務上常用的判定方法，即所謂的三倍中誤差準則。但是包含李德仁等人（2005）及 Lehmann(2013)都提出這個準則並不是相當完備，因為並未考慮到第 II 類誤差的機率（也就是(1-檢驗功效β)的機率），即錯誤資料誤判為正確資料的機率。前一小節所提到的非中心參數（式(20)）內涵其實已經包含檢驗功效，α=0.0027、c=3，非中心參數δ=4 的情況下，第 II 類誤差的機率為 16%，其計算的公式如式(24)（李德仁等人， 2005）： β(α, δ) = 1 −. 1 √2𝜋. ∫. 𝐾𝛼 −𝛿. 𝑢2. 𝑒 − 2 𝑑𝑢 +. −∞. 1 √2𝜋. ∫. −(𝐾𝛼 +𝛿). 𝑢2. 𝑒 − 2 𝑑𝑢. −∞. (24). 政治大. ‧ 國. 立. 學. 其中𝐾𝛼 代表在指定的信心水準下的臨界值。於是資料探測法需要同時考慮三個參數，使得使用者可以控制決策錯誤的機率在可接受的範圍，Baarda 是取α=0.1%（第 I 類誤差機率 0.1%）、β=80%（第 II 類誤差機率 20%）而δ=4.13。資料探測法的計算步驟如下：. ‧. 1. 依照目標函數計算最小二乘平差法。. sit. y. Nat. 2. 計算各個觀測量的標準化殘差。. io. al. er. 3. 依給定的信心水準與檢驗功效來決定非中心參數大小。. v. n. 4. 判斷各個觀測量的標準化殘差的絕對值是否超過非中心參數大小，並依此定位粗差。. Ch. engchi. i n U. 但資料探測法在統計假說上的備選假說是假設平差系統中僅只有一個粗差，在粗差偵測的實務上，是將標準化殘差絕對值最大的觀測量降權或是剔除。在實際的情況下，測量觀測量中的粗差觀測量大於 1 個的機率很高，進行一次資料探測法並不能完全定位出所有粗差觀測量。理論上資料探測法可以處理多個粗差的狀況，也就是使用迭代的資料探測法。操作方法是，在第一次迭代時，備選假設是觀測量中僅存在一個粗差，將標準化殘差的絕對值超過臨界值的一筆或多筆觀測量中，絕對值最大的觀測量剔除，也就是資料探測法，並進行下一差迭代。第二次迭代過程中，同樣基於備選假設是觀測量中僅存在一個粗差，同樣進行資料探測法剔除一個粗差後再進行下一次迭代。直到所有的觀測量的標準化殘差絕對值都低於臨界值後，迭代的資料探測法結束，同時也剔除了多個粗差。 19.

(34) 但是直接將可能的粗差剔除再進行計算的方法，理論上會有個盲點。在觀測量含有多個粗差的前提下，最小二乘平差的原理會使得各個觀測量的改正數受到多個粗差的影響，粗差並不一定會反應在改正數上；也有可能正常觀測量的改正數因為受到周圍粗差觀測量的影響而變大。如果因為粗差影響而使得改正數偏高的正常觀測量，其標準化殘差不小心超過了臨界值，在迭代的資料探測法中可能會被刪除，導致後續的參數估計產生錯誤。三、選擇權迭代法迭代的資料探測法每次迭代過程直接將可疑資料刪去的作法並不是相當理想，可能會誤刪正確資料，而解決的方法可以從觀測量的權值處理，也就是本節所謂的選擇權迭代法。. 政治大. 選擇權迭代法的基本作法是，根據每一次迭代計算後的改正數、標準化殘差或是其他統計量，依據事先選擇的權函數來計算下一次各個觀測量的權值。最簡單的權函數如式(25)： 𝑃𝑖𝑖𝑛+1 =. 學. ‧ 國. 立. 1 (𝜈𝑖𝑛 )2. ‧. (25). n. al. er. io. sit. y. Nat. 其中，𝑃𝑖𝑖𝑛+1代表第 n+1 次迭代時，權矩陣的第 ii 個元素；𝜈𝑖𝑛 代表第 n 次迭代時，第 i 個觀測量的改正數。這個方法是典型的利用平差後的改正數平方倒數評估觀測量的權，觀測量的改正數越大則權越小。在實務上，式 (25)分母應該加入常數項，以免改正數大小為零或幾乎為零時，計算機無法進行運算。. Ch. engchi. i n U. v. 丹麥法（Danish Method）是其中一種選擇權迭代法，設計的初衷是 Krarup 針對光束法平差的粗差偵測所設計的，其運算的基本方法為（Juhl, 1984）： 1. 在第一階段，按照基本的最小二乘平差法計算。由於原文是針對光束法平差，於是這個步驟通常需要進行 2 至 3 次的迭代計算。 2. 在第二階段，導入設計的權函數計算各個觀測量的權值，如式(26)： 𝑎. 𝑘 ∙ √𝑎𝑝 ∙ 𝑟 𝑝𝑖 = exp(−0.05 ( ) ) 𝜎0 (26) 其中，𝑝𝑖 代表第 i 個觀測量的權值；k 為常數，在第二及第三階段分別為 20.

(35) 1.0、0.6；ap 為先驗的權值，通常設定為 1；r 為殘差；𝜎0 為前一次迭代的單位權中誤差；a 為常數，在第二及第三階段分別為 4.4 及 6.0。 3. 在第三階段，同樣依照式(26)的權函數進行迭代計算，但由於參數設定的關係，在本階段只有少量的觀測量會被降權。丹麥法權函數的概念是，在第二階段時盡量嚴格的搜尋可能的粗差，可能會有 25%的殘差較大的觀測量會被降權；而在第三步驟時，由於參數的設定會使得僅有少數觀測量被降權（0.5%~2%）；理論上，含有粗差的觀測量會在整個過程中持續不斷地降權。不同於前面描述的迭代的資料探測法，是將每次超出臨界值的觀測量剔除，丹麥法計算一開始盡量將所有觀測量保留，而逐步將所有可疑的粗差降權，可以避免一開始錯誤的判斷而將正確觀測量剔除。. 政治大. 值得注意的是，包含丹麥法在內與其他有關大地測量或攝影測量有關的權函數，其函數的建立方式都是依長年的實驗數據所歸納，並沒有嚴謹的理論支持。. 立. ‧ 國. 學. 於是在運用選擇權迭代法時，應該針對不同的測量實務選擇合適的權函數；或是根據前人提出的方法與概念自行建立一套合適的權函數來使用。. Nat. y. ‧. 第四節粒子群演算法原理. er. io. sit. 一、基本粒子群演算法理論. 粒子群演算法（Particle Swarm Optimization）由 Kennedy 與 Eberhart 在 1995 年所提出，是一個基於群體智慧的演算法。粒子群演算法利用多個粒子組成一個粒子群，並在給定的搜尋空間（search space）中搜尋最佳. n. al. Ch. engchi. i n U. v. 解。在粒子群演算法中，所有粒子所在的位置都被視為該目標函數（objective function）的可能解，藉由在每次迭代步更新粒子的速度及位置，粒子群演算法試著找到該目標函數的最佳解。利用圖 4 可以簡單說明粒子群演算法搜尋最佳位置的概念。所有粒子為了使自身及群體在搜尋空間中是處於最佳的位置，藉由計算每次迭代次數本身的位置適應值（即函數值）及全體最佳位置適應值，並比較與上一次迭代次數的適應值來決定如何移動下一步。理論上如果計算成功，最終粒子群將會收斂到極端值（extrema）。基本粒子群演算法的缺陷在於它只能保證收斂到區域的極端值，並不能保證能夠找到全域的極端值。. 21.

(36) 政治大圖 4 粒子群演算法粒子搜尋極端值的過程立. ‧ 國. 學. 二、基本粒子群演算法算法. ‧. 假設在 d 維的搜尋空間中有 M 個粒子群，第 i 個粒子的位置及速率分. er. io. sit. y. Nat. 別表示為：𝑋𝑖,1 , 𝑋𝑖,2 , … 𝑋𝑖,𝑑 , i = 1,2, … , M 及𝑉𝑖,1 , 𝑉𝑖,2 , … 𝑉𝑖,𝑑 , i = 1,2, … , M。每個粒子在第 t 次的迭代次數則利用式(27)及式(28)更新其速率及位置：（Eberhart, ea al, 1998） 𝑉𝑖𝑑 (𝑡 + 1) = 𝑉𝑖,𝑑 (𝑡) + 𝑐1 ∙ 𝑟1,𝑖,𝑑 (𝑡) ∙ (𝑃𝑖,𝑑 (𝑡) − 𝑋𝑖,𝑑 (𝑡) + 𝑐2 ∙ 𝑟2,𝑖,𝑑 ∙ (𝐺𝑑 (𝑡) − 𝑋𝑖,𝑑 (𝑡))). n. al. Ch. engchi. i n U. v. (27). 𝑋𝑖,𝑑 (𝑡 + 1) = 𝑉𝑖,𝑑 (𝑡 + 1) + 𝑋𝑖,𝑑 (𝑡) (28) 其中𝑉𝑖𝑑 代表第 i 個粒子在 d 維搜尋空間中的速率；𝑃𝑖𝑑 代表每個粒子的最佳位置；𝐺𝑑 代表粒子群的最佳位置；𝑟1,𝑖,𝑑 及𝑟2,𝑖,𝑑 是介於 0 和 1 間的隨機數列；而𝑐1及𝑐2 則是學習因子，用來控制粒子在移動上對個體最佳位置及群體最佳位置的參考比例。在每次迭代中，利用式(29)來更新個體最佳位置：（J. Sun, et al., 2004）. 22.

(37) 𝑃𝑖,𝑑 (𝑡), 𝑃𝑖,𝑑 (𝑡 + 1) = { 𝑋𝑖,𝑑 (𝑡 + 1),. 𝑖𝑓 𝑓[𝑋𝑖,𝑑 (𝑡 + 1)] ≥ 𝑓[𝑃𝑖,𝑑 (𝑡)] 𝑖𝑓 𝑓[𝑋𝑖,𝑑 (𝑡 + 1)] < 𝑓[𝑃𝑖,𝑑 (𝑡)] (29). 接著利用式(30)找到粒子群中最好的粒子位置，並用式(31)來更新群體的最佳位置： g = arg min {𝑓[𝑃𝑖,𝑑 (𝑡)]} , g ∈ {1,2, … , M} 1≤𝑖≤𝑀. (30) 𝐺𝑑 (t) = 𝑃𝑔,𝑑 (𝑡) (31). 政治大. ‧ 國. 立. 學. 如果下一次迭代的新的群體最佳位置優於上一次的群體最佳位置，則更新之，否則維持原狀。演算法會持續計算直到迭代結束，如果計算結果是收斂的，各粒子將會集中到解空間的最佳位置上；如果計算結果不收斂，粒子會散布在各空間中，這時需要提高迭代次數或是增加粒子數量來解決。圖 5 為基本粒子群演算法的演算過程。迭代次數 For t=1 to n. 計算每個粒子目前的適應值. n. er. io. sit. y. Nat. 初始化粒子群. al. ‧. 演算法開始. Ch. 各粒子的適應值是否優於上一次. engchi. i n U 否. v. 粒子的位置維持不變. 是. 粒子更新到新的位置. 更新群體的最佳值. 利用速度及位置更新公式來更新粒子的速度與位置. 否. 是否達到最大迭代次數. 輸出成果. 是. 演算法結束. 圖 5 基本粒子群演算法演算的流程. 23.

(38) 三、量子行為粒子群演算法理論量子行為粒子群演算法（Quantum-behaved Particle Swarm Optimization, QPSO），是粒子群演算法的改良算法，由孫俊等人於 2004 年所發表。本演算法的發展重點在於改良基本粒子群演算法只能保證收斂在區域極端值的缺點，嘗試以量子力學的角度將粒子群演算法移動的方式機率化，從而使得算法能夠收斂到全域極端值上（global extrema）。QPSO 與基本粒子群演算法的最大差異是粒子搜尋的方式不同。基本粒子群演算法給每個粒子速度及方向，粒子以直線運動方式搜尋最佳位置；QPSO 中的粒子的運動狀態以波函數描述，波函數ψ給定的概率密度函數來確定粒子在某個時間某個位置出現的機率。（孫俊等人，2011）方式出現在整個搜尋空間中，擴展了粒子移動的範圍。 QPSO 在參數設定上除了粒子群數量、迭代次數及問題維度外，還有收縮擴張係數α（constraction-expansion coefficient）。根據式(34)，α可以做為粒子的平均最好位置與目前位置間的距離加權量；當α越大，代表下. 政治大. 立. ‧. ‧ 國. 學. 一步的粒子搜尋範圍越大。α參數在控制上可以用線性遞減或是固定值方式給值；如果採線性遞減方法，代表粒子的搜尋策略是從一開始的大範圍搜尋到末期的小範圍搜尋，線性遞減公式則如式(35)。一般而言，最大及最小的α通常設為 1.0 及 0.5，但會因問題的複雜性而有不同的給值。. sit. y. Nat. 四、量子行為粒子群演算法算法. n. al. er. io. 假設在 d 維的搜尋空間中有 M 個粒子群，第 i 個粒子的位置及速率分別表示為：𝑋𝑖,1 , 𝑋𝑖,2 , … 𝑋𝑖,𝑑 , i = 1,2, … , M 及𝑉𝑖,1 , 𝑉𝑖,2 , … 𝑉𝑖,𝑑 , i = 1,2, … , M。在 QPSO 中，每個粒子依照式(32)、式(33)及式(34)更新其下一次迭代的位置：（J. Sun, et al, 2005）. Ch. engchi. i n U. v. 𝑀. 1 mbest 𝑑 = ∑ 𝑃𝑖,𝑑 𝑀 𝑖=1. (32) 𝑝𝑖,𝑑 = 𝜑 ∗ 𝑃𝑖,𝑑 + (1 − 𝜑) ∗ 𝐺𝑑 , φ = rand(0,1) (33) 1 𝑋𝑖,𝑑 = 𝑝𝑖,𝑑 ± 𝛼 ∗ |𝑚𝑏𝑒𝑠𝑡𝑑 − 𝑋𝑖,𝑑 | ∗ ln ( ) , 𝑢 = 𝑟𝑎𝑛𝑑(0,1) 𝑢 (34). 24.

(39) α = (𝛼𝑚𝑎𝑥 − 𝛼𝑚𝑖𝑛 ) ∗ (𝑡𝑚𝑎𝑥 − 𝑡𝑛𝑜𝑤 ) ÷ 𝑡𝑚𝑎𝑥 + 𝛼𝑚𝑖𝑛 (35) 其中，mbest 𝑑 為在 d 維下所有粒子最好位置的平均位置；𝑃𝑖,𝑑 代表第 i 個粒子在 d 維空間中的最佳位置；𝐺𝑑 為群體的最佳位置；𝑝𝑖,𝑑 為在𝐺𝑑 與𝑃𝑖𝑑 之間的吸引子（local attractor）；φ為介於 0 和 1 之間常態分佈下的隨機數； 𝑢是另一個介於 0 與 1 之間常態分佈下的隨機數；𝛼為收縮擴張係數； 𝛼𝑚𝑎𝑥 與𝛼𝑚𝑖𝑛 分別為最大及最小的收縮擴張係數；𝑡𝑚𝑎𝑥 與𝑡𝑛𝑜𝑤 分別為最大的迭代次數及目前的迭代次數。圖 6 為 QPSO 的演算流程圖，而 QPSO 的詳細演算過程則臚列於下（J. Sun, et al., 2004）（J. Sun, et al, 2005）： Step1：初始化一組隨機位置的粒子群。. 政治大. Step2：利用式(32)在每一維度下決定該維的平均最好位置mbest 𝑑 。. 立. Step3：利用式(29)決定目前各個粒子的最佳位置𝑃𝑖,𝑑 。. ‧ 國. 學. Step4：利用式(30)及式(31)來決定各維度的群體的最佳位置𝐺𝑑 。. ‧. Step5：比較下一次迭代的群體最佳位置，如果優於上一次迭代的值，則取代之；否則維持原位置。. y. Nat. sit. Step6：利用式(33)計算一個介於𝐺𝑑 與𝑃𝑖𝑑 間一個隨機的吸引子。. n. al. er. io. Step7：利用 Step6 計算的吸引子與式(34)來計算每個粒子新的位置。. i n U. v. Step8：重複 Step2~Step7，直到滿足最大的迭代次數。. Ch. engchi. 25.

(40) 演算法開始. 迭代次數 For t=1 to n. 初始化粒子群. 計算每個粒子目前的適應值. 計算平均最好位置 mbest. 各粒子的適應值是否優於上一次. 否. 粒子的位置維持不變. 政治大是. 立. 粒子更新到新的位置. ‧ 國. 學. Nat. 計算隨機的吸引子. n. al. er. io. sit. y. ‧. 更新群體的最佳值. Ch. engchi. 計算每個粒子新的位置. 否. i n U. v. 是是否已達最大迭代次數. 輸出成果. 演算法結束. 圖 6 量子行為粒子群演算法演算的流程. 第五節最佳化權矩陣一、基本構想最佳化權矩陣是本文所提出的一個針對粗差偵測的構想，是透過最佳化演算法來計算各個觀測量最佳的權重，使得平差後的各個觀測量檢定統 26.

(41) 計量絕對值能夠最小化；而這樣的運算也能透過權重值與統計量來判斷何者為粗差觀測量，達到粗差定位的效果。在本章第二節所談論的粗差偵測方法中，選擇權迭代法是一種有效的方法；透過事前決定的權函數，依照每次迭代後的統計量對可疑的觀測量進行降權，達到粗差定位的效果；但是權函數並沒有一個統一的選擇標準。也就是不同的權函數找到的粗差可能並不相同，原因在於實務上的權函數常常是透過經驗法則所訂定，不同測量方法下使用權函數並不相同，針對大地測量網型平差的粗差偵測操作方法不會和攝影測量光束法平差的操作方法相同。本文對粗差偵測方法不完全使用選擇權迭代的想法，而是採用其概念：從調整觀測量的權重來定位粗差，並使用最佳化演算法來計算觀測量最佳的權重值。最佳化演算法的功能在於能夠找到目標函數的極端值，對於粗差偵測來說，目標函數可以設定為平差後的中誤差、改正數總合或是各個觀測量的標準化殘差；在最佳化演算法後，計算觀測量的最佳權重值使得目標函數能夠達到設定的極端值。在這樣的構想下，理論上一組觀測量會存在一組對於目標函數來說，最好的權矩陣；而這個最好的權矩陣理論上可以藉由權重值定位出粗差。. 立. 政治大. ‧. ‧ 國. 學. 二、演算方法. Nat. y. sit. n. al. er. io. 本小節討論最佳化權矩陣的實際演算方法，本演算法將針對曲面擬合大地起伏進行粗差偵測，而在進行粗差偵測前，必須先決定哪種觀測量是需要進行粗差偵測。以曲面擬合大地起伏中，觀測量包含平面坐標與垂直坐標，平面坐標是已知水準點的 TWD97 坐標；而垂直坐標是已知水準點正高與橢球高的差異量，即大地起伏。從內政部的一等水準測量作業規範來看，並未要求平面坐標的精度，但對於正高資料有高標準的品質要求，而水準點的平面坐標及橢球高主要透過衛星定位測量所得，從內政部（2003）公告之台灣一等二級水準網之水準測量、衛星定位測量及重力測量成果說明為例，可以得到正高的點位中誤差為 8.47 公厘、經緯度的標準偏差平均值皆為±1.9 公分，而橢球高的標準偏差高達±5.5 公分。從以上資料可以發現，在先驗的資料品質上來看，品質較差的是橢球高資料、其次為平面資料，品質最好的是正高資料（達釐米級）。所以在討論曲面擬合大地起伏的粗差偵測時，應該把重心放在橢球高上，因為該資料的品質最差，最有可能產生粗差。. Ch. engchi. i n U. v. 於是本研究的粗差偵測應該將重心放在橢球高的粗差偵測，但本研究在平差觀測量上是大地起伏，與正高有函數關係，不能單獨討論。大地起伏是正高與橢球高的差異量，所以基於誤差傳播定律，大地起伏的中誤差. 27.