使用波形重定提升Geosat/GM測高精度與重力異常：以台灣近海為例

全文

(1)國立交通大學. 土木工程學系碩士論文. 使用波形重定提升 Geosat/GM 測高精度與重力異常：以台灣近海為例 Improving Accuracy of Geosat/GM Satellite Altimetry and Gravity Anomaly by Waveform Retracking:a Case Study Around Taiwan. 研究生：劉祐廷指導教授：黃金維中華民國九十四年六月.

(2) 使用波形重定提升 Geosat/GM 測高精度與重力異常：以台灣近海為例 Improving Accuracy of Geosat/GM Satellite Altimetry and Gravity Anomaly by Waveform Retracking:a Case Study Around Taiwan. 研究生：劉祐廷. Student：Yu-Ting Liu. 指導教授：黃金維. Advisor：Dr. Cheinway Hwang. 國立交通大學土木工程學系碩士論文. A Thesis Submitted to Department of Civil Engineering College of Engineering National Chiao Tung Unversity in Partial Fulfillment to the Requirements for the Degree of Master In Civil Engineering June 2005. Hsinchu, Taiwan, Republic of China. 中華民國九十四年六月.

(3) 使用波形重定提升 Geosat/GM 測高精度與重力異常：以台灣近海為例研究生：劉祐廷. 指導教授：黃金維. 國立交通大學土木工程學系. 摘要. 此研究焦點在於使 Geosat/GM 衛星測高資料藉由波形重定提升精度，而實驗區域為台灣海域。與開放海域之測高波形比較，海域之中接近陸地之波形較為複雜。本研究使用改良型門檻值演算法定義閥門改正量，隨之可計算距離改正量與海水面高改正量。重定之成果是藉由大地起伏比較海水面高來評估，該大地起伏是由台灣周圍高解析度大地水準面模型而得。對全區比較，由重定改正之海水面高與大地起伏之間提升 8%。在離海岸線 10km 以內之海域比較則改善了 17%且可使用之海水面高數量增加。計算重力異常之方法為 inverse Vening Meinesz 公式與最小二乘配置法，兩種方法同時使用原始海水面高與重定改正後之海水面高做比較。由測高計算得到之重力異常與船載重力異常比較結果顯示，相較於原始海水面高來說，重定改正後之海水面高得到較好的結果，最好之範例可以達到改善 12%。而最好之範例為最小二乘配置法使用重定改正後之海水面高結合空載重力資料計算得到之重力異常。. i.

(4) Improving Accuracy of Geosat/GM Satellite Altimetry and Gravity Anomaly by Waveform Retracking: a Case Study Around Taiwan Student：Yuting Liu. Advisor：Dr. Cheinway Hwang. Institute of Civil Engineering National Chiao Tung University. Abstract This research focuses on improvement of Geosat/GM satellite altimeter data by waveform retracking and the experimental area is over Taiwan waters. In comparison to altimeter waveforms over open oceans, waveforms over waters close to land are complicated. This work uses an improved threshold algorithm of waveform retracking to determine gate correction and then to compute range correction and sea surface height (SSH) correction. The performance of retracking is assessed by comparing SSH with geoidal heights from a high resolution geoid model around Taiwan. Over the entire region of comparison, the agreement between retracking-corrected SSH and geoidal heights is improved by 8%. Over waters less than 10 km to the coasts, the improvement is 10% and the number of valid SSH increases. The methods of inverse Vening Meinesz formula and least square collocation are used to compute gravity anomalies from both raw SSH and retracking-corrected SSH. The results show that, in comparison to raw SSH, retracking-corrected SSH yield a better agreement between altimeter-derived. ii.

(5) gravity anomalies and shipborne gravity anomalies, with a best case of 12% improvement.. Least-squares collocation is used to compute gravity anomalies from a. combined set of retracking-corrected SSH and airborne gravity data, which agree best with the shipborne gravity anomalies in all cases.. iii.

(6) 誌謝. 求學的這段路程，如今終於告一段落，今後將以新鮮人的身分邁入這個社會，心中有些許不捨與期待，最後這兩年的碩士生活，讓我體驗了很多人生經驗，相信這對於往後我在社會上的日子有很大的助益。. 首先感謝吾師黃金維教授於課業上不吝指教，於老師身上學到許多求學問與待人處事的道理，在此向恩師敬上最崇高的敬意。感謝陳春盛教授、史天元教授與李振燾教授於課業上用心的指導。感謝口試委員張嘉強教授與崔國強教授對於本論文精闢的見解與寶貴的意見使本論文更臻完備。. 感謝郭金運老師與鄧曉麗博士於論文研究過程中提供建議與協助，使本論文能順利完成。感謝博士班學長姊成機、欣瑩、榮寬、豫麒、大綱、宇伸、自強、廷融、志敏、弘基、亘昶以及畢業的學長姐大雄、米粉、宜珊、BOSS、介嵐、勇者、阿福、福利與小支，你們在課業上的指導使我受益良多。感謝這兩年來一起學習的夥伴展鵬、貓哥、佩珊、粘、阿達、印松、世青，多虧有你們，在新竹的兩年過的多彩多姿。感謝研究室學弟阿伯的幫忙。還有其他平時鼓勵我、幫助我的好朋友們，除了謝謝還是謝謝。. 最後要感謝我的家人，你們平時願意聽我訴苦、給我意見，你們的支持是我最大的動力。. iv.

(7) 目. 錄. 中文摘要………………………………….…………………………………………. i. 英文摘要……………………………………………………………………………. ii. 誌謝…………………………………………………………………………………. iv. 目錄…………………………………………………………………………………. v. 表目錄………………………………………………………………………………. vii. 圖目錄………………………………………………………………………………. viii. 第一章. 緒論………………………………………………………………………. 1. 1-1. 研究動機…………………………………………………………………. 1. 1-2. 文獻回顧…………………………………………………………………. 2. 1-3. 研究方法…………………………………………………………………. 2. 1-4. 論文架構…………………………………………………………………. 3. 衛星測高原理與 Geosat 測高衛星介紹…………………………………. 4. 2-1. 衛星測高原理……………………………………………………………. 4. 2-2. 誤差與改正………………………………………………………………. 6. 2-2-1. 軌道誤差……………………………………………………………. 6. 2-2-2. 測高儀誤差…………………………………………………………. 6. 2-2-3. 大地水準面算起至瞬間海水面之距離……………………………. 8. 2-3. Geosat 測高衛星介紹……………………………………………………. 8. 2-4. 使用資料…………………………………………………………………. 9. 第二章. 第三章 3-1. 波形與波形重定…………………………………………………………. 16. 介紹波形………………………………………………………………. 16. v.

(8) 波形重定方法介紹………………………………………………………. 19. β參數演算法…………………………………………………….. 20. 3-2-2 OCOG 演算法………………………………………………………. 24. 3-2-3. 門檻值演算法………………………………………………………. 25. 3-2-4. 改良型門檻值演算法………………………………………………. 26. 3-3. 選擇波形重定方法………………………………………………………. 28. 3-4. 波形重定前後成果比較…………………………………………………. 31. 3-4-1. 單一軌跡比較………………………………………………………. 31. 3-4-2. 全區比較……………………………………………………………. 35. 重力異常計算……………………………………………………………. 41. 計算重力異常方法………………………………………………………. 41. 4-1-1. 高斯濾波……………………………………………………………. 44. 4-1-2. 去除-回復法………………………………………………………. 45. 4-1-3. IVM………………………………………………………………. 47. 4-1-4. 最小二乘配置法……………………………………………………. 48. 成果分析…………………………………………………………………. 52. 結論與建議………………………………………………………………. 72. 參考文獻……………………………………………………………………………. 74. 作者簡歷……………………………………………………………………………. 77. 3-2. 3-2-1. 第四章 4-1. 4-2 第五章. vi.

(9) 表. 目. 錄. 表 2-1 GDRs 光碟資訊……………………………………………………………. 10. 表 2-2 GDRs 資料紀錄格式………………………………………………………. 11. 表 2-3 WDR 資料記錄格式………………………………………………………. 15. 表 3-1. 測高衛星資訊……………………………………………………………. 18. 表 3-2. 演算法解算波形成功率統計表…………………………………………. 30. 表 3-3. 改正後海水面高與大地水準面精度統計表……………………………. 31. 表 3-4. 不同距離對波形重定前後精度比較……………………………………. 39. 表 3-5. 全部區域分成不同區域對波形重定前後精度比較……………………. 40. 表 3-6. 沿岸 20 公里分成不同區域對波形重定前後精度比較…………………. 40. 表 3-7. 沿岸 10 公里分成不同區域對波形重定前後精度比較…………………. 40. 表 4-1. 以 IVM 計算波形重定前後重力異常與船測重力異常之精度分析表…. 53. 表 4-2. 以最小二乘配置法計算波形重定前後重力異常與船測重力異常之精度分. 析表………………………………………………………………………………… 表 4-3. 54. 以最小二乘配置法加空載計算波形重定前後重力異常與船測重力異常之. 精度分析表…………………………………………………………………………. vii. 54.

(10) 圖. 目. 錄. 圖 2-1. 衛星測高幾何關係圖………………………………………………………. 5. 圖 3-1. 波形概要圖………………………………………………………………. 17. 圖 3-2. 脈衝接觸海水面之狀況示意圖…………………………………………. 17. 圖 3-3. Geosat/GM 海洋波形………………………………………………………. 19. 圖 3-4. Geosat 靠近陸地波形………………………………………………………. 20. 圖 3-5. β-5 參數示意圖……………………………………………………………. 22. 圖 3-6. β-9 參數示意圖……………………………………………………………. 23. 圖 3-7. 尖錐波形示意圖…………………………………………………………. 23. 圖 3-8. OCOG 示意圖………………………………………………………. 25. 圖 3-9. 改良型門檻值演算法流程圖……………………………………………. 28. 圖 3-10. Geosat/GM 衛星 85206 地面軌跡圖……………………………………. 29. 圖 3-11. 改正後海水面高與大地水準面關係圖…………………………………. 30. 圖 3-12. Geosat/GM 衛星 85282 地面軌跡圖……………………………………. 32. 圖 3-13. 由 85282 軌跡取波形 1-20，左上角之數字表示圖 3-9 箭頭順序……. 33. 圖 3-14. 由 85282 軌跡取波形 21-40，左上角之數字表示圖 3-9 箭頭順序……. 34. 圖 3-15. 85282 波形重定前後沿軌跡之海水面高比較…………………………. 34. 圖 3-16. 實驗區域與衛星地面軌跡圖……………………………………………. 36. 圖 3-17. 原始海水面高與大地水準面差值圖……………………………………. 37. 圖 3-18. 波形重定後海水面高與大地水準面差值圖……………………………. 38. 圖 4-1. 以 IVM 計算重力異常流程圖……………………………………………. 42. 圖 4-2. 以最小二乘配置法計算重力異常流程圖………………………………. 43. 圖 4-3. 高斯濾波前後示意圖……………………………………………………. 45. viii.

(11) 圖 4-4. 長波長大地起伏等高線圖………………………………………………. 50. 圖 4-5. 長波長重力異常等值線圖………………………………………………. 51. 圖 4-6. 以最小二乘配置法計算測高資料波形重定前之重力異常等值線圖…. 55. 圖 4-7. 以最小二乘配置法計算測高資料波形重定後之重力異常等值線圖…. 56. 圖 4-8. 以 IVM 計算波形重定前重力異常分布圖………………………………. 57. 圖 4-9. 以 IVM 計算波形重定後重力異常分布圖………………………………. 58. 圖 4-10. 以最小二乘配置法計算波形重定前重力異常分布圖…………………. 59. 圖 4-11. 以最小二乘配置法計算波形重定後重力異常分布圖…………………. 60. 圖 4-12. 以最小二乘配置法加空載計算波形重定前重力異常分布圖…………. 61. 圖 4-13. 以最小二乘配置法加空載計算波形重定後重力異常分布圖…………. 62. 圖 4-14 IVM 計算波形重定前重力異常與船載重力異常差值圖………………. 63. 圖 4-15 IVM 計算波形重定後重力異常與船載重力異常差值圖………………. 64. 圖 4-16. 最小二乘配置法計算波形重定前重力異常與船載重力異常差值圖…. 65. 圖 4-17. 最小二乘配置法計算波形重定後重力異常與船載重力異常差值圖…. 66. 圖 4-18. 最小二乘配置法加空載計算波形重定前重力異常與船載重力異常差值. 圖…………………………………………………………………………………… 圖 4-19. 67. 最小二乘配置法加空載計算波形重定後重力異常與船載重力異常差值. 圖……………………………………………………………………………………. 68. 圖 4-20 以 IVM 計算之重力異常與船測重力異常差值過大分布圖……………. 69. 圖 4-21. 以最小二乘配置法計算之重力異常與船測重力異常差值過大分布圖. 70. 圖 4-22. 以最小二乘配置法加空載計算之重力異常與船測重力異常差值過大分布. 圖……………………………………………………………………………………. ix. 71.

(12) 第一章. 緒論. 1-1 研究動機. 衛星測高（satellite altimetry）發展迄今已有 30 多年，針對觀測量進行改正方法已臻成熟，如海潮改正模式、固體潮改正模式、電離層改正、乾濕對流層改正等。而衛星測高主要觀測量是藉由衛星之精確軌道以及衛星至天底點瞬間海水面（at-nadir instantaneous sea surface）垂直距離而得到之瞬間海水面（instantaneous sea surface），最後再經由地球物理改正以及環境改正後得到海水面高（sea surface height, SSH）。. 台灣地處歐亞大陸板塊及菲律賓板塊交界處，西邊靠近中國大陸之海底地形為大陸棚，地形較為平坦且深度為兩百公尺以內；而東邊靠近太平洋處為一海溝，海底地形於離岸十多公里處則開始劇烈變化，海洋深度可到四公里左右。又以台灣西岸為沙岸；東岸為沿岸之特殊地形，而且台灣本島附近還有澎湖、綠島、蘭嶼等小島，所以台灣及鄰域為一具有豐富實驗性之區域。. 於深海中，海水面可以近似為一平面，所以經海水面反射而回傳之脈衝波形（waveform）非常良好，故測距精度較高，但是當衛星接近陸地約 20 公里時，其足跡（footprint）會由於接觸到陸地或是較為複雜之海水面而導致回傳之脈衝波形被污染，進而影響衛星測距之精度，所以藉由波形重定（waveform retracking）求得近岸衛星測高改正後之海水面高。本研究除了以台灣最新大地水準面為判斷改正後海水面高精度是否提升，也利用計算重力異常來作為判斷波形重定前後是否對於精度上有所提升。. 1.

(13) 1-2 文獻回顧. 自從 Brown （1977）提出平均脈衝回傳理論後，就有許多人針對測高衛星回傳之波形進行研究並發展具有物理意義或是數學統計意義之演算法來密合波形，如 Wingham et al （1986）發展之 OCOG （off center of gravity）演算法、Martin et al. （1983）發展之β-參數演算法、以及 Zwally （1996）所發展之門檻值（threshold）演算法等，都是針對冰面反射之雷達波形而衍生。陳松安（2001）使用β-5 參數演算法於 ERS-1 測高衛星上，彭敏峰（2003）針對 TOPEX/POSEIDON 測高衛星所經過大陸內陸湖泊之軌跡（track），也使用β-5 參數演算法對其波形進行處理。 Deng （2004）對於 ERS-1 及 ERS-2 測高衛星資料進行波形重定，並且於台灣沿海以及澳洲沿岸得到非常好之改善。. 重力異常對於大地應用上非常重要，可以以重力異常推估垂線偏差，也可以計算大地起伏，吾人以海水面高轉換重力異常，於高豫麒（1997）使用一維快速傅立葉轉換（1D Fast Fourier Transform, 1D FFT）計算全球 2 ′ × 2 ′ 網格重力異常， Hwang（1998）的 inverse Vening Meinesz formula （IVM），使用核函數（kernel function）梯度將垂線偏差轉換為重力異常，張榮傑（2001）也使用了 IVM，將東西、南北垂線分量轉換為重力異常。. 1-3 研究方法. 本研究取得 Geosat/GM 之 geophysical data records （GDRs）及 waveform data record （WDR）資料，由於 GDRs 與 WDR 所記錄之時間格式不同，藉由 GDRs 與 WDR 時間換算公式，可以將 WDR 疊合到 GDRs 上，並且藉由重定 WDR 波形得到時間改正量隨之轉換為距離改正量，最後可針對 GDRs 所提供海水面高進行. 2.

(14) 改正而得到改正後海水面高（corrected SSH）。取得改正後海水面高後，將與台灣最新大地水準面模型（geoid model）進行比較，並且以不同距離、區域作為分析準則。. 得到改正後海水面高，則分別利用 IVM 與最小二乘配置法（least squares collocation, LSC）求得重力異常，並且於船載實測重力做精度比較分析。. 1-4 論文架構. 本論文共分五章，各章節內容為：第一章：說明研究動機、使用方法、文獻回顧與章節介紹。第二章：介紹衛星測高原理與 Geosat/GM 資料。第三章：詳談波形重定原理與波形重定後得到之改正後海水面高與大地水準面比較及分析。第四章：介紹所使用之計算重力異常方法，可分為逆范寧梅尼茲公式與最小二乘配置法兩種，並由計算得之重力異常分析其結果。第五章：對於本研究做結論與建議。. 3.

(15) 第二章衛星測高原理與 Geosat 測高衛星介紹. 2-1 衛星測高原理. 衛星測高基本原理為測高衛星上設置一雷達測高儀（radar altimeter）並向海水面發射微波脈衝（microwave pulses），經由海水面反射脈衝，然後由衛星上天線接收此一反射脈衝，經由計算發射至接收之時間差，即可計算出衛星至海水面之距離（Seeber, 1993）。. R = c×. Δt 2. （2-1）. 其中 R 是衛星至海水面之距離，即測高衛星量測距離， c 是光速， Δt 是衛星發射脈衝至接收訊號中間經過之時間。衛星測高基本觀測量為海水面高：. SSH = h − R. （2-2）. 其中 h 為衛星高度值，為該時刻衛星至參考橢球面之距離。. 衛星測高詳細之幾何關係可以用下列觀測方程式表達（參考圖 2.1）：. h = N + H + ΔH + R + d. （2-3）. N 為大地起伏（geoidal height），由橢球至大地水準面（geoid）之距離（Heiskanen and Moritz, 1967）。 H 為海面地形（sea surface topography, SST），由平均海水面高（mean sea surface. height, mssh）起算至大地水準面之距離，該值可以達到 1-2 公尺。. 4.

(16) ΔH 為瞬間效應，如潮汐變化、電離層效應、對流層效應等。 R 為衛星測高所得到之距離。. d 為計算軌道與真實軌道差。. 圖 2-1 衛星測高幾何關係圖. 5.

(17) 2-2 誤差與改正. 衛星測高誤差主要來自三個方面（Seeber, 1993），這三方面都關係到某程度之誤差：（1）軌道誤差：計算軌道與真實軌道差值。（2）測高儀誤差：沿著訊號傳播路徑之影響。（3）大地水準面起算至瞬間海水面之差值。. 2-2-1 軌道誤差. 軌道誤差主要是由下列原因所引起：－使用有限解析度與精度之地球重力場來計算軌道。－追蹤站座標之誤差。－追蹤系統精度上之限制。－軌道計算之模式錯誤。. 最主要之影響通常是來自地球重力場，因此需要發展適合之重力場模型。其次為追蹤系統之影響，以操作目的來說，全天候與輻射追蹤系統為最佳，如都卜勒（Doppler）追蹤系統。即使在最順利的情況與最新追蹤系統下，其餘兩項軌道誤差還是比測高儀測距精度大，因此需要更高精度之軌道計算方式。. 2-2-2 測高儀誤差. 沿著訊號路徑影響可分為：－儀器誤差。. 6.

(18) －訊號傳播誤差。最主要儀器誤差影響列在下面：－衛星質量中心與雷達天線相位中心之距離。－測高衛星電子傳播延遲。－量測系統時間誤差。. 這些影響在製造測高衛星時，可以被降至最低，所有儀器誤差影響在精密量測檢定場時即可求得並定義測高衛星之率定。對 SEASAT-1 來說，其率定值與零並無顯著不同（Kolenkiewicz et al, 1982）。另一個儀器誤差為天底誤差（nadir error），即脈衝光束與垂直方向之差，此受脈衝長度與波寬影響，因為足跡為脈衝長度所定義，所以可以忽略脈衝長度對天底誤差之影響。. 訊號傳播誤差是由電離層與對流層對於傳播速度之影響，訊號頻率 14GHz 於電離層折射影響下約 5 至 20 公分，而當衛星上搭載兩種頻率之儀器時，則可模擬電離層折射影響，如 TOPEX 衛星上裝有雙頻測高儀（13.6GHz 與 5.3GHz）。而垂直量測於對流層折射影響下大約 2.3 公尺，當測高衛星搭載輻射計沿軌跡同時觀測水蒸氣與輻射資訊時，所有對流層影響都可以用適當反射模型求至公分級精度，如 SEASAT、ERS-1、TOPEX 都可達到該精度，而 GEOSAT 衛星沒有搭載輻射計，所以其對流層影響之精度由氣象衛星觀測全球水蒸氣場所提供，精度大約為 2 公分（Cheney et al. 1991）。. 現實海洋狀態（sea state）影響也被包含在傳播誤差之中。反射區域大小依賴著海水表面之粗糙程度，粗糙程度越大造成之反射區域也越大。由於脈衝能量由波頂回傳較早，所以回傳波形前緣斜率會隨著波高增加而會呈現衰敗情形。. 7.

(19) 2-2-3 大地水準面算起至瞬間海水面之距離. 大地水準面起算至瞬間海水面之距離可以分為包含固定時間 H 與變動時間 ΔH ，海水面產生變化主要原因為海潮及大氣壓力效應，測高儀觀測時已經將隨著. 時間變化部分以平均海水面改正。. 全球海水面變形取自於大氣壓力，約可達到 10 至 50 公分，地區性影響則只有數公分並於測高儀觀測改正時被忽略（Monka, 1984），如果測高儀精度持續增加，則需要考慮地區性大氣壓力影響。. 海水面主要變形部分是由海潮（ocean tides）引起，於開放海域，海潮引起之振幅可以達到約 1 公尺，當靠近海岸或淺灘時，其振幅更可以達到 2 公尺，以某些特定例子來說會振幅可能更大。許多全球海潮改正模式精度可以達到 10 公分，甚至有模型可達到精度 3 公分（Wagner, 1991）。而固體潮（solid earth tides）會引起地球本體變形導致高度上變化，最多可達 50 公分，而目前固體潮改正之精度約可達到 1 公分。. 2-3 Geosat 測高衛星介紹. 美國海軍於 1985 年 3 月 12 日發射 Geosat（geodetic satellite）衛星，飛行計畫預計 5 年，1985 年 3 月 30 日開始接收測高衛星資料，1990 年 1 月 5 日結束任務。 Geosat 測距儀觀測精度為±3.5 公分（Seeber, 1993），衛星傾角為 108 度，飛行高度大約為 800 公里，衛星地面軌跡大約每秒 6-7 公里，測高儀所發射之微波寬度（microwave beamwidth）為 2 度，頻率為 13.5GHz，每秒發射 1020 個脈衝（1020Hz），由衛星上設置之 AGC（automatic gain control）接收自海水面回傳之. 8.

(20) 脈衝能量，並將每 100 個回傳能量累加平均後始紀錄其平均值，進而得到每秒 10 筆（10Hz）資料。Geosat 依照時間先後共執行兩個任務，分別為大地任務（geodetic mission, GM）以及重複軌道任務（exact repeat mission, ERM）：. 大地任務執行時間由 1985 年 3 月 31 日至 1986 年 9 月 30 止，主要目的在於求得高解析度之海洋大地水準面，該任務資料於初期並不公開，直到 1992 年 6 月才公佈所有資料。大地任務執行期間，其衛星軌跡非常密集並採取不重複軌道觀測。於緯度 60 度，衛星軌跡間距為 2-3 公里，赤道附近衛星軌跡與下一點之間距大約 10 公里，而赤道附近之地面相鄰軌跡（cross track）平均間距為 4 公里，因此可以提供精確之海洋重力場以及解析度。. 大地任務終止之後，衛星即調動軌道為每 17 天重複一次地面軌跡，並且與 SEASAT 衛星之地面軌跡非常相符，稱為重複軌道任務。重複軌道任務為不保密並且可自由使用於大地測量以及海洋學等。ERM 之相鄰地面軌跡於緯度 60 度時大約為 75 公里，遠比大地任務為寬。重複軌道任務期間大約三年，其資料廣泛地應用於各海洋科學，如海面變化。. 2-4 Geosat/GM 資料. 本研究所使用資料為 Geosat/GM，包含 GDRs 部分及 WDR 部分。GDRs 為美國國家海洋及大氣管理局（national oceanic and atmospheric administration, NOAA）提供，一共 10 片光碟（參考表 2-1），儲存資料以一天為單位，分別是 1-4 片光碟紀錄 Geosat/GM 資料，自 1985 年第 90 天開始至 1986 年第 273 天止，共 549 天； 5-10 片光碟紀錄 Geosat/ERM 資料，自 1986 年第 312 天開始至 1989 年第 364 天止，. 9.

(21) 共 1107 天。其中 Geosat/GM 之儲存格式（參考表 2-2），為每秒資料記錄有 34 個參數、由 78 個位元（byte）儲存，而 GDRs 地面軌跡資料本身為不完整覆蓋，所以一天大約只有 5 萬筆資料。. 表 2-1 GDRs 光碟資訊光碟. 任務. 開始時間. 結束時間. 1. GM. 1985-090. 1985-227. 2. GM. 1985-228. 1985-365. 3. GM. 1986-001. 1986-138. 4. GM. 1986-139. 1986-273. 5. ERM. 1986-312. 1987-082. 6. ERM. 1987-083. 1987-235. 7. ERM. 1987-236. 1988-040. 8. ERM. 1988-041. 1988-210. 9. ERM. 1988-211. 1989-031. 10. ERM. 1989-032. 1989-364. 10.

(22) 表 2-2 GDRs 資料紀錄格式項目. 參數. 單位. 位元. 描述. 1. UTC. sec. 4. UTC time since 1985/01/01.. 2. UTC. microsec. 4. UTC time, microseconds part.. 3. LAT. microdeg. 4. Latitude, microdegrees N.. 4. LON. microdeg. 4. Longitude, microdegrees E. JGM-3 orbit relative to reference ellipsoid. 5. ORB. ㎜. 4 ae = 6378136.3 m; 1/f = 298.257. 1-second average sea height relative to. 6. H. ㎝. 2 reference ellipsoid. Standard deviation of the 10/sec values. 7. SIG_H. ㎝. 2 about the 1-second H. Mean sea surface height from the Ohio State. 8. MSSH. ㎝. 2 MSS95 model. 10/sec sea height values. To derive time. 9-18. H1-H10. ㎝. 10*2 tags.. 19. SWH. ㎝. 2. Significant wave height Wind speed at 10m height,from Freilich and. 20. WS. ㎝/sec. 2 Challenor（1994）model.. 21. SIG_0. 0.01 dB. 2. 22. SSB. ㎜. 2. Sigma naught,radar backscatter coefficient Sea state bias derived by Gaspar, Ogor, and Hamdaoui（1996）.. 23. L_TID. ㎜. 2. Load tide from CSR 3.0 model.. 24. FLAGS. －. 2. See table below for flag bit definitions.. 11.

(23) 25. H_OFF. m. S_TID. ㎜. 2. H offset to be added to all heights over land. Solid tide from T/P “TIDPOT” algorithm,. 26. 2 based on Cartwright & Edden（1973）. 27. O_TID. ㎜. 2. 28. WET_NCEP. ㎜. 2. Ocean tide from CSR 3.0 model. Wet troposphere correction from NCEP/NCAR reanalysis model. Wet troposphere correction from NASA. 29. WET_NVAP. ㎜. 2 NVAP climatology. Dry troposphere correction from. 30. DRY_NCEP. ㎜. 2 NCEP/NCAR reanalysis model.. 31. IONO. ㎜. 2. Ionosphere correction from IRI95 model. Wet. 32. WET_T/S. ㎜. troposphere. correction. from. 2 TOVS/SSMI observations. Dry troposphere correction from ECMWF. 33. DRY_ECMWF. ㎜. 2 model.. 34. ATT Total. 0.01 deg. 2. Attitude（Space craft off-nadir orientation）.. 78. Bytes. 12.

(24) 表 2-2 中，第 9 項至第 18 項為海水面高，還需要加入其他地球物理改正量，才能應用在計算全球重力異常上：. H = H − (WET _ NCEP + DRY _ NCEP + IONO + O _ TID + S _ TID + L _ TID + SSB + IB + glo _ ib + hcal + uso) ÷ 10. （2-4）. 其中 H ：經地球物理改正後之海水面高。。 WET _ NCEP ：濕對流層改正。使用 NCEP/NCAR 模式（Kalnay et al., 1996）。 DRY _ NCEP ：乾對流層改正。使用 NCEP/NCAR 模式（Kalnay et al., 1996）. IONO ：電離層改正。 O _ TID ：海潮改正。為 CSR3.0 改正模式。 S _ TID ：固體潮改正。 L _ TID ：海潮改正。. SSB ：海洋狀態偏差改正。 IB ：地方性逆氣壓改正。. glo _ ib ：全球性逆氣壓改正。. hcal ：為震盪器內部因為溫度變化引起零點飄移以及儀器內部率定改正。. uso ：測高儀頻率因時間而變化之改正。上述所有改正，前面七項為主要改正，都紀錄在 GDRs 裡面，後面三項為次要改正，其變化都非常緩慢，而第八項改正 IB 是由乾對流層計算而得：. P(mbar) =. − DRY_NCEP(mm) (2.277 × (1.0 + 0.0026 × cos(2.0 × φ ))). IB(mm) = −9.948 × (P(mbar) − 1013.3). 13. （2-5）. （2-6）.

(25) 其中 φ 是衛星所在之緯度。. WDR 是由 Johns Hopkins University 應用物理研究室（ applied physics laboratory）提供，儲存資料以半日為單位儲存，所以一天有兩個檔案，其資料儲存格式（參考表 2-3）為每筆資料 660 位元，一天大約近 9 萬筆資料，每筆資料記錄著時間以及 10Hz 波形。. GDRs 時間紀錄格式為協調世界時（Coordinate Universal Time, UTC），自 1985 年 1 月 1 日開始，與當日資料起始之 FC0 ，而 WDR 時間紀錄格式為 FC （frame count），由下列公式計算而得：. FC = 32 × MFC + mFC. （2-7）. 其中 MFC （major frame count）與 mFC （minor frame count）由檔案直接讀出，本研究是將 GDRs 裡提供之協調世界值轉換並加上當日起始之 FC0 後即得到 FC_GDR ，可由此與 WDR 提供之 FC 相比，如 FC_GDR 與 FC 相同則可説是同一. 筆資料，下面為 FC_GDR 計算公式：. ⎛ ⎞ s85 - epoch_s85 - 0.005 FC_GDR = int ⎜⎜ FC0 + − 4.4 ⎟⎟ sec_per_fc ⎝ ⎠. 其中 FC0 為 GDRs 提供每日起始之 FC 。. s85 為 GDRs 提供之協調世界時。 epoch_s85 為 GDRs 提供每日起始之協調世界時。 sec_per_fc 為每筆 FC 所佔之時間。. 14. （2.8）.

(26) 表 2-3 WDR 資料記錄格式項目. 參數. 數據範圍. 型式. 開始位元. 1. MFC. 0Æ（224-1）. R3. 1. 2. mFC. 0Æ31. R1. 4. 3. Mode. 0Æ（230-1）. R4. 5. 4. Flagword. 0Æ（230-1）. R4. 9. 5. WS（1,1）. 0Æ255. R1. 13. 6. WS（2,1）. 0Æ255. R1. 14. ···. ···. ···. ···. ···. 634. WS（63,10）. 0Æ255. R1. 642. 635. SCALE（1）. 1,2,4. R1. 643. ···. ···. ···. ···. ···. 644. SCALE（10）. 1,2,4. R1. 652. 645. SPARES. NA. 8（X）. 653. Total. 660 Bytes. 15.

(27) 第三章波形與波形重定. 3-1 波形. 波形（圖 3-1），分為三部分：熱雜訊（thermal noise）部份、上升區部份與下降區部份，上升區通常以前緣（leading edge）稱之。波形是由測高儀天線發射微波脈衝，經海水面反射脈衝能量後，由測高儀上裝置 AGC 以不連續方式接收並以值域 0-255 儲存，將接收到之值連線則形成波形，圖 3-2 展示脈衝接觸海水面與波形形成相對應之情況，圖上方三個為側視圖、下方三個為俯視圖，時間 Ta 為脈衝剛接觸至海水面，其接觸範圍為一個點，對應至波形上則為能量正要上升，而時間 Tb 為脈衝接觸海水面一段時間，其接觸範圍為一個面，該面稱為足跡，對應至波形上則為接收能量至最大值，最後時間 Tc 為時間 Tb 之後，脈衝能量開始下降而與海水面接觸為一個環狀面。. 表 3-1 為所有測高衛星資訊，以 Geosat 為例，測高衛星距離觀測為發射時間至回傳脈衝最大振幅之一半，於圖 3-1 中即為上升區之一半，通常稱之為前緣中點，而電腦定該時間點為波形第 30.5 個閥門（gate），稱為預設閥門（tracking gate），參考圖 3-3，為 Geosat/GM 於深海中之波形，以圖中虛線表示預設閥門，前後各推 30 個閥門後即可得到一個波形，所以其閥門個數共 60 個，閥門與閥門之時間間隔為 3.125 奈秒（nanosecond, ns），可由此時間間隔乘上閥門個數與光速即為距離：. ΔR =. c × ΔGa 2. （3-1）. 其中 ΔGa 為一個閥門時間間隔。 ΔR 為一個閥門代表之距離，此計算得 0.46875 公尺。. 16.

(28) Ta Tb Tc. 圖 3-1 波形概要圖. 圖 3-2 脈衝接觸海水面之狀況示意圖. 17.

(29) 表 3-1 測高衛星資訊衛星名稱. 發射時間. 傾角 (度). 高度. 脈衝頻率. (km). 天線波寬 (度). 閘門時間. 儲存波形頻率. (ns). (Hz). SEASAT. 1978.6.27. 108. 800. 1.6. 30.5. 3.125. 10. GEOSAT. 1985.3.12. 108. 800. 60. 30.5. 3.125. 10. ERS-1. 1991.7.17. 98. 1020. 64. 32.5. 3.03(海洋模式) 12.12(冰面模式). 20. ERS-2. 1995.4.21. 1.3. 1020. 64. 32.5. 3.03(海洋模式) 12.12(冰面模式). 20. TOPEX(Ku). 1334. 1. 4500. 128. 32.5. 3.125. 10. 66. 1334. 2.7. 1200. 128. 35.5. 3.125. 5. 1992.8.10. 66. 1334. 1.1. 1700. 60. 29.5. 3.125. 20. GFO. 1998.2.10. 108. 800. 1.6. 1020. 128. 32.5. 3.125. 10. Jason-1(Ku). 2001.12.7. 66. 1334. 1.3. 1800. 104. 32.5. 3.125. 20. Jason-1(C). 2001.12.7. 66. 1334. 3.4. 300. 104. 32.5. 3.125. 20. 閘門個數. 預設閥門. 1020. 60. 2.1. 1020. 784. 1.3. 98. 784. 1992.8.10. 66. TOPEX(C). 1992.8.10. POSEIDON. (Hz). 18.

(30) 圖 3-3 Geosat/GM 海洋波形. 3-2 波形重定方法介紹. 由於有些脈衝不是發射在開放海域（open ocean）之中，所以產生之波形就不像海洋波形一樣有著清楚海洋波形形狀。會產生與海洋波形不同之波形，通常是因為脈衝接觸至海水面時同時接觸陸地或接觸到複雜海水面都有可能導致天線收到之脈衝有著複雜訊號。圖 3-4 為測高衛星地面軌跡接近澎湖產生之波形，由圖中可以很明顯看出，測高衛星預設閥門並不在前緣中點，而且有兩個斜面（ramp），所以吾人需要以波形重定重新定出前緣中點，如圖 3-4 共定出兩個前緣中點（以黑色虛線表示），經由重定後即可求得測高衛星觀測海水面之正確時間，再由（3-1）式換算成距離改正量，進而改正海水面高。下面則為波形重定方法介紹：. 19.

(31) 圖 3-4 Geosat 靠近陸地波形. 3-2-1 β參數演算法（β-parameter algorithm）. β參數演算法是由 Martin et al.（1983）發展，為第一個針對由冰面反射脈衝形成之波形進行波形重定演算法，其原理參考 Brown（1977）平均脈衝反射原理， β參數演算法需要有良好初始值，以進行最小二乘平差，再重複迭代後求得β參數，可分為β-5 參數（圖 3-5）與β-9 參數（圖 3-6）兩種，兩者非常類似，不過差別在於β-5 參數只可求解單一斜面波形，β-9 參數可以求解雙斜面波形，但是如果波形呈現尖錐（peak）形狀（如圖 3-7）時，使用β參數演算法將會使計算產生發散而無法得到計算成果，詳細β-5 參數計算方法可參照陳松安（2001），而β -9 參數計算方法可由β-5 參數做延伸。下面介紹β-5 參數之作法：. 20.

(32) ⎛ t − β3 ⎞ ⎟⎟ y (t ) = β1 + β 2 (1 + β 5 Q )P⎜⎜ ⎝ β4 ⎠. 0 ⎧ Q=⎨ ⎩t − (β 3 + 0.5β 4 ). P( x ) = ∫. x. −∞. （3-2）. for t < β 3 + 0.5β 4 for t ≥ β 3 + 0.5β 4. （3-3）. ⎛ − q2 ⎞ ⎟⎟dq exp⎜⎜ 2π ⎝ 2 ⎠ 1. （3-4）. 而β-9 參數為：. ⎛ t − β7 ⎛ t − β3 ⎞ ⎟⎟ + β 6 (1 + β 9 Q2 )P⎜⎜ y (t ) = β1 + β 2 (1 + β 5 Q1 )P⎜⎜ ⎝ β4 ⎠ ⎝ β8. 0 ⎧ Q1 = ⎨ ⎩t − (β 3 + 0.5β 4 ). for t < β 3 + 0.5β 4. 0 ⎧ Q2 = ⎨ ⎩t − (β 7 + 0.5β 8 ). for t < β 7 + 0.5β 8. ⎞ ⎟⎟ ⎠. for t ≥ β 3 + 0.5β 4. for t ≥ β 7 + 0.5β 8. 其中 y(t ) ：第 t 個閥門值，以 Geosat 而言，t 由 1 至 60。. β1 ：熱雜訊。 β 2 及 β 6 ：回傳訊號振幅。. β3 及 β 7 ：前緣中點，可與預設閥門差值計算為距離改正量。 β 4 及 β8 ：上升時間，由接收到能量至最大振幅一半。. 21. （3-5）. （3-6）. （3-7）.

(33) β5 及 β9 ：下降區斜率。. 圖 3-5 β-5 參數示意圖. 22.

(34) 圖 3-6 β-9 參數示意圖. 圖 3-7 尖錐波形示意圖. 23.

(35) 3-2-2 OCOG 演算法. Wingham et al.（1986）發表 OCOG 演算法，主要概念可參考圖 3-8，是以數值方式統計出波形振幅（amplitude）、寬度（width）與重心位置（center of gravity, COG），而 Deng （2004）修改公式，主要目的是為了降低熱雜訊影響，下面為 OCOG 公式：. A=. 60 − n. 60 − n. ∑ Pi 4 (t ). ∑ P (t ) 2. i =1+ n. i =1+ n. ⎛ 60−n ⎞ W = ⎜ ∑ Pi 2 (t )⎟ ⎝ i =1+ n ⎠. COG =. 2. i =1+ n. LEG = COG −. （3-8）. 60 − n. ∑ P (t ) 4. i =1+ n. 60 − n. ∑ iPi 2 (t ). i. i. （3-9）. 60 − n. ∑ P (t ) 2. i =1+ n. i. （3-10）. W 2. （3-11）. 其中 Pi (t ) 為閥門值。 A 為振福。. W 為寬度。 COG 為波形之重心。 LEG 為前緣中點。以 OCOG 計算得到之前緣中點，可以直接使用在改正海水面高，也可以當成β參數初始值，但是使用 OCOG 得到之海水面高起伏較震盪，精度也較低，所以大都用來求得初始值為主。. 24.

(36) 圖 3-8 OCOG 示意圖. 3-2-3 門檻值演算法. 門檻值演算法（Davis, 1997）是以 OCOG 為基礎做計算，為純統計演算法，利用不同門檻值乘上振幅得到一個閥門值，利用該值與前緣內插可以求得閥門： 5. PN =. ∑P i =1. i. 5. （3-12）. Tl = ( A − PN ) ⋅ Th + PN. Gr = Gk − 1 +. （3-13）. Tl − Pk −1 Pk − Pk −1. （3-14）. 25.

(37) 其中 Pi 與 Pk 都為該閥門之閥門值。 PN 為波形前 5 個閥門之閥門值平均值。. Th 為門檻值，本研究因考慮前緣中點為振幅ㄧ半之原則而使用 50﹪。 Gk 為第 k 個閥門值大於 Tl 之閥門。. Gr 為前緣中點。如遇到 Pk 與 Pk −1 之閥門值相同時，則 k 以 k+1 代入。門檻值演算法比較適合單一斜面之波形，如果為雙斜面波形，其計算而得之前緣中點通常以第一個閥門值大於 Tl 為主，非常之不客觀，且雙斜面如果以門檻值計算，通常無法得到該斜面預設門檻值。. 3-2-4 改良型門檻值演算法（improved threshold algorithm）. 上一節提到門檻值演算法，由於某些複雜波形無法計算出較準確前緣中點，所以在此提供了改良型門檻值演算法。圖 3-9 為改良型門檻值演算法流程圖，主要概念是先由相隔一個閥門之兩個閥門值比較，如果相減平均後大於 ε 1 則表示衛星開始收到能量，接著由相鄰閥門值相減，如果大於 ε 2 則表示衛星持續接收到能量，則將該閥門算入次波形（sub-waveform），並繼續檢查下一個閥門，直到小於 ε 2 則停止檢查，接著將剛剛由檢查 ε 2 得到之閥門範圍前後各加 4 個閥門，則形成次波形，只取 4 個閥門是為了減少過多閥門而產生之雜訊。將此次波形依序帶入 OCOG 演算法與 50﹪門檻值演算法即可得到前緣中點，由於此方法可以得到不只一個前緣中點，所以必須於最後加入一項判斷，本研究是以上一筆觀測資料之海水面高為基準，尋求最接近之海水面高，即可確定哪個前緣中點為吾人所需之時間點，如果遇到上一筆資料有問題者，則以上上筆資料為基準，依此類推。此外，不以先驗大地水準面高度為準則，是由於所觀測之海水面高與大地水準面高度有時會. 26.

(38) 存在一差值，如遇上波形重定後之海水面高反而更靠近大地水準面時，會使觀測量無法落在一連續面上，導致資料產生更大之誤差。. 3-3 選擇波形重定方法. 選擇適合波形重定方法非常重要，吾人針對 Geosat/GM 測高衛星波形以下列兩點原則來選擇最佳波形重定演算法：（1）解算波形成功率：能以該波形重定方法成功處理之波形數目越多，則可以說該波形重定方法越好。（2）沿軌跡波形重定後海水面高與大地水準面之比較：波形重定可得到距離改正量而得到改正後海水面高，如果改正後海水面高較平滑且與大地水準面較吻合者，則可認為該波形重定方法較好。. 吾人以 1985 年第 206 日（以 85206 表示）之下降軌道（如圖 3-10）為例來探討上述兩點原則，以此選擇本研究所使用之波形重定方法，由於 OCOG 演算法於計算後發現本身精度不高，所以不列入下列分析之中。以前述（1）之選擇來說，參考表 3-2，說明 85206 之衛星地面軌跡所有波形個數、波形重定可處理個數以及波形重定解算成功率，可以看出β-5 參數演算法大約只能處理 70﹪之波形，門檻值演算法可以處理波形之比例達到 100﹪，而改良型門檻值演算法能處理波形比例也接近 100﹪。接著分析（2），參考圖 3-11，為三個演算法計算得到之海水面高與大地水準面關係圖，大地水準面以黑色線段表示，而以線段與三角形表示者為β-5 參數演算法，以線段與正方形表示者為門檻值演算法，以線段與菱形表示者為改良型門檻值演算法。由圖中可以看出，由於β-5 參數演算法解算成功之比例較小，且大部分不能解算之波形都較靠近陸地，雖然於深海中有些海水面高較接近大地. 27.

(39) IF i = 50. 結束. i = i +1. 由第i個閥門開始，i預設為5 If. Pi + 2 − Pi > ε1 2. NO. YES Flag1 = i. k = k +1. 相鄰兩個閥門值相減 If Pi +1 − Pi > ε 2 i=Flag1+Flag2 k=0. NO Flag2 = k. 由i–4至i+k+4 組成次波形. 將次波形帶入OCOG與 Threshold得到前緣中點. 圖 3-9 改良型門檻值演算法流程圖. 28. YES.

(40) 圖 3-10 Geosat/GM 衛星 85206 地面軌跡圖. 29.

(41) 水準面，可是不符合本研究之需求，所以不考慮之。而門檻值演算法雖然能夠處理波形比例高達 100﹪，但是由其改正後海水面高可以發現較接近陸地處，其海水面高起伏較為劇烈，也不適合使用，而改良型門檻值演算法雖然可以處理波形之比例較門檻值演算法為低，可是由改良型門檻值演算法改正後之海水面高可以發現起伏不似門檻值演算法為大，最後可以數值方式呈現演算法成果好壞，參考表 3-3，為改正後海水面高與大地水準面之精度比較，由於β-5 參數演算法能解算之比例較小，所以用β-5 參數演算法能解算海水面高之數據為比較依據，可發現以改良型門檻值演算法為最佳，所以本研究決定採用改良型門檻值演算法。. 表 3-2 演算法解算波形成功率統計表全部個數. 可處理個數. 波形可處理比例. β-5. 285. 201. 70.53%. Threshold. 285. 285. 100.00%. Improved Threshold. 285. 283. 99.30%. 28 26. 海水面高(m). 24 22 20 Improve Threshold 18. Threshold Beta-5. 16 14 23.8. geoid 24. 24.2. 24.4. 24.6. 24.8. 25. 25.2. 緯度(degrees). 圖 3-11 改正後海水面高與大地水準面關係圖 30. 25.4. 25.6.

(42) 表 3-3 改正後海水面高與大地水準面精度統計表，單位：公尺標準偏差 Beta-5. 0.48. Threshold. 0.46. Improved Threshold. 0.26. 3-4 波形重定前後成果比較. 在取得資料比較前，吾人先將 10Hz 資料擬合成 2Hz 資料，擬合成 2Hz 是為了不讓 10Hz 資料對於後續應用之計算量產生太大負擔，而且 10Hz 對於重力計算來說，其雜訊會太多。擬合成 2Hz 還有一項好處，沿著軌跡（along-track）解析度大約為 7 公里，而相鄰軌跡解析度於赤道空間分布上大約為 4 公里，如果取 2Hz 可以使沿著軌跡與相鄰軌跡之精度相當。10Hz 資料於地面軌跡大約為 7 公里，而距離 7 公里可能橫跨海與陸地，由於衛星測高得到之海水面高，在進入陸地後，其高程會突然上升，如果將 10Hz 資料以二次曲線擬合成 2Hz 有可能會造成擬合曲線曲度過小，導致 2Hz 資料都不能用，所以吾人以 5 筆資料平均為 1 筆，並以三倍中誤差為準則去除異常值，其目的在於提高靠近陸地可使用資料量。. 3-4-1 單一軌跡比較. 吾人取 1985 年第 282 日（以 85282 表示）上升軌道（如圖 3-12），衛星地面軌跡由台灣進入台灣海峽，圖中軌跡為 2Hz，依圖上箭頭方向按照順序共取 40 個波形（圖 3-13 與圖 3-14）作為展示，圖 3-13 為第 1 至第 20 個波形，第 1 個波形為衛星地面軌跡剛由陸地進入海洋，圖 3-14 為第 21 個至第 40 個波形，第一個波形離陸地最近而第 40 個波形離海岸最遠，由第 1 至第 3 個波形可發現脈衝能量因. 31.

(43) 圖 3-12 Geosat/GM 衛星 85282 地面軌跡圖. 32.

(44) 接觸地面而發散，以至於衛星沒有接收到足夠脈衝能量，第 4 至第 12 波形由於陸地或者是複雜海水面影響，波形有緩慢上升趨勢，而第 13 個坡形之後即可看出有海洋波形形狀出現。參考圖 3-15，由左至右為台灣至台灣海峽方向，顯示波形重定前後之海水面高比較，黑色線為台灣最新大地水準面（呂誌強，2004），藍色線段為 GDRs 提供之原始資料，紫色線段為原始資料經由波形重定後得到之海水面高，三者皆為 2Hz 海水面高資料，可以看出原始資料較接近陸地處，海水面高開始向上跳躍，而經由波形重定後可以矯正海水面高並且使海水面高較接近大地水準面。. 圖 3-13 由 85282 軌跡取波形 1-20，左上角之數字表示圖 3-9 箭頭順序. 33.

(45) 圖 3-14 由 85282 軌跡取波形 21-40，左上角之數字表示圖 3-9 箭頭順序. 24 retracked 22. raw. 海水面高(m). geoid 20 18 16 14 12 23.2. 23.4. 23.6. 23.8. 24. 24.2. 緯度(degrees). 圖 3-15 85282 波形重定前後沿軌跡之海水面高比較. 34. 24.4. 24.6.

(46) 3-4-2 全區比較. 吾人以東經 119.5 至 122.5 度、北緯 21.5 至 25.5 度為本實驗範圍，包含了台灣本島、澎湖、綠島以及蘭嶼，範圍內所有衛星地面軌跡資料可以參考圖 3-16，共使用 165 條地面軌跡與 9055 個測高觀測數據，本實驗使用改良型門檻值演算法，並分別針對全區、離海岸線 20 公里、離海岸線 10 公里資料分析波形重定前後與大地水準面之比較，最後還針對區域以不同距離作分析。. 圖 3-17 為原始海水面高與大地水準面差值圖，圖 3-18 為波形重定後之海水面高與大地水準面之差值圖，由這兩張圖可以看出，波形重定不只對於近海岸附近之海水面高有一定程度改善，對於深海中由於測高衛星本身觀測資料不甚良好者，也可以進行改善。. 接著以統計方式來顯示結果，本研究以下列公式計算改善率：. IMP =. STD Retracked − STD Raw × 100% STD Raw. （3-15）. IMP 為改善百分比。. STD Retracked 為波形重定後海水面高與大地水準面高度之標準偏差。 STD Raw 為波形重定前海水面高與大地水準面高度之標準偏差。. 表 3-4 顯示波形重定前後，離陸地不同距離之海水面高與大地水準面比較，而圖 3-19 以直方圖形式取代表 3-3，由表可以發現，全區比較之精度改善 7.8﹪，離陸地 20 公里以內精度改善 16.8﹪，而距離陸地 10 公里內精度更可以改善 17.71﹪，表示波形重定對於測高資料距離陸地越近者，其改善程度越大。. 35.

(47) 圖 3-16 實驗區域與衛星地面軌跡圖. 36.

(48) 圖 3-17 原始海水面高與大地水準面差值圖，單位：公尺. 37.

(49) 圖 3-18 波形重定後海水面高與大地水準面差值圖，單位：公尺. 38.

(50) 表 3-4 不同距離對波形重定前後精度比較，單位：公尺 STD 全區. 重定前. 重定後. 改善率. 0.744. 0.686. 7.80%. 沿海 20 公里. 1.274. 1.06. 16.80%. 沿海 10 公里. 1.92. 1.58. 17.71%. 2. STD(m). 1.5. raw. 1. retracking. 0.5. 0 全區. 沿海20公里. 沿海10公里. 區域. 圖 3-19 不同距離對波形重定前後精度比較直方圖. 最後再以分區方式統計成果，將全部區域分為右上、左上、右下及左下。表 3-5 為全區統計表，可以看出右上方經由波形重定後改善最多，左上方呈現結果較差，有可能為該區資料過於稀疏，或是資料本身就不好之影響，表 3-6 顯示測高資料距離陸地 20 公里內之精度比較，一樣是右上方改善程度最大，而改善程度最小為左下方區域，表 3-7 表示資料距離陸地 10 公里內之精度比較，改善程度與 20 公里以內相同，估計原因可能為澎湖島附近資料經波形重定後也無法將其海水面高校正至正確海水面所引起。. 39.

(51) 表 3-5 全部區域分成不同區域對波形重定前後精度比較，單位：公尺全區. 重定前. 重定後. 改善率. 東北方. 0.603. 0.477. 20.90%. 西北方. 0.534. 0.553. -3.56%. 東南方. 0.56. 0.515. 8.04%. 西南方. 0.989. 0.873. 11.73%. 表 3-6 沿岸 20 公里分成不同區域對波形重定前後精度比較，單位：公尺 20 公里以內. 重定前. 重定後. 改善率. 東北方. 1.284. 0.95. 26.01%. 西北方. 0.624. 0.5. 19.87%. 東南方. 0.856. 0.72. 15.89%. 西南方. 2.112. 1.813. 14.16%. 表 3-7 沿岸 10 公里分成不同區域對波形重定前後精度比較，單位：公尺 10 公里以內. 重定前. 重定後. 改善率. 東北方. 2.009. 1.476. 26.53%. 西北方. 0.828. 0.625. 24.52%. 東南方. 1.358. 1.119. 17.60%. 西南方. 2.819. 2.435. 13.62%. 40.

(52) 第四章重力異常計算. 4-1 計算重力異常方法. 本文呈現三種計算重力異常結果，第一種為 IVM 計算得重力異常，其流程參考圖 4-1，為先移去長波長大地水準面梯度，再將得到之殘餘梯度組成南北、東西垂線偏差分量，隨之以一維快速傅立葉轉換進行計算並考慮最內圈影響（innermost zone effects）而得到殘餘重力異常，最後將殘餘重力異常、最內圈重力異常與長波長重力異常相加即得最後重力異常，該方法可詳見 Hwang（1998）。第二種為最小二乘配置法搭配去除-回復法（remove-restore method），先將測高觀測數據計算之大地水準梯度去除長波長大地水準梯度而求得殘餘大地水準梯度，隨之以最小二乘配置法將殘餘大地梯度轉換為殘餘重力異常，最後再加上長波長重力異常後，即可得到重力異常，該方法可詳見 Hwang and Parsons（1995）。第三種與第二種非常類似，不過運用了最小二乘配置法可以同時計算不同資料之原理，同時處理衛星測高觀測量與空載重力異常進而得到重力異常，第二種與第三種之流程圖可參考圖 4-2。. 本文使用之長波長重力位模式為 GGM02C（詳見 GRACE Home Page）加上 EGM96（詳見 CDDIS Home Page）模式。其中 GGM02C 球諧係數是由第 2 階至 200 階，EGM96 球諧係數為 201 階至 360 階。欲使用上述方法，則需搭配高斯濾波（Gaussian filter），下面介紹這些方法：. 41.

(53) 測高衛星觀測量 SSH. 全球大地位模式 GGM02C+EGM96. 計算大地水準梯度 e. elong eres = e − elong. Δg long. 組成南北、東西網格. IVM. Δg res. add. 重力異常. 圖 4-1 以 IVM 計算重力異常流程圖. 42. 計算最內圈效應.

(54) 測高衛星觀測量 SSH. 全球大地位模式 GGM02C+EGM96. 空載重力異常 Δga. Δgalong. 計算大地水準梯度 e elong. Δga res = Δga − Δgalong. Δg long. eres = e − elong. LSC. Δg = Δg long + Δg res. Δg res. 圖 4-2 以最小二乘配置法計算重力異常流程圖. 43. LSC.

(55) 4-1-1 高斯濾波. 由於 Geosat/GM 測高衛星觀測量經波形重定後較為震盪，所以先以高斯濾波平滑觀測量，下列為高斯濾波公式：. ⎡ 1 ⎛ ds ⎞ 2 ⎤ w(ϕ , λ ) = exp ⎢− ⎜ ⎟ ⎥ ⎣⎢ 2 ⎝ σ ⎠ ⎦⎥. 其中 Rs =. ds ≤ Rs. （4-1）. D 為一半搜尋視窗。 2. D 為搜尋視窗大小。. 而σ =. D 。 6. 搜尋視窗大小 D 影響著資料處理過後結果，如果搜尋視窗太大會造成資料過於平滑，失去真實觀測值，如果搜尋視窗太小則對原始觀測量沒有濾波作用也無法偵測到粗差。經由實作經驗，本研究以搜尋視窗大小 24 公里效果為最佳，參考圖 4-3，菱形為測高觀測量與大地水準面之差值，線段與三角形為測高觀測量與大地水準面之差值經由高斯濾波處理過後之結果，可以看出高斯濾波平滑測高觀測量之結果。. 44.

(56) 1.4 高斯濾波前. 與大地水準面差值 (m). 1.2. 高斯濾波後. 1 0.8 0.6 0.4 0.2 0 21. 21.5. 22. 22.5. 23. 23.5. 24. 緯度(degrees). 圖 4-3 高斯濾波前後示意圖. 4-1-2 去除-回復法. 去除-回復法為整合積分法、最小二乘配置法、球諧展開式之方法，可用於計算大地起伏。依照 Hwang（1997）可以將大地水準梯度分為三部份，長波長大地水準梯度、短波長大地水準梯度以及殘餘大地水準梯度，長波長大地水準梯度可以由全球大地位模式求得，短波長大地水準梯度為地形效應，殘餘大地水準梯度可由 4-2 式計算而得：. e res = e − elong − e RTM. （4-2）. 其中 eres 為殘餘大地水準梯度。. 45.

(57) e 為由觀測量計算而得大地水準梯度。 elong 為長波長大地水準梯度。. eRTM 為短波長大地水準梯度。利用衛星軌跡經緯度換算得觀測量之間距離而求得大地水準梯度，接著扣掉由長波長大地起伏（如圖 4-4）所計算得長波長大地水準梯度，再由公式 4-2 計算得殘餘大地水準梯度，由於海水面之短波長大地水準梯度非常微小且趨近於零，所以忽略不加以計算。. 重力異常也可以分為三部份，分別為長波長重力異常、短波長重力異常及殘餘重力異常。長波長重力異常由全球擾動位模式計算而得，短波長重力異常為地形效應所引起，殘餘重力異常可由殘餘大地水準梯度以最小二乘配置法求得，所以計算重力異常可用下列公式：. Δg = Δg long + Δg RTM + Δg res. （4-3）. 其中 Δg 為重力異常。 Δg long 為長波長重力異常。. Δg RTM 為短波長重力異常。 Δg res 為殘餘重力異常。求得殘餘重力異常後即可帶入公式 4-3 以得到重力異常 Δg，而長波長重力異常（如圖 4-5）之計算公式為：. Δg long. GM = 2 a. N. n. n=2. m =0. ∑ (n − 1)∑ (C. nm. cos mλ + S nm sin mλ )Pnm (sin φ ). （4-4）. 而短波長重力異常由於資料位於海水面而趨近於零，所以忽略不予以計算。. 46.

(58) 4-1-3 IVM. IVM 是將垂線偏差以最小二乘配置法組成南北、東西分量，下一小節會介紹最小二乘配置法，隨之以一維快速傅立葉轉換計算為殘餘重力異常，最後再加上參考重力場與最內圈效應影響後即可得重力異常，圖 4-2 為流程圖。詳細方法可參考 Hwang（1998），其計算公式為：. Δg ( p ) = −. γ0 4π. ⎛ ∂ψ pq. ∫∫ σ H ′⎜⎜. ⎝ ∂φ q. ∂ψ pq ⎞ ⎟ ⋅ (ξ qη q )dσ q cos φ∂λ q ⎟⎠. γ0 σ H ′(ξ q cos α qp + η q sin α qp )dσ q 4π ∫∫ γ = 0 ∫∫ σ H ′ε qp dσ q 4π =. ψ ⎛ ⎜ sin 3 pq 1 2 H (ψ pq ) = + log⎜ ⎜ ψ pq ψ pq sin ⎜ 1 + sin 2 2 ⎝. ⎞ ⎟ ⎟ ⎟ ⎟ ⎠. （4-5）. （4-6）. 其中 ε qp 為單位球上 q 點沿方位角 α qp 垂線偏差之分量。. ξ q 與η q 為 q 點上垂線偏差南北、東西分量。. γ0 =. GM 為平均重力。 R. φ 為緯度。 H (ψ pq ) 為一核函數。. 為求計算快速，本研究使用一維快速傅立葉轉換，該方法可將空間或時間轉換至頻率相乘計算，以加快計算速度。其公式為：. 47.

(59) Δg φ p = γ 0. ΔφΔλ φn λn ∑ ∑ H ′(Δλqp )(ξ cos φ cos α qp + η cos φ sin α qp ) 4π φq =φ1λq =λ1. ⎫⎪ ΔφΔλ −1 ⎧⎪ φn F1 ⎨ ∑ F1 (H ′(Δλ qp )cos α qp )F1 (ξ cos φ ) + F1 (H ′(Δλ qp )sin α qp )F1 (η cos φ ) ⎬ = γ0 4π ⎪⎩φq =φ 1 ⎪⎭. [. ]. （4-7）. 其中 Δλ qp = λ q − λ p 。 Δφ 與 Δλ 為緯度與經度方向上之網格間隔。. F1 為 1D FFT。. 4-1-4 最小二乘配置法. 最小二乘配置法原理為根據觀測數據以最小二乘方式擬合而求得一近似函數，進而推估未觀測資料之值（胡明城等，1993）。也是一種可以將不同觀測量結合在一起計算重力異常之方法（Moritz, 1980），可是最小二乘配置法需要根據先驗精度或是觀測值與估計值之間關係來計算協變方矩陣以估算未知值。其公式為：. Δg res = C Δge (C ee + C nn ) eres −1. （4-8）. 其中 Δg res 為殘餘重力異常。 C Δge 為殘餘大地梯度與殘餘重力異常協變方矩陣，是由 GGM02C 與 EGM96 模型. 計算得到，由交大提供計算好之協變方表格。 C ee 為殘餘大地梯度之協變方矩陣，由交大提供計算好之協變方表格。 C nn 為雜訊值，為一對角線矩陣，雜訊值越大則計算過後之結果越平滑（黃金維， 2003）。C nn 本應由衛星測高觀測量之標準偏差帶入，但由於測高數據本身較震盪，. 48.

(60) 所以全部以實作經驗所得 2.0 弧秒（arc-second）取代。. 最小二乘配置法可以結合不同觀測量，所以本文還使用空載資料與衛星測高資料一併進行計算，其公式為：. Δg res = (C ΔgΔga. ⎛ C Δga + C mm C Δge )⎜⎜ ⎝ C eΔga. ⎞ ⎟ C ee + C nn ⎟⎠ C Δgae. −1. ⎛ Δga r ⎞ ⎜⎜ ⎟⎟ ⎝ eres ⎠. （4-9）. 其中 Δgar 為空載殘餘重力異常，直接扣掉高度平均 5156 公尺之長波長重力異常而得。 C Δgae 、 C eΔga 為空載殘餘重力異常與殘餘大地水準面梯度之協變方矩陣，由交大提. 供計算好之協變方表格。 C mm 為雜訊值，為一對角線矩陣，其值以本身精度 3.0mgal 帶入。 C ΔgΔga 為殘餘空載重力異常與殘餘重力異常協變方矩陣，是由 GGM02C 與 EGM96. 模型計算得到，由交大提供計算好之協辯方表格。由於最小二乘配置法理論上需要對於待求解之重力場訊號為零，所以需要配合去除-回復法使用。對於所有重力場觀測資料，組成 C Δge 與 C ee 之前需要先移去長波長分量與短波長分量，再進行最小二乘配置運算而得到殘餘重力異常，隨之再對殘餘重力異常進行回復動作，即可得到重力異常。. 49.

(61) 圖 4-4 長波長大地起伏等高線圖，等高線間隔 0.5 公尺. 50.

(62) 圖 4-5 長波長重力異常等值線圖，等值線間隔 10mgal. 51.

(63) 4-2 成果分析. 利用 4-1 節方法計算之重力異常，以最小二乘配置法為例，圖 4-6 為原始測高資料計算之重力異常等值線圖，圖 4-7 為資料經波形重定後計算之重力異常等值線圖，由圖中對照可發現圖右上角之牛眼（Bull eye）消失不見，台北縣以北處、宜蘭外海龜山島附近及宜蘭沿海之重力異常密集也較為降低，澎湖重力異常密集處也較為緩和，嘉義以南至台南以北沿海部分與恆春南邊沿海其重力異常也有明顯改善，大部分改善都分布於沿岸，表示測高衛星接近陸地時，其量測能力大幅降低，可由波形重定將測得原始海水面高修正為正確海水面高。不過還是有些部份無法進行改正，如屏東西部沿海、澎湖與綠島部份還是存在著牛眼之情形，顯示經由波形重定無法修正該處測高觀測量，而宜蘭花蓮外海也存在著較為複雜之等值線分布，上述這些部分在在都顯示著波形重定並無法有效改善不好之觀測量。. 接著以二種方法計算得三種結果，分別為 IVM、最小二乘配置法、最小二乘配置法同時計算衛星測高資料與空載資料（以最小二乘配置法加空載替代），每個結果各別可繪出波形重定前後重力異常分布圖，圖 4-8 與圖 4-9 為使用 IVM 計算得到之重力異常分布圖，圖 4-10 與圖 4-11 為使用最小二乘配置法計算得到之重力異常分布圖，圖 4-12 與圖 4-13 為最小二乘配置法加空載計算得到之重力異常分布圖，可發現以 IVM 計算得到之重力異常，其最大最小值之絕對值較其餘兩種為小且全區之重力異常可看出較為平滑，不過於東岸重力異常值變化較大處，其顯示結果較差，而以最小二乘配置法加空載計算得到之重力異常與單純只用最小二乘配置法得到之重力異常來比較，大致上於重力異常之絕對值較小處比較沒有影響，而重力異常之絕對值較大處則有些微不同。. 接著與中央大學船測資料做比較，此處一樣是每個方法繪出波形重定前後之. 52.

(64) 成果，圖 4-14 與圖 4-15 為使用 IVM 計算得到之重力異常與船測重力異常差值圖，圖 4-16 與圖 4-17 為使用最小二乘配置法計算得到之重力異常與船測重力異常差值圖，圖 4-18 與圖 4-19 為最小二乘配置法加空載計算得到之重力異常與船測重力異常差值圖，單就以波形重定前後差值圖來比較，可以發現經由波形重定後，其改正處皆與上面提到之等值線圖改正部份大同小異，而對於不同方法上來討論，則發現以 IVM 計算得到之重力異常值與船測重力異常差異最大，而最小二乘配置法次之，最小二乘配置法加空載差異為三者中最小。. 最後以數值方式表達三種方法與船測重力比較結果，表 4-1 為圖 4-14 與圖 4-15 之值統計結果，表 4-2 為圖 4-16 與圖 4-17 之值統計結果，表 4-3 為圖 4-18 與圖 4-19 之值統計結果，表 4-1 之最大最小值之絕對值為三者最小，其原因為 IVM 計算得到之重力異常較為平滑，以改善率來說，IVM 只有提升 3﹪，而最小二乘配置法提升了 11﹪，最小二乘配置法加空載提升約 12﹪，所以最小二乘配置法加空載計算得到之重力異常與船測重力異常比較之精度為最佳。不過有一點值得注意，表 4-1、表 4-2 與表 4-3 之最大值與最小值，顯示著波形重定後所計算之重力異常與船測差值超過 100mgal，吾人將差值大於 50mgal 者挑出並分別繪製圖 4-20、圖 4-21 與圖 4-22，可發現與船測重力異常差值較大者全都在宜蘭外海處，其原因為測高資料經波形重定後無法改正該位置海水面高導致，一方面可以由此結果對於改正後海水面高作檢查，一方面也可以對波形重定後還是不好之資料作篩選。. 表 4-1 以 IVM 計算波形重定前後重力異常與船測重力異常之精度分析表，單位： mgal 最大值. 最小值. 標準偏差. 平均值. 均方根. 重定前. 118.54. -110.57. 20.06. -0.77. 20.07. 重定後. 113.02. -107.87. 19.42. 0.49. 19.42. 53.

(65) 表 4-2 以最小二乘配置法計算波形重定前後重力異常與船測重力異常之精度分析表，單位：mgal 最大值. 最小值. 標準偏差. 平均值. 均方根. 重定前. 143.22. -114.26. 13.48. -1.56. 13.57. 重定後. 137.83. -112.329. 11.98. -0.91. 12.01. 表 4-3 以最小二乘配置法加空載計算波形重定前後重力異常與船測重力異常之精度分析表，單位：mgal 最大值. 最小值. 標準偏差. 平均值. 均方根. 重定前. 134.05. -111.40. 12.58. 1.62. 12.69. 重定後. 130.09. -109.32. 11.05. 1.59. 11.17. 54.

(66) 圖 4-6 以最小二乘配置法計算測高資料波形重定前之重力異常等值線圖，等值線間隔 10mgal. 55.