測高與重力衛星資料應用於重力與水文變化

 

國 立 交 通 大 學

土 木 工 程 學 系

博 士 論 文

測 高 與 重 力 衛 星 資 料

應 用 於 重 力 與 水 文 變 化

Gravity and hydrological changes from satellite

altimetry and gravimetry

研 究 生 : 高 豫 麒

指 導 教 授 : 黃 金 維

 

測 高 與 重 力 衛 星 資 料

應 用 於 重 力 與 水 文 變 化

Gravity and hydrological changes from satellite

altimetry and gravimetry

研究生:高豫麒

Student : Yu-Chi Kao

指導教授:黃金維

Advisor : Dr. Cheinway Hwang

國 立 交 通 大 學

土 木 工 程 學 系

博 士 論 文

A Dissertation

Submitted to Department of Civil Engineering

College of Engineering

National Chiao Tung University

in Partial Fulfillment of the Requirements

for the Degree

of Doctor of Philosophy

In

Civil Engineering

June 2010

  _I

測高與重力衛星資料應用於重力與水文變化

研究生:高豫麒 指導教授:黃金維

國立交通大學土木工程學系

摘 要

本文的主要研究課題為(1)使用衛星測高(Satellite Altimetry，SA)與熱容(steric) 兩種資料計算海水面高度及質量之時變量，並探討海水質量變化對地球質心、扁 率及極運動之影響量。(2)使用 SA 與 GRACE (Gravity Recovery and Climate Experiment)重力衛星資料研究監測湖泊水位高度變化之可行性。本文以 1992 年 10 月至 2007 年 12 月間的 SA 與 steric 資料，獲得以下研究結果：(1)海水面高度 異常(Sea Leavel Anamaly，SLA)、steric 與改正後海水面高度異常(Corrected SLA， CSLA)之全球海水面高度年變率(Global Sea Level Trend，GSLT)結果，分別為 2.94±0.18 mm/yr、0.32±0.28 mm/yr 及 2.47±0.14 mm/yr；(2)在多個海域中 CSLA 的海水面高度年變率與SLA 及 steric 的結果相反，顯示此些海域受到 steric 的效 應影響極大；(3)由 CSLA 反演計算的ΔC₀₀係數年變率為 1.16±0.07×10-10，振幅 為 5.13×10-10，相位為- 81.2°；(4)Δx、Δ 與 zy Δ 地球質心偏移量年變率分別為 -0.105±0.015 mm/yr、0.011±0.019 mm/yr 及-0.234±0.015 mm/yr ；(5) J2年變率為

0.57±0.08×10-11/yr，振幅為 1.98×10-11/yr，相位為 -127.21°；(6) ΔC₂₁的變化量約 為ΔS₂₁係數的3 倍，即極移運動中 x 方向變動量大於 y 方向約 3 倍。

本文使用SA 測算的湖水位、每月的 GRACE 重力觀測量與美國氣候預測中

心(Climate Prediction Center，CPC)水文模式的等效水位高度(Equivalent Water Height，EWH)變化量等資料，推估 Baikal 與 Balkhash 湖集水域內，流入這兩個 湖泊的水量，分別為 60% 和 30%。在相位變化方面， Baikal 湖區 GRACE 結 果與 SA 及 CPC 模式相差約 7 個月，此結果差異主要影響因素為湖泊座落位置 的氣候環境與人為干擾量。本文中比較 Baikal 與 Balkhash 湖各兩個不同時間段 的月均溫，其結果均顯示後時間段的月均溫均較前時間高，分別高約0.64 °C 與 0.7 °C ，其中 Baikal 湖的年均溫由零度以下，上升至零度以上，此結果將造成 原永凍土的解融，增加注入 Baikal 湖的水量。對照兩湖區年均溫及地球扁率 J2 係數的變化趨勢，兩者均在1997-1998 El Niño 事件發生後變化轉折。

  _II

Gravity and hydrological changes from satellite

altimetry and satellite gravimetry

Student : Yu-Chi Kao Advisor : Dr. Cheinway Hwang

Institute of Civil Engineering

National Chiao Tung University

Abstract

The main research topics of this study are : (1) to determine changes of sea surface heights and oceanic mass using satellite altimetry (SA) and steric heights, and to determine the contributions of such changes to variations of geocenter, J2, and

polar motions. (2) to determine lake level changes using SA and Gravity Recovery and Climate Experiment (GRACE) gravimetry data. Using SA and steric heights over October 1992 to December 2007, the following results are obtained : (1) the rates of sea level anomaly (SLA), steric height anomaly and corrected SLA (CSLA) are 2.94 ± 0.18, 0.32 ± 0.28 mm/yr and 2.47 ± 0.14 mm/yr. (2) The trends of CSLA, SLA and steric heights are conflicting over some oceanic areas, suggesting that large uncertainties of the underlying data may exist here. (3) The rate, annual amplitude and phase of ΔC₀₀, which is the oceanic mass-induced variation of the zero-degree geopotential coefficient, are 1.16 ± 0.07 × 10-10, 5.1 3 × 10-10 and -81.2º. (4) The rates of the geocenter variations in x, y, and z are -0.105 ± 0.015, 0.011 ± 0.019 and -0.234 ± 0.015 mm/yr. (5) The rate and annual amplitude and phase of J2 are 0.57 ±

0.08 × 10-11/yr, 1.98 × 10-11/yr and -127.21 . (6) The∘ ΔC₂₁ magnitude of is 3 times greater than that of ΔS₂₁, suggesting that the oceanic mass contributes more to the x component of the polar motion than the y component.

  _III

Using the variations of lake level from SA and the equivalent water heights (EWH) from the monthly GRACE gravity fields and the Climate Prediction Center (CPC) hydrological model, it is estimated that about 60% and 30% of water in the catchment areas of Lakes Baikal and Balkhash enter these two lakes. The phase of the annual variation of the Baikal lake level from GRACE differ by 7 months from the phases obtained with satellite altimetry and the CPC model. This difference is due to a climate factor and some man-made causes. Monthly temperature data over two time periods around Lakes Baikal and Balkhash show that the mean temperature in the earlier period is larger than that in the later period. On average, temperatures over Baikal and Balkhash increase by 0.64° and 0.7°C. Over these two periods the mean temperature around Baikal changes from sub-zero to above-zero, suggesting that the permafrost here may thaw to increase the amount of water entering Baikal. A common feature in the trends of J2 and the lake levels of Baikal and Balkhash is that

  _IV

誌 謝

畢業了！感觸真的很多！心中想感謝的人真的是太多了，首先必須感謝我的 指導教授黃金維老師，因為他的鼓勵與協助，才能讓我順利畢業拿到學位。其次 是咱們溫馨快樂的研究團隊，我退伍後這兩年半與各位學弟妹們，一起討論問 題，一起揪團吃飯，一起打屁聊天，絕對是我人生歷程中一段值得回憶的生活歷 練，在此必須將咱們團隊人員的姓名列出，以免日後忘記，畢竟自知對人名的記 憶力太差了！鄭景中博士林廷融、、鍾易達、曾子榜、施亘昶、高瑞其、黃鉅富、 連紫漪、黃啟訓、魏祥鴻、Natt(陳英雄) 、李宜珊、許宏銳、姚鐙凱、蘇洵頡、 陳逸如、彭千惠、孫佳龍、高綉雯、楊恩銘、陳彥杕。 最後要感謝的當然是我的老婆大人，畢竟多次的心情轉折，均是老婆大人 的撫慰與支持，才讓我重拾信心與鬥志。老婆大人，謝謝妳！

  _V

目 錄

中 文 摘 要… … … . . … … … … Ⅰ 英 文 摘 要… … … . . … … … Ⅱ 誌 謝… … … . . … … … . . … … … Ⅳ 目錄……… V 表目錄………...…Ⅷ  圖目錄……….. Ⅸ  第 1 章 緒 論 … … … . … 1 第 1 - 1 節 研 究 緣 起 … … … . … . … 1 第 1 - 2 節 研 究 方 向 及 流 程 … … … . . … … 3 第 1 - 3 節 各 章 節 概 述 … … … . . … 6 第 2 章 球 諧 函 數 係 數 的 解 析 與 整 合 … … … . . … 7 第 2 - 1 節 章 節 緒 論 … … … 7 第 2 - 2 節 完 全 正 規 諧 函 數 係 數 與 普 通 諧 函 數 係 數 之 關 係 式 . . … . . 7 第 2-3 節 水文高度及大氣壓力異常值重力位係數計算式推導……9 第 2 - 4 節 F F T 計 算 技 巧 於 球 諧 係 數 解 析 之 應 用 … … … … . . … … 11 第 2 - 5 節 F F T 計 算 技 巧 於 球 諧 係 數 整 合 之 應 用 … … … . . … 1 4 第 2 - 6 節 Wa v e l e t 轉 換 計 算 分 析 與 應 用 … … … . . 1 5 第 2 - 7 節 本 章 小 結 … … … . . … 1 9 第 3 章 資 料 來 源 與 處 理 … … … . … … 2 0 第 3 - 1 節 章 節 緒 論 … … … . … … … . … 2 0 第 3 - 2 節 AV I S O 衛 星 測 高 資 料 … … … . … … . . … … … 2 0 第 3 - 3 節 s t e r i c 熱 容 效 應 資 料 … … … . … … … . … … 2 4 第 3 - 4 節 G R A C E 資 料 … … … … . … … … . … … … … … … . … … 2 6 第 3 - 5 節 C P C - L D A S 水 文 模 式 資 料 … … … . … . … 2 9

  _VI 第 3 - 6 節 本 章 小 結 … … … . . … … … … . … . … 3 0 第 4 章 衛星測高資料成果分析………..……….………32 第 4 - 1 節 章 節 緒 論 … … … . … … 3 2 第 4 - 2 節 海 水 面 高 度 變 化 相 關 時 變 量 之 成 果 探 討 … … … . . … … 3 3 第 4 - 3 節 海 水 面 高 度 異 常 值 與 E l N i ñ o 之 關 係 … … … 4 4 第 4 - 4 節 海 水 質 量 時 變 量 對 低 階 重 力 係 數 之 影 響 … … … . … … 4 9 第 4-4-1 節 海水質量時變量:C0 0 係數………...………49 第 4-4-2 節 地球質心坐標變率：C1 0、C11、S11 係數…..……52 第 4-4-3 節 J2 係數的變化 …….……….………..……58 第 4-4-4 節 極運動：C2 1 及 S2 1 係數……….……62 第 4 - 5 節 本 章 小 節 … … … . … … … 6 6 第 5 章 測 高 及 重 力 衛 星 分 析 水 文 變 化 … … … . … … 6 9 第 5-1 節 章節緒論………..…………69 第 5-2 節 湖泊簡介………..…70 第 5-3 節 衛星測高資料應用於湖水位高度變化監測之研究分析…73 第 5-4 節 GRACE 測算湖水位高度之理論與方法………77 第 5-5 節 GRACE 湖水位高度變化結果分析探討………...85 第 5 - 6 節 J2 地 球 扁 率 變 化 與 湖 水 位 高 度 變 化 之 關 連 性 探 討… … . 8 9 第 5-7 節 本章小結……….………...………..91 第 6 章 結 論 與 建 議 … … … . … … … . 9 3 參 考 文 獻… … … . . … … … . . 9 8 附 錄 A : 球諧係數解析程式 sha.f90 簡介與應用………111 附 錄 B : 球諧係數整合程式 syn.f90 簡介與應用………...…115

  _VII

附 錄 D : 本 文 所 使 用 資 料 下 載 網 址 … … … . . … 1 2 1 附 錄 E : 英 文 縮 寫 對 照 … … … . … … … . . … 1 2 2 作 者 簡 歷… … … . . … 1 2 4 學 術 著 作 目 錄… … … . . … 1 2 5

  _VIII

表 目 錄

表3-1：1997 年各月份 AVISO DT-MSLA 1/4°網格及本文重新取樣為 1°×1°網格後之 SLA 標準偏差表………22 表3-2：CPC 及 GRACE 2007 年 6 月資料之 PFM 及高斯濾波平滑化 測試成果表………28 表4-1：1992 年 10 月至 2007 年 12 月間 SA、steric 與 CSLA 之 GSLT 成果………35 表4-2：1992 年 10 月至 2007 年 12 月間 CSLA 之全球海面區域重力異常 及大地起伏年變率成果………..………41 表4-3：ΔC₀₀與轉換後海水質量時變相關計算結果………50 表4-4：地球質心年振幅與年周期變化結果………54 表 4-5：不同研究使用 SLR 觀測量計算 J2 變率統計表………….………58 表 4-6：J2 時變計算成果...60 表 5-1：Baikal 與 Balkhash 湖的相關地理氣候環境參考資料………72 表 5-2：Baikal 與 Balkhash 湖於所選期間內的湖水位年變率………76

  _IX

圖 目 錄

圖3-1：太平洋亞熱帶地區 (N20°~40°，E100°~240°) 2007 年 1 月 18 日的 SLA 觀測量與其誤差量的百分比精度評估圖(摘自 AVISO)…………23 圖3-2：1997 年 11 月全球 SLA 圖， (a)及(b)的資料網格分別為 1/4°及 1°×1°…23 圖3-3：1997 年 11 月 steric 效應之海水面高度異常………..………25 圖3-4：高斯平均濾波值曲線圖，計算階數 60 階….…….………27 圖 3-5：CPC 模式與 GRACE 2007 年 6 月全球 EWH 異常圖………….………28 圖3-6：2002 年 10 月至 2007 年 12 月間 CPC 水文模式之 EWH 變化率………30 圖4-1：1992 年 10 月至 2007 年 12 月全球海水面高度年變率。(a)SLA、 (b)steric、(c)CSLA ………37 圖4-2：1992 年 10 月至 2007 年 12 月全球海水面年相位變化。(a)SLA、 (b)steric、(c)CSLA ………39 圖4-3：1992 年 10 月至 2007 年 12 月全球海水面年震幅變化。(a)SLA、 (b)steric、(c)CSLA ………40 圖4-4：全球主要洋流流向及位置圖(引摘自MiraCosta College) ………41 圖4-5：1992 年 10 月至 2007 年 12 月 CSLA 重力異常年變率圖…………..….42 圖4-6：1992 年 10 月至 2007 年 12 月 CSLA 重力異常年相位變化圖…………42 圖4-7：1992 年 10 月至 2007 年 12 月 CSLA 重力異常年振幅變化圖…………42 圖4-8：1992 年 10 月至 2007 年 12 月 CSLA 大地起伏年變率圖………...……43 圖4-9：1992 年 10 月至 2007 年 12 月 CSLA 大地起伏年相位變化圖…….……43 圖 4-10：1992 年 10 月至 2007 年 12 月 CSLA 大地起伏年振幅圖………43 圖4-11：太平洋 NINO_3.4 區之海水表面溫度異常指標圖………46 圖4-12：1992 至 2007 年間 El Niño 現象發生之年月海水面高度異常變化圖…47 圖4-13：1992 至 2007 年間 La Niña 現象發生之年月海水面高度異常變化圖…48 圖4-14：海水質量時變折線與 wavelet 小波頻譜分析圖………51 圖4-15：地球質心Δx時變與wavelet 小波頻譜分析圖……..………55

  _X 圖 4-16：地球質心Δ 時變與 wavelet 小波頻譜分析圖………56 y 圖 4-17：地球質心Δ 時變與 wavelet 小波頻譜分析圖………57 z 圖 4-18：J2 時變圖………60 圖 4-19：CSLA 計算得之 J2係數時變與 wavelet 小波頻譜分析圖………61 圖 4-20：ΔS₂₁係數變化與小波分析圖………64 圖 4-21：ΔC₂₁係數變化與小波分析圖……….………65 圖5-1：Baikal 及 Balkhash 湖泊的地形和河流系統………72

圖5-2：T/P 衛星通過 Baikal (a)及 Balkhash (b) 湖的地面軌跡………75

圖5-3：Baikal (a) 及 Balkhash (b) 湖水位高度變化曲線………76

圖5-4：2007 年 6 月西伯利亞地區 GRACE (a)與 CPC (b)資料計算之等水位 高度變化量比較圖………80

圖5-5：2002 年至 2009 年間 Baikal 及 Balkhash 湖區域重力異常年變率………81

圖 5-6：GRACE 資料所計算之年振幅 (a) 與年相位圖 (b)………82

圖 5-7：Baikal (a) 與 Balkhash (b) 湖之集水域範圍與網格分割圖………83

圖5-8：Baikal (a) 與 Balkhash (b)之 SA、GRACE 與 CPC 模式等 3 種資料 計算之湖水位高度變化曲線……….………84

圖5-9：Baikal (a) 與 Balkhash (b)湖月平均降水量折線圖………87

圖 5-10：Baikal (a) 與 Balkhash (b)湖不同時間段的月均溫比較圖………88

附錄圖1：1997 年 11 月 SLA 圖……….……118

附錄圖2：1992 年 10 月至 2007 年 12 月海水質量變化所造成的 C20 係數 wavelet 頻譜分析圖………...………120

  ₁

第 1 章 緒 論

第 1-1 節 研究緣起

重力場時變的問題及其機制早已為多位學者所提出研究探討（Yoder et al.， 1983 ； Rubancam，1984 ； Gegout and Cazenave，1993 ； Nerem et al .，1993 ； Eanes，1995 ； Dong et al.，1996 ； Cheng and Tapley，1999），其機制是因地 球內部及外部質量交換流動分佈變化所引起，例如地幔 (mantle) 內質量移動； 海水、大氣質量的變化；地潮 (earth tide) 及海潮 (ocean tide) 的變化；後冰河 期反彈 (Post Glacial Rebound，PGR) 及南極 (Antarctic)、格陵蘭 (Greenland) 冰 層融化的變化；地下水位、湖水位及降雪量變化等多種因素所造成。這些質量時 序變化與地球氣候環境的變遷是息息相關的，瞭解重力場時變的情形將有助於吾 人明瞭過去與現今的地球氣候變化情形，更甚而可用以評估預測未來環境的變 化，提供人類繼續生存於地球的應變之道。 瞭解及偵測現今GSLT 是研究全球氣候變遷的一項重要科學研究。在近 20 年來透過SA 方法與潮位站觀測資料，已經明確的證明全球海水面是呈上升趨勢

(Cazenave et al.，2004 ； Beckley et al.，2007)。目前應用於監測海水面高度變 化的觀測技術中，SA 結果是最直接的，而實際上也是現今偵測 GSLT 最有效的

方法。各研究機構及學者使用各種 SA 資料研究 GSLT 多有結果發表，例如

Leuliette et al. (2004) 與 Nerem et al. (2006) 使用 T/P 及 Jason-1 觀測資料精確計 算出1993 至 2003 年這 10 年間的 GSLT 為 2.8 ± 0.4 mm/yr； Kuo (2005) 則根據 1985 至 2002 年約 530 個驗潮站資料分析計算結果顯示 GSLT 分別為 2.9±0.5 mm/yr； Beckley et al. (2007) 使用重新架構在 ITRF2005 地形框架系統下的 T/P 及 Jason-1 觀測資料，所計算得 1993 至 2007 年間共 14 年之 GSLT 為 3.36 ± 0.41mm/yr。

  ₂

地球上各種生物的活動都與水有關，湖泊蓄積了大量的淡水，常常是一個

區域的生活中心，探討湖泊的水文循環 (Hydrologic cycle) 演化及利用情形，在

現今全球水資源逐漸匱乏的時刻，亦是一個值得重視的研究課題。有關GRACE

重力衛星觀測資料應用於重力時變與水儲存量變化之研究現象，例如Chen et al.

(2004) ； Davis et al. (2004) ； Wahr et al. (2004) ； Andersen and Hinderer (2005) 與 Swenson and Wahr (2009) 等。在研究一個湖泊集水域內的水文平衡變化，主 要為湖泊水位變化及集水域儲水量變化這兩個研究課題 (Bedient and Huber， 2002；Gordon and Famiglietti，2004)，充分的探討及瞭解這兩項水文變化課題， 將有助於水資源管理及應用。根據Wahr et al. ( 2004) 針對 GRACE 觀測資料的

研究，在數百平方公里及數星期的時空序列下，能有效的監測陸地EWH 變化，

其監測精渡在陸地可達1.5 cm rms，因此面積達數百平方公里的湖泊或水庫，理

論上應可利用GRACE 的每月重力觀測數據，經過計算推演後，監測這些湖泊的

水位高度。傳統監測湖水位 (Lake Level，LL) 高度變化，大多採用水位計 (stage

gage) 直接量測。然自SA 技術的發展與應用後，SA 已是一個用來監測湖泊水位

高度變化的有效方法 (Hwang et al.，2005；Brikett，1998)。SA 量測計算湖泊水

位變化的方法如同水位計一樣，屬直接量測。理論上SA 資料可用於陸地和海洋，

在海洋面上大多用以測算海洋環流、海潮、海洋重力及海水面高度變化等；而在

陸地上的兩個典型應用則為湖水位高度變化 (Lake level change，LLC) (Birkett， 1995 and 1998；Hwang et al.，2005) 與地面高程變化監測 (Berry et al.，1997)， 相關此類的應用與研究報告可參考Hwang et al. (2008)。 綜上所述，本文的研究目的有二：(1) 期透過 SA 技術測算全球海水面高度 變化量，並以海水質量變化量來研究其對地球質心、扁率及極運動的貢獻量；(2) 以GRACE 所觀測到的重力時變係數資料，經由第二章之理論計算方法，探討監 測LLC 之可行性，西伯利亞的貝加爾湖 (Baikal) 及中亞的巴爾喀什湖 (Balkhash) 是本研究所選定的兩個測試湖泊。

  ₃

第 1-2 節 研究方向及流程

基於上述二個研究目的，本文的研究方向及流程有四，分別為：

1. 推導 SA 所觀測到的海水面高度變化量計算重力位係數之公式，並利用快速 傅立葉轉換 (Fast Fourier Transform，FFT) 計算法，有效率的處理這些龐雜 的時空資料。相關的公式結合FFT 計算法撰寫成兩個 FORTRAN 程式，sha.f90 為可將大地起伏、重力異常值、水文高度及大氣壓力變化量計算轉換為相對 應的重力位係數程式；syn.f90 則相反的可由重力位係數計算出相對應之大地 起伏、重力異常值、水文高度及大氣壓力變化量，用以分析全球之重力時變 研究課題。研究地球時變重力場的課題，當其資料時間段越長及空間域的範 圍越大，將更能說明重力時變的情形及其相對應的環境變遷影響。在本研究

中所使用的資料主要為SA、steric 及 GRACE Level-2 重力係數資料。其中 SA 資 料 來 自 於 AVISO (Archiving ， Validation and interpretation of satellite Oceangraphic data)，資料期間段為 1992 年 10 月至 2007 年 12 月，而 GRACE 資料則分別來自於3 個資料處理中心，各為 GFZ (GeoForschungs Zentrum)、 CSR (Center for Space Research)及 JPL (Jet Propulsion Laboratory)，經簡單的平

均處理後，以平均值為計算資料，資料期間段則為 2002 年 10 月至 2007 年 12 月。資料的使用空間範圍，SA 資料為配合 steric 效應資料，所使用的資 料空間範圍為經度：0° ~ 360°、緯度：-70° ~ 70°間的海面區域。 GRACE Level-2 資 料 本 身 已 是 重 力 位 係 數 ， 資 料 範 圍則涵蓋全球。上述這些時空資料均非常巨量，即使依本研究以月為單位檔 案，相關處理計算量亦非常繁雜費時，因此一有效率的計算方法與程式是絕 對必要的。 2. 使用長期的 SA 資料，研究 GSLT 及其全球重力異常與大地起伏之變化速率。 影響海水面高度變化的主要因素可綜歸為3 項：(1)因溫度及鹽度的改變所造 成的海水膨脹或縮收效應稱為steric 效應， steric 效應並非質量的變化；(2) 因

  ₄

潮汐 (tides)、洋流 (current)、渦流 (eddy)及海洋、陸地、大氣間的水文循環 所造成的海水質量分佈變化；(3)冰河地殼均衡調整效應 (Glacial Isostatic Adjustment，GIA)。 SA 的觀測資料為海水面高度 (Sea Surface Height，SSH)， 當SSH 扣除某一參考平均海水面 (Mean Sea Level，MSL)後，可得 SLA，即 SLA = SSH－MSL。而依據理論 SLA 是包含了海水質量變化、steric 效應與 GIA 效應這三種影響分量，所以在解算僅海水質量變化所引起 GSLT 時，必

須扣除 steric 與 GIA 效應影響量，即改正後之海水面高度異常值 CSLA

(Corrected SLA) = SLA－steric－GIA。不過須特別提出的是，本文中暫未使

用任何的GIA 模式資料，所以最終所計算獲得之 CSLA 及後續衍生成果，均

內含GIA 效應。如各學者之研究指出 (Minster et al.，1999；董曉軍及黃城， 2000)，由於全球溫度升高所造成的暖化效應，除了造成海水面因 steric 效應 升高外，更造成大量融雪、融冰所化成的淡水 (fresh water) 流入海洋，而使 海水面高度上升。海水面高度的上升，相對應的是陸地面積的減少，尤其各 大洋中的小島，例如馬爾地夫 (Maldives)在未來數十年內會因海水面的上升 而消失，因此研究及監測GSLT 是近 20 年來重要的研究課題，尤其自 SA 技 術的發展後，提供了高頻率高密度及幾乎涵蓋全球海域的海水面高資料，令

其不僅能研究GSLT 問題，亦能由 SLA 顯而易見的看到 El Niño、La Niña 現

象、洋流與渦流之變化。本研究除透過長期的 SA 資料研究 GSLT 外，亦利 用程序1 之球諧係數解析及整合程式，由 CSLA 時變值計算出相關之重力位 係數，再由重力位係數整合計算出海水質量變化對全球重力場之變化影響量。 3. 藉由低階球諧係數之時變量，分析研究海水質量、地球質心坐標、極運動及 扁率等之變化量。重力位低階係數變化量，其意義為相對應之地球物理量的 變化。零階的重力位係數表示整體的質量，對整個地球而言，質量可假設為 守恆不變，零階係數變化量應為零，因此一旦零階係數變化量不為零，則代 表質量有所變動增減。以 SA 資料為例，CSLA 僅為海水區域範圍之觀測數 據，理論而言海水質量是變動的，是以可透過CSLA 所算得的零階重力位係

  ₅

數，推導得相應時間內之海水質量變動量。在質量守恆的假設下，如已知海

水質量變化量即可算得陸地及大氣水質量的變動量。一階係數的變化量

10

C

Δ 、ΔC₁₁、ΔS₁₁分別代表地球質心與地心地固框架( Earth Fixed Reference Frame，EFRF) 原點間的偏移量Δ 、z Δx、 Δ (Cretaxu et al.，2002)。針對y

這一階係數有兩種處理的方式，第一種為假設地球瞬時質心與地球參考重力

場重心一致，則ΔC₁₀ =ΔC₁₁ =ΔS₁₁ =0，GRACE 的重力位係數的結果即為此

種假設；第二種則是假設參考重力場重心與 EFRF 原點重合，允許一階係數

值不為 0，而這主要源由是來自於水質量時空分佈變化 (Cretaxu et al.，

2002) 。二階係數中常用的J₂代表地球扁率的變化量，即地球形狀的變化情

形。SLR (Satellite Laser Range)是一個公認可靠觀測J2變化量的一項技術，其

觀測數值代表整個地球扁率的總變化量。而同樣經由CSLA 資料所計算得之 2 J ，可瞭解海水質量變動對地球扁率之影響量。C21及S21係數與極運動有關， 吾人可經由C21及S21係數的變量ΔC₂₁及ΔS₂₁，計算得相對應的Δ 、x_p Δ (極y_p 運動分量)，而瞭解極軸的變動情形。 4. 湖水位高度變化量的計算監測。GRACE 衛星自 2002 年開始正常運作後，其 所提供的觀測資料，對研究大地測量、地球物理、水文學、大氣學等有關地 球環境相關學科有著極大的助益，尤其資料時間頻率的提升至以月為單位， 更是研究環境子系統間質量變化時間頻普的好資料，例如Velicogna 與 Wahr (2006) 利用 GRACE 衛星資料，估計 2002-2005 年間南極大陸冰層的融蝕年 速率約152±80 km3，相當於造成全球海平面年上升率為0.4±0.2 mm，而這流 失的質量大部分來自於南極大陸西部的冰層。理論上，小於數百年的時變重

力週期的主要變化因子來源是水質量的變動 (Wahr et al.，1998)，因此 GRACE

被大量且有效的用來研究有關水文系統的課目，例如 Chen et al. (2006、

2008)；Swenson et al. (2003、2006)等。理論上一個水文流域內，其水文系統 應處於平衡的狀態之下，即輸入的水質量會等於截流積存下來的水質量 (水

  ₆

庫、湖泊蓄水、地下水增補、窪蓄、土壤含水量等) 及輸出水質量 (蒸發、 蒸散、河流外流等) 之和，因此在封閉且面積在數百平方公里以上之水文系

統內的湖泊，應可利用 GRACE 資料來監測計算湖泊水位高度的變化量。本

文第四研究程序即使用GRACE 資料、SA (T/P 及 Jason-1) 與 CPC 水文模式 研究探討計算 Baika 及 Balkhash 的湖水位高度變化量。

第 1-3 節 各章節概述

本文共分成六大章節，各章節內容概述如下： 第1 章節為緒論：闡明本研究之研究緣起、方向與程序。 第2 章節為球諧函數係數的解析與整合：闡明推導 SA 所觀測到的海水高度變化 量計算重力位係數之公式，並結合 FFT 之原理推導球諧函數係數解析與整 合計算之公式，以有效率地處理及計算大量的SA、steric 及 GRACE 資料。 第3 章為資料來源與處理：闡述本文所使用的四種主要資料來源與處理程序，以 為後續第四、第五章節之成果計算分析之用。 第4 章為 SA 資料成果分析：分析 SA 及 steric 資料所計算之 GSLT 結果，並分 析全球大地起伏、重力異常之年變率、相位及振幅等之成果。另探討 SLA 與El Niño 之關係，且藉由海水質量時變量計算出之重力係數來分別探討低 階係數與相關地球物理意義之影響量。 第5 章為以測高及重力衛星分析水文變化：本章在於探討 GRACE 衛星重力觀測 量、SA 觀測量與 CPC 全球水文模式用於監測湖水位高度變化之可行性 第6 章為結論與建議：綜整本論文之研究成果，並對研究中不足之處提出說明與 建議。

  ₇

第 2 章 球諧函數係數的解析與整合

第 2-1 節 章節緒論

本章節是闡明球諧函數係數解析 (analysis) 與整合 (synthesis) 計算之公式 推導，並結合 FFT 之原理，以有效率地處理及計算大量的觀測資料。各小節內 容簡述如下：第2-1 節綜合簡述本章節及各小節之內容；第 2-2 節由球面上任意 函數 f(θ,λ)的球諧函數級數展開式推導起，含括正交函數及完全正規化理論， 推導說明普通諧函數係數與完全正規化函數係數之關係式。第 2-3 節則由依據

Heiskanen and Moritz (1985) 地球外部引力位 V 的球諧函數級數展開式推導起，

在考量地球非剛體變形負載Love (kn)係數下，推導大氣壓力異常值及水文高度異 常值之重力位係數計算式。第2-4 節是推導當球諧係數以複數型態表達時，該複 數型態可推演為FFT 的計算表達式，進而可使用 FFT 計算技巧於球諧係數分析 之應用；第 2-5 節 則為推導說明 FFT 計算技巧於球諧係數整合之應用。第 2-6 節則簡介wavelet 小波分析法的理論及其應用。第 2-7 節則為綜整結論本章內容。

第 2-2 節 完全正規諧函數係數與普通諧函數係數之關係式

根據理論 (Heiskanen and Moritz，1985)，在球面上任意的函數 f(θ,λ)都可 以球諧函數展開式來表示： )] , ( ) , ( [ ) , ( ) , ( 0 θ λ 0 0 θ λ θ λ λ θ nm nm n n n m nm nm n a R b S Y f = _∑∞ = _{∑ ∑} + = ∞ = = (2-1) 上式中a 、_nm b 為球諧係數，n、m 分別為階數 (degree) 及次數 (order)，而 _nm λ θ λ θ P m R_nm( , ) = _nm(cos )cos (2-2) λ θ λ θ P m S_nm ( , ) = _nm (cos )sin (2-3) 由於正交函數的定義是〝在單位球面上任何兩個不同的函數的乘積為零〞，所以 當n≠i 或 m≠j 時

∫∫

_σ R_nm (θ,λ)R_ij(θ,λ)dσ = 0 (2-4)

  ₈ ∫∫_σ S_nm(θ,λ)S_ij(θ,λ)dσ = 0 (2-5) 在任意情況下

∫∫

_σ R_nm(θ,λ)S_ij(θ,λ)dσ =0 (2-6) 而當n=i 且 m=j 時

[

]

∫∫

_σ θ λ σ = π₊ 1 2 4 ) , ( 2 0 n d R_n ，m=0 (2-7)

[

]

[

]

∫∫

_σ θ λ σ=

∫∫

_σ θ λ σ= π_{+ −}+ )! ( )! ( 1 2 2 ) , ( ) , ( 2 2 m n m n n d S d R_{nm nm} ，m≠0 (2-8) 在上面的公式中

∫∫ ∫ ∫

_σ= _λ2π_{=0 θ}π₌₀ 為單位球面積分的簡式，dσ =sinθdθdλ為單位球 體的面元素或稱立體角元素，此立體角元素被定義為單位球所對應的面。如將 (2-1) 式的兩端都乘以R_ij(θ,λ)，並對整個單位球進行積分，則(2-1) 式可改寫為

[

]

∫∫

_σ f(θ,λ)R_ij(θ,λ)dσ =a_ij

∫∫

_σ R_ij(θ,λ)2dσ (2-9) 依據正交函數定義的性質，(2-9) 式中除了當 n=i 且 m=j 的一項外，其餘項均等 於零。等號右邊分別以 (2-7) 及 (2-8) 式代入，可計算求得a 。同樣的計算方_ij 法，將 (2-1) 式的兩端都乘以S_ij(θ,λ)，則可求得

b

_ij，因此

∫∫

+ = σ σ θ λ θ π f P d n a_n ( , ) _n(cos ) 4 1 2 0 ，m=0 (2-10) σ λ θ λ θ λ θ π σ d S R f m n m n n b a nm nm nm nm

∫∫

⎭ ⎬ ⎫ ⎩ ⎨ ⎧ + − + = ⎭ ⎬ ⎫ ⎩ ⎨ ⎧ ) , ( ) , ( ) , ( )! ( )! ( 2 1 2 ，m≠0 (2-11) (2-10) 及 (2-11) 式為一般球諧係數的解算式，計算式除較複雜外，且當 m=0 或 m≠0 時計算公式亦不同，如將R 及_nm S 代以一相差一個常數因子的函數，即通_nm

稱之完全正規化諧函數R 及_nm S ，則依據 Heiskanen and Moritz (1985) 所示，_nm

nm R 、S 和_nm R 、_nm S 的關係式如下： _nm ) (cos 1 2 ) , ( 1 2 ) , ( _{0 0} 0 θ λ n θ λ n θ n n R n P R = + = + (2-12) ⎭ ⎬ ⎫ ⎩ ⎨ ⎧ + − + = ⎭ ⎬ ⎫ ⎩ ⎨ ⎧ ) , ( ) , ( )! ( )! ( ) 1 2 ( 2 ) , ( ) , ( λ θ λ θ λ θ λ θ nm nm nm nm S R m n m n n S R ，m≠0 (2-13)

  ₉ 正交關係對完全正規化諧函數，其式變為 1 4 1 4 1 2 ₌ 2 ₌

∫∫

∫∫

σ π σ π _σ Rnmd _σ Snmd (2-14) 因此 (2-1) 如改用完全正規化諧函數來展開，則 )] , ( ) , ( [ ) , ( ) , ( 0 θ λ 0 0 θ λ θ λ λ θ nm nm n n n m nm nm n a R b S Y f = _∑∞ = _{∑ ∑} + = ∞ = = (2-15) σ λ θ λ θ λ θ π _σ S d R f b a nm nm nm nm

∫∫

⎭ ⎬ ⎫ ⎩ ⎨ ⎧ = ⎭ ⎬ ⎫ ⎩ ⎨ ⎧ ) , ( ) , ( ) , ( 4 1 (2-16)

式中R_nm(θ,λ)=P_nm(cosθ)cosmλ，S_nm(θ,λ)=P_nm(cosθ)sinmλ，而P_nm(cosθ)為 完全正規化勒建德函數 (fully normalize legendre function)。綜整 (2-10)、(2-11)、 (2-12)、(2-13) 及 (2-14) 式，完全正規諧函數的係數和普通諧函數係數的關係為 1 2 0 0 + = n a a n n (2-17) ⎭ ⎬ ⎫ ⎩ ⎨ ⎧ − + + = ⎭ ⎬ ⎫ ⎩ ⎨ ⎧ nm nm nm nm b a m n m n n b a )! ( )! ( ) 1 2 ( 2 1 _{，m≠0 (2-18)}

第 2-3 節 水文高度及大氣壓力異常值重力位係數計算式推導

依據Heiskanen and Moritz (1985) 物理大地測量學乙書內的理論，地球外部

的引力位V (potential)可以球諧係數級數展開如下式：

(

)

⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ + ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + = _{∑ ∑}∞ =0 =0 ) (cos sin cos 1 n n m nm nm nm n E _{a m b m P} r a r GM V λ λ θ (2-19) 式中 GM 是引力常數與地球質量的乘積，r 為待求點至地心的距離，

a

_E為地球 平均半徑， a 、_nm b_nm則為待求之球諧係數，θ 為地心餘緯(co-latitude)，λ 為地 心經度（向東為正）。式中第一項 r GM 是表示齊次式 (homogeneous) 的球體及其 相對應的刻卜勒運動的引力位和，而第二項次雙重累加式，則是表示由擾動量所 形成的擾動位能T，即

  ₁₀

(

)

⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ + ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ =

_{∑ ∑}

∞ =0 =0 ) (cos sin cos n n m nm nm nm n E _{a m b m P} r a r GM T λ λ θ (2-20) 而距離l 倒數的球諧函數表示式為 )] , ( ) , ( ) , ( ) , ( [ 1 2 1 1 1 0 0 λ θ λ θ λ θ λ θ _{nm nm nm} n m nm n E n S S R R r a n r l ⎟⎠ ′ ′ + ′ ′ ⎞ ⎜ ⎝ ⎛ + = ∑ ∑ = ∞ = (2-21) 將上式代入 (2-20) 式，則a 、_nm b 的計算式為 _nm dM m m P r Ma n b a nm n n E nm nm

∫∫∫

⎭ ⎬ ⎫ ⎩ ⎨ ⎧ + = ⎭ ⎬ ⎫ ⎩ ⎨ ⎧ λ λ θ sin cos ) (cos ) 1 2 ( 1 (2-22) qds r d d q r d d dr r dV dM_{=ρ⋅ =ρ⋅} 2 _{sinθ θ λ=} 2 _{sinθ θ λ=} 2 _(2-23) 式中dM 為質量單元，dV 為體積單元，dS 為面積單元，ρ 為密度，q 為表面積單 元質量。將 (2-23) 式代入 (2-22) 式中，且設

a

_E＝r，則 ds m m P q M n a b a s nm E nm nm

∫∫

⎭ ⎬ ⎫ ⎩ ⎨ ⎧ + = ⎭ ⎬ ⎫ ⎩ ⎨ ⎧ λ λ θ sin cos ) (cos ) 1 2 ( 2 (2-24) 考量地表單位面積質量負載為 w w s _h g P q= =ρ (2-25) (2-25) 式中P 為大氣壓力異常值，_s h 為水文高度異常值，ρ_w w 表水密度，g 為正 常重力。將 (2-25) 式代入 (2-24) 式，則大氣壓力異常值及水文高度異常值所影 響的重力位係數計算式如下：

∫∫

_⎭⎬⎫ ⎩ ⎨ ⎧ + = ⎪⎭ ⎪ ⎬ ⎫ ⎪⎩ ⎪ ⎨ ⎧ s nm w w E h nm h nm _ds m m P h M n a b a w w λ λ θ ρ sin cos ) (cos ) 1 2 ( 2 (2-26)

∫∫

_⎭⎬⎫ ⎩ ⎨ ⎧ + = ⎪⎭ ⎪ ⎬ ⎫ ⎪⎩ ⎪ ⎨ ⎧ s nm s E P nm P nm _ds m m P P Mg n a b a s s λ λ θ sin cos ) (cos ) 1 2 ( 2 (2-27) 如再考慮地球非剛體的彈性變形的負載Love 係數 kn，則 (2-26)、(2-27)可表示如 (2-28)、(2-29) 式

  ₁₁

∫∫

_⎭⎬⎫ ⎩ ⎨ ⎧ + + = ⎪⎭ ⎪ ⎬ ⎫ ⎪⎩ ⎪ ⎨ ⎧ s nm w n w E h nm h nm _ds m m P h M n k a b a w w λ λ θ ρ sin cos ) (cos ) 1 2 ( ) 1 ( 2 (2-28)

∫∫

_⎭⎬⎫ ⎩ ⎨ ⎧ + + = ⎪⎭ ⎪ ⎬ ⎫ ⎪⎩ ⎪ ⎨ ⎧ s nm s n E P nm P nm ds m m P P Mg n k a b a s s λ λ θ sin cos ) (cos ) 1 2 ( ) 1 ( 2 (2-29) 如所知，球的體積公式為 3 3 4 R π ，而地球質量 M=ρ_EV = 3 3 4 R E π ρ (2-30) (2-30) 式中ρ 為地球平均密度，約為 5517 kg/m_E 3_{。將 (2-30) 式代入 (2-28)、(2-29)} 式，則 (2-31)、(2-32) 式之水文高度異常值及大氣壓力異常值所產生的重力位係 數變化計算式即為 Whar et al. (1998) 所推導應用之公式：

∫∫

_⎭⎬⎫ ⎩ ⎨ ⎧ + + = ⎪⎭ ⎪ ⎬ ⎫ ⎪⎩ ⎪ ⎨ ⎧ s nm w E E n w h nm h nm ds m m P h n a k b a w w λ λ θ πρ ρ sin cos ) (cos ) 1 2 ( 4 ) 1 ( 3 ( 2- 31 )

∫∫

_⎭⎬⎫ ⎩ ⎨ ⎧ + + = ⎪⎭ ⎪ ⎬ ⎫ ⎪⎩ ⎪ ⎨ ⎧ s nm s E E n P nm P nm _ds m m P P g n a k b a s s λ λ θ πρ sin cos ) (cos ) 1 2 ( 4 ) 1 ( 3 (2-32)

第 2-4 節 FFT 計算技巧於球諧係數解析之應用

在第 2-2 小節中之 (2-16) 式為正交球諧係數的連續積分求解式，如將該公 式以複數形式來表達，則 λ θ θ θ λ θ π λ σ d d e P f b i a C_{nm nm nm} ( , ) _nm(cos ) im sin 4 1

∫∫

⋅ = + = (2-33) 在實際的觀測量上，函數 f(θ,λ)都為不連續的觀測點，所以 (2-33) 式一般而言 是無法直接計算的，而必須將 (2-33) 式改寫成不連續的公式。現若給定在 θ、λ 方向等間隔之網格資料 (Δθ×Δλ)，則 (2-33) 式的不連續公式可寫為 (Hwang et al.，2006)

  ₁₂

∑ ∑

∫

∑ ∑

∫

∫

∑ ∑

∫ ∫

− = − = − = − = − = − = + + + + + ⋅ = = = 1 0 1 0 1 0 1 0 1 0 1 0 1 1 1 1 1 ) , ( 4 1 ) ( ) , ( 4 1 ) (cos ) , ( 4 1 M k N l im k nm l k n M k N l t t im nm l k n M k N l im nm l k n nm l l k k l l k k l l d e P I f q d e dt t P f q d d e P f q C λ λ λ λ λ λ θ θ λ λ λ λ λ θ π λ λ θ π λ θ θ λ θ π (2-34) 式中M =(π/Δθ+1)為緯度方向的網格數，N=(2π/Δλ+1)為經度方向的網格 數，t =cos

θ

、dt=sinθdθ 、t_k =cos((k−1)Δθ)，λ_l =(l−1)Δλ， k

nm

P

I 為締和勒

建德函數 (Associated Lengender Function) 的積分 (Paul，1978)， f(θ_k,λ_l)為平

均值，

q

_n為平滑因子，主要是為使得 (2-34) 式不連續公式的結果更趨近於積分 式的結果 (Rapp，1989)。根據 Colombo (1981，p. 76) 的建議，平滑因子

q

_n的設 定機制為 L n q L n L q L n q n n n n n > = < < = ≤ ≤ = , 1 3 / , 3 / 0 , 2 β β (2-35) 其中N 為所展開最大的階數，而β 為 Pellinen 平滑運算元，其定義為 _n

[

(cos ) (cos )

]

1 2 1 cos 1 1 0 1 0 1 0 ψ ψ ψ βn = ₋ ⋅ _n+ Pn− −Pn+ (2-36) 式中ψ 為球帽之半徑，其面積與₀ f(θ_k,λ_l)網格面積相同。 在(2-34)式中

_∫

l+1 l d eim λ λ λ _{λ 可及其分為 m=0 及 m≠0 兩種情形討論。當 m=0 時}

λ

λ

λ

λ λ λ

₌

_⋅

_Δ

₌

_Δ

∫

+

1

1 l l

d

e

im (2-37) 當m≠0 時

  ₁₃

[

]

N iml im iml im iml l im im im im im

e

m

e

i

e

m

e

i

e

e

m

i

e

e

m

i

e

im

d

e

l l l l l l / 2 ) 1 (

)

1

(

)

1

(

)

(

1

_{1 1} 1 π λ λ λ λ λ λ λ λ λ λ λ λ λ

_λ

−

−

=

×

−

−

=

−

−

=

−

−

=

=

Δ Δ Δ Δ Δ + + + +

∫

( 2 - 3 8 ) 於 (2-38) 式中

[

sin( ) (1 cos( ))

]

1 ) 1 sin (cos ) 1 ( λ λ λ _{λ λ} Δ − − Δ = − Δ − Δ − = − − Δ m i m m m m i m i m e i im (2-39) 整合 (2-37)、(2-38) 及 (2-39) 式，並引入符號gm，且

[

]

⎪⎩ ⎪ ⎨ ⎧ ≠ Δ − − Δ = Δ = 0 , )) cos( 1 ( ) sin( 1 0 , m m i m m m g_{m λ λ} λ (2-40) 結合 (2-34) 式及 (2-40) 式，則 (2-34) 式可改寫為

[

]

∑

∑ ∑

∑∑

− = − − − = − = − = − = − + = ⋅ = ⋅ ⋅ ⋅ = 2 / ) 1 ( 0 1 0 1 0 / 2 1 0 1 0 / 2 ) ( ) 1 ( ) ( 4 ) , ( 4 ) , ( 4 1 M k k M m n k k nm n m M k N l N iml l k k nm n m M k N l N iml m k nm l k n nm m f m f P I q g e f P I q g e g P I f q C

π

λ

θ

π

λ

θ

π

π π (2-41) 上式中 f(θ_k,λ_l)是取計算區塊四個頂點值的平均值，計算式如下： )] ) 1 ( , ( ) ) 1 ( , ) 1 ( ( ) , ) 1 (( ) , ( [ 4 1 ) , ( λ θ λ θ λ θ λ θ λ θ Δ + Δ + Δ + Δ + + Δ Δ + + Δ Δ = l k f l k f l k f l k f f _{k l} (2-42) 令k=N-1 為沿經度方向的區塊數，p=l-1 所以

∑

− = = 1 0 / 2 ) , ( ) ( k p k mp i k k k m f e f θ λ π ，m=0,…,k-1 (2-43) 比較 (2-43) 與 (2-31 )式，本研究即利用上式之 FFT 演算法撰寫一解析球諧係數 之 FORTRAN 程式：sha.f90，程式中是設定同時由最北及最南呈帶狀往赤道計 算，有關本程式的相關功能與用法請參考附錄A。

  ₁₄

第 2-5 節 FFT 計算技巧於球諧係數整合之應用

上一小節是推導球諧係數解析的 FFT 計算公式，反之當已有球諧係數時， 既可以使用逆的FFT 計算式整合係數，反演計算相關觀測量。根據(2-16)式，利 用已知球諧係數計算網格點上之函數值公式可表示為(設網格大小為Δθ×Δλ):

∑ ∑

= = Δ Δ + Δ Δ = Δ Δ l n n m nm nm nm nmR k l b S k l a l k f 0 0 )] , ( ) , ( [ ) , ( θ λ θ λ θ λ (2-44) 式中k=0,…M，l=0,…N-1。需注意的是，上式中當經度為 0°與 360°時，f 函數值 相等 ，即 f(kΔθ,NΔλ)= f(kΔθ,0)。為了能使用 FFT 的計算技巧，(2-44) 式必 須改寫成 (2-45) 式：

(

)

1 ,...., 0 , ,.... 0 ) sin( ) cos( ) sin( ) ( ) cos( ) ( ) , ( 0 0 − = = Δ + Δ = Δ ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ _Δ + Δ ⎩ ⎨ ⎧ ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ Δ = Δ Δ

∑

∑

∑ ∑

= = N l M k ml S ml C ml b k P ml a k P l k f L nm nm L m n nm nm L L m n nm nm λ λ λ θ λ θ λ θ (2-45) 在上式中必須注意S0=0，令 m=L+1,…,N-1 時 Cm=Sm=0，且 2 m m m iS C B = − (2-46) 則 1 ,..., 0 ), Re( 2 Re 2 ) , ( 1 0 / 2 − = = ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ = Δ Δ ∑− = N l P e B l k f n m N mlk i m π λ θ _(2-47) 式中 Re 表示複數形式的實數部分，而 P 則可利用 FFT 演算法沿平行圈網格計 算而得，是以於計算過程中另需注意Δθ是否能整除90°，即 90°須為Δθ 的倍數。 此外，因締合勒建德函數滿足 (2-48) 式， ) ( ) 1 ( ) ( t P t P n m _nm nm − = − − (2-48) 因此於計算締合勒建德函數值時，只須計算北半球之值，而南半球之締合勒建德 函數值則可直接引用 (2-48) 式求得 (Colombo，1981)。於實際計算時，本研究 先將大地起伏與重力異常值計算式改寫成 (2-49) 及 (2-50) 式

  ₁₅

[

]

∑ ∑ ∑ ∑ ∑ ∑ = Δ Δ = = = = ∞ = + = ⎪⎭ ⎪ ⎬ ⎫ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + ⎥ ⎥ ⎦ ⎤ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ⎩ ⎨ ⎧ ⎢ ⎣ ⎡ = + × = Δ L m N m N m L m n nm nm n E n nm nm n E L m L m n n nm nm nm nm n m n n E n m S m C m K P r a rg GM m J P r a rg GM S K R J r a rg GM r N 0 0 0 0 )) sin( ) cos( ( ) sin( ) ( ) cos( ) ( ) , ( ) , ( ) ( ) , , ( λ λ λ θ λ θ λ θ λ θ λ θ (2-49)

[

]

∑ ∑ ∑ ∑ ∑ ∑ = Δ Δ = = = = ∞ = + = ⎪⎭ ⎪ ⎬ ⎫ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ − + ⎥ ⎥ ⎦ ⎤ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ⎩ ⎨ ⎧ ⎢⎣ ⎡ ₋ = + × − = Δ L m g m g m L m n nm nm n E nm nm n E L m L m n nm nm nm nm n m n n E m S m C m K P r a n r GM m J P r a n r GM S K R J r a n r GM r g 0 2 0 2 0 0 2 )) sin( ) cos( ( ) sin( ) ( ) 1 ( ) cos( ) ( ) 1 ( ) , ( ) , ( ) )( 1 ( ) , , ( λ λ λ θ λ θ λ θ λ θ λ θ (2-50) 再利用此二式之 FFT 演算法撰寫一整合球諧係數之 FORTRAN 程式：syn.f90， 而有關本程式的相關功能與用法請參考附錄B。

第 2-6 節 Wavelet 轉換計算分析與應用

小波轉換 (Wavelet Transform，WT) 分析應用於地球科學領域，最早可推 溯自1980 年地震訊號的分析 (Morlet et al.，1982)，隨後由 Grossmann and Morlet (1984) 與 Goupillaud et al. (1984)給予正式定義，並由 Meyer (1992)、Mallat (1989)、 Daubechies （1992）、 Chui ( 1992)、 Wornell (1995) 與 Holschneider (1995) 等學者應用於其他領域的研究。近年小波轉換已成為地球科學領域分析 局部變化的有效工具，例如大氣冷鋒的研究 (Gamage and Blumen，1993)、熱帶 對流變化之探討 (Weng and Lau，1994)、ENSO (El Niño–Southern Oscillation， ENSO) 周期變化之分析 ( Wang and Wang，1996)、海洋溫度變化 (Flinchem and Jay，2000)、海水面高度變化 (Percival and Mofjeld，1997)、陸地溫度長期氣候 變遷 (Sharp，2003)等研究。

  ₁₆ 傅立葉轉換適用於一般常規性訊號頻譜轉換分析，對局部異常訊號的解析 並不適切，其往往無法針對特點給出特別訊號，例如降水在時間和許多其他地球 物理現象的相互影響過程訊號。小波轉換基本上是使用兩種方法來解析研究地球 物理過程或信號 (Kumar，1995）： (1)以積分為核心，分析擷取出有關訊號。(2) 依據擷取之訊號，解釋資料間之特徵。以小波轉換作為分析信號核心的優點，在 於它使信號的局部特徵與所使用的時間序列間隔能夠相匹配對應，即可同時具有 廣泛的大尺度和小尺度精細的特點，即其優點在於分析視窗可變動，可用較長的 空間域或時間來獲得低頻訊號，或以較短的空間域或時間來取得高頻訊號。這種 特性對非穩態、短週期分量及不同尺度特徵或訊號之分析相當有效。 小波轉換函數以積分形式 f(t)可定義為(Kumar，1995）：

∫

−∞∞ Ψ > = ( ) ( ) 0 ) , (λ t f u _λ, u du λ Wf _t (2-51) 其中 ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ − Ψ = Ψ λ λ λt u u t 1 ) ( , (2-52) 即為小波函數的母函數。在式中λ 是比例參數，t 是平移參數，Ψ_λ,t(u)是Ψ_λ,t(u) 的共軛參數。改變λ 值會有擴張( λ >1)或縮減( λ <1) )Ψ(t 的效應，更改 t 值則 具有分析不同點 t 周圍函數 f(t)的效果。 當λ 比例增加時，小波變得更分散 與長時間 f(t) 相關，反之亦然。 因此小波轉換類似于變焦鏡頭的功能，能提供 有彈性的時間尺度視窗。 上式小波轉換因為比例λ 和時間參數 t 均為連續值， 故稱為連續小波轉換。 在(2-52)式中， λ 1 為確保Ψ(u)與 ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ − Ψ λ t u 之範數(norm)相等，所以小波母 函數Ψ(u)必須具有收斂性，滿足容許條件(admissibility condition)及同值條件 (equivalently condition) (Mallat，1998)，即為下列二式:

∫

∞ ∞ = Ψ -0 ) ( duu (2-53)

  ₁₇ 0 ) ( 2 < Ψ

∫

∞ ∞ − ω ω ω d (2-54) 式中ω 為頻率變數。在使用及選擇小波函數時，必須考慮以下幾個要素(Andrew

and Gao，1996 、Torrence and Compo，1998):

(1) 平滑性：小波函數之平滑性是由其存在的導函數數目決定，數目越多則越平 滑，才能有效分析訊號的特性。 (2) 時間及空間的局部性：小波函數的有效頻譜寬度與其時間及空間之局部解析 能力有關，有效頻譜寬度越小，其局部解析能力越佳，且其與平滑性成正比， 即寬度越大越平滑。 (3) 正交性：對大尺度分析而言，非正交小波分析是多餘的，而在正交小波分析 中，每一尺度的捲積數與其小波基底寬度成比例，最能完全重現訊號，且是 用於平滑及連續變化之時序訊號分析。 (4) 對稱性：對稱性的小波函數可避免相位位移，即小波係數相對於原始訊號式 無移位的。 (5) 力矩因素：力矩消散次數越多，平滑性越低，但越能表示高次多項式的訊號。 含有力矩消散(vanishing moments)M 之小波母函數Ψ(u)，必須滿足下式: ∞

_∫

∞ − = = Ψ -1 ..., 2 , 1 , 0 0 ) (u du m M um (2-55) (6) 複數或實數:複數小波函數能同時獲得振幅與相位之資訊，其是用於擷取周期 震盪性訊號。實數小波函數僅能獲得單一訊號分量，適用於偵測尖峰值或不 連續點訊號。 (7) 訊號形狀:不同資料的解析訊號特徵不同，為能充分了解資料所代表的內涵， 選擇適當的小波函數有其必要性，例如訊號波動劇烈或為階梯狀，則可考慮 採用Harr 小波函數。

  ₁₈

現今常用的連續小波轉換為 Mexican hat 及 Morlet ( Daubechies，1992)，

其中Morlet 因其實數部具對稱性，頻率域為呈高斯分佈曲線，具有較佳的局部 化功能，對頻率峰值的解析極有效益。Morlet 小波函數定義為(Kumar，1995）： 5 ) 1 ( ) ( 14 2 2 ₀ 2 0 ≥ − = Ψ _{t π}− _{e t e}−iwt_e−t _{w (2-56)} 式中ω₀為時間域上單位寬度高斯曲線內的震盪波數，ω₀值愈大，則於時間域上 之單位寬度高斯曲線內的震盪波數增加，局部化能力減弱，解析力降低。而頻率 域其視窗向高頻移動，平滑化效果減低，局部化能力增加，解析力提升。反之， 若ω₀值愈小，則空間域解析能力增強，而頻率域解析能力降低。通常於選擇時 是以母函數Ψ(u)的前兩個最大波峰值之比值近似0.5 為原則來選定ω₀值，即 3364 . 5 2 ln 2 12 0 ⎟ ≅ ⎠ ⎞ ⎜ ⎝ ⎛ =π ω (2-57) 於實際應用時大多取ω₀ ≥5，此即為(2-56)式中 Kumar (1995）所定義應用之值。 在(2-56)式中的第二項非常小，於實際應用上可以忽略 (Daubechies，1992)。Morlet 為複數形態的小波函數，在本文中所需分析的訊號大多為實數型態，且為偵測間 峰值及不連續點，故於實際應用時僅取其實部，即(2-56)式可改寫為下式 (Michel et al.，1996) : ) 5 cos( ) ( 14 2 2 t e t = − −t Ψ π (2-58) 綜合(2-57)與(2-58)式，並根據 Torrence and Compo (1998) 之推導，Morlet

小波函數尺度參數λ 與傅立葉轉換周期 P 之關係，可計算如下: λ πλ ω ω πλ 1.232 5 2 5 4 2 4 P 2 2 0 0 ≅ + + = + + = (2-59) 本文即使用(2-59)式之比例參數作為分析解釋本文之結果。於實際應用上， 吾人是使用 MATLAB 軟體內的 Morlet 小波模式來進行資料週期分析，並結合 GMT 軟體繪製輸出小波頻譜分析圖。相關計算操作程序請參考附錄 C。

  ₁₉

第 2-7 節 本章小結

綜整結論本章內容如下： (1) 本章節最重要的論述與貢獻在於球諧係數解析與整合程式的推導及撰寫，此 兩個程式均結合FFT 之原理，能有效率地處理及計算大量的觀測資料，極大 助益本文的研究。基本上本文的各項成果均源於此兩程式，第4 章及第 5 章 成果的分析均構基於此。 (2) 本文在考量地球非剛體變形負載係數下，推導因水文高度異常值及大氣壓力 異常值所影響之重力為係數計算式，此一計算公式，已納入syn.f90 程式中， 為第5 章 GRACE 監測湖水位高度成果分析計算之用。 (3) 以小波函數作為分析信號核心的優點，在於分析視窗可變動，可用較長的空 間域或時間來獲得低頻訊號，或以較短的空間域或時間來取得高頻訊號。這 種特性對非穩態、短週期分量及不同尺度特徵或訊號之分析相當有效。

  ₂₀

第 3 章 資料來源與處理

第 3-1 節 章節緒論

本研究所使用的資料種類主要有四種，分別為SA、steric、GRACE 重力位 係數及 CPC-LDAS 水文模式資料。本章節內容即在闡述說明這四種資料的來源 及處理程序，以期有完善可靠的資料品質，以供後續研究之用。各小節內容如下： 第3-1 小節仍是概述本章各小節之內容；第 3-2 及第 3-3 節則分別說明 SA 與 steric 資料的處理方法與程序。因為 SA 結果是含括了 steric 效應，而 steric 效應是無 關質量變化的，因此本文研究海水質量所引起的海水面高度變化及後續低階重力 位係數之研究，均須使用扣除steric 效應後的 SA 改正資料 CSLA。第 3-4 小節是 討論GRACE 資料的處理情形。由於原有 GRACE 資料品質並不佳，各學者的研 究指出，GRACE 資料除含有偶然誤差外，尚有著明顯的系統性誤差 (Chamber， 2007)，因此使用 GRACE 資料前必須經過前處理，方能應用於後續之重力場變 化研究。第3-5 節 CPC-LDAS 水文模式在本研究中是用以和 GRACE 計算成果比 較驗證之用，節內是說明該資料特性。第3-6 節為本章之綜合結論。

第 3-2 節 AVISO 衛星測高資料

SA 資料的處理技術，經多年來的發展已相當成熟，相關的官方資料中心均 供 應 各 種 衛 星 及 不 同 格 式 的 SA 資 料 ， 例 如 JPL 的 PO.DAAC (Physical Oceanorgphy Distributed Active Archive Centre) 及 AVISO。本研究所使用的 SA 資 料 是 來 自 於 AVISO ， 該 資 料 中 心 不 僅 提 供 個 別 衛 星 觀 測 資 料 ( 含 括 Topex/Poseidon 、 GEOSAT 、 ERS-1 、 ERS-2、 Jason-1、 GFO 、 Envisat

等) ，亦提供多顆衛星之整合性資料，其中的 DUACS (Data Unification and Altimeter Combination System) 即為提供整合型資料的系統。此系統提供近即時 ((Near) Real Time，(N)RT) 及延時 (Delayed Time，DT) 兩大類資料，便於不同 需求之單位或個人應用。DUACS 系統資料處理的詳細程序與方法請參閱該中心

  ₂₁ 文件。

本 研 究 中 吾 人 使 用 DT-MSLA (Maps of Sea level Anomalies) 的 Ref (Reference) 1/4° 全 球 網 格 測 高 資 料 ， 該 資 料 是 使 用 兩 顆 衛 星 (Jason-2 / Envisat or Jason-1 / Envisat or Topex/Poseidon / ERS)的觀測資料經同質化 (homogeneous) 後，再將這些多衛星任務資料整合所得之成果。每一個檔案為涵 蓋7 天時間的多衛星整合型資料，資料涵蓋範圍全球，使用時間段自 1992 年 10 月至2007 年 12 月，共計下載 794 個檔案。在本研究中所需求的最小時間段單位 為〝月〞，網格大小均為1°×1°，是以所有使用到的資料均必須歸化計算為以 1°×1° 網格及月為儲存單位的檔案。由 AVISO 所下載之 DT-MSLA 資料，吾人計採 用兩個簡單的程序將其改算為所需之檔案形式: (1) 首先將同月份時間內的檔案，利用一般算術平均法歸算為月平均資料。 (2) 其次使用 GMT 軟體內的 Gauss 濾波法重新取樣，將原 1/4°之網格取樣為 1°×1°網格。 經此兩程序重新處理後，共計獲得183 個以月為時間單位的檔案，資料期間段仍 為1992 年 10 月至 2007 年 12 月。 AVISO 均有提供其各項成果的誤差分析量，其中本研究所使用的 DT-MSLA 成果，於經多項誤差改正(儀器誤差改正、環境誤差改正、衛星軌道改正等) 及 使用spline 演算法平滑濾波去除長波長誤差後，經評估在深海區 SLA 的精度約 3 cm，近岸地區或淺海地區的精度約 7 cm。圖 3-1 是 AVISO 所提供的太平洋亞熱 帶地區(N20°~40°，E100°~240°) 2007 年 1 月 18 日的 SLA 精度評估圖。在圖中的 不同顏色是表示觀測量與其誤差量相較的百分比，其中深海地區其比值差量在5 ％～6％以下，而近岸區可達 20％以上。本研究所使用的 SA 資料，如上所述是 重新取樣過的，為評估新的SA 資料品質，吾人以 1997 年整年 12 個月的資料來 測試評估，評估測試結果如表3-1。表中數據為各月份內部精度結果，結果顯示

  ₂₂ 格之 SLA 標準偏差約 0.4 cm。不過由於 1°×1°的資料是經由 Gauss 濾波平滑化 過，因此，此差異量並不能用以說明重新取樣後的資料品質優於原資料，只能說 明重新取樣後的1°×1°資料內部精度並未降低。圖 3-2 (a)、(b) 分別為 1997 年 11 月 1/4°×1/4°及 1°×1°網格全球 SLA 圖，圖中顯示兩者並無明顯差異性，此亦能 證明經處理後之1°×1°網格 SLA 資料之內部精度並未降低。另外，由圖中亦顯示 1997 年 11 月南美洲厄瓜多爾及秘魯西部海岸外的聖嬰現象的 SLA 值可高於 50 cm 以上。 表3-1：1997 年各月份 AVISO DT-MSLA 1/4°網格及本文重新取樣為 1°×1° 網格後之SLA 標準偏差表 Month 1 2 3 4 5 6 7 8 9 10 11 12 1/4°網格 Std(cm) 5.78 5.90 5.97 5.54 5.53 5.80 5.95 6.37 6.91 7.12 7.53 7.37 1°×1°網格 Std(cm) 5.38 5.48 5.58 5.15 5.11 5.38 5.56 5.99 6.54 6.79 7.19 7.05 較差(cm) Std 0.40 0.42 0.39 0.39 0.42 0.42 0.39 0.38 0.37 0.37 0.34 0.32

  ₂₃ 圖3-1：太平洋亞熱帶地區 (N20°~40°，E100°~240°) 2007 年 1 月 18 日的 SLA 觀測量與其誤差量的百分比精度評估圖(摘自 AVISO) (a) (b) 圖 3-2：1997 年 11 月全球 SLA 圖， (a)及(b)的資料網格分別為 1/4°×1/4°及 1°×1°

  ₂₄

第 3-3 節 steric 熱容效應資料

steric 是造成海水面上升的主要因素之一。 SA 的觀測數據內是包含著此種

效應，因此要研究純海水質量增減所造成的海水面高度變化量，就必須扣除steric

效應影響量，有關 steric 效應的影響量及研究多有學者著墨(Cabanes et al.， 2001a； Antonvo et al.，2002；Levitus et al.，2005；Ishii et al.，2005)。在衛星 觀測技術尚未發展成熟前，海水溫度及鹽度的量測記錄，大多依靠船測或浮標等 方式，所獲得的資料量及區域範圍性，相對於海洋面積及時間頻率而言，是相當 匱乏的，因此steric 效應的計算成果就有其侷促性，是以使用於改正 SA 資料， 反而可能造成CSLA 品質不可靠。然近 20 年來，由於衛星觀測技術的發展，再 配合原有船隻等觀測技術及觀測量，令steric 的計算結果越加穩定可靠。在本研 究中，吾人是參考使用國立成功大學郭重言博士(private communication) 所提供 的steric 資料，用以獲得 CSLA，以供後續研究之使用。 根據理論steric 效應為海水溫度效應(thermosteric)及鹽度變化效應(halosteric) 的綜合影響量。海水溫度變化會造成海水體積的膨脹或收縮，即當海水溫度升高 時會造成海水面升高，反之當海水溫度降低時，將造成海水面降低，這種效應是 無關質量變化的。鹽度的變化同樣會造成海水體積的膨脹或縮收，當鹽度增加 時，海水密度增加，造成海水面高度降低；而當海水鹽度降低時，海水密度減少，

此時造成海水面升高。計算steric 海水面高度 (Steric Sea Level，SSL) 變化量， 非僅使用海水表面觀測數據，應較嚴謹的考量海水垂直體積的整體變化總和量。

本研究所使用的SSL 資料為郭重言博士所提供，其所使用的海水表面資料(Ocean

Surface Data ， OSD) 來 自 於 NOAA(National Oceanic and Atmospheric Administration)的 WDA (World Data Atlase) 資料庫，由 Levitus 等人所編輯。而 SSL 的計算式如下： dz S T z t S T z t t h SSL [ ( , , , , , ) ( , , , , , )] 1 ) , , ( 0 λ ϕ ρ λ ϕ ρ ρ λ ϕ η η

_∫

− − = Δ (3-1)

  ₂₅ 式中η 為海面高度，h 為最大海水深度，ρ₀為參考密度 (1027 kg/m3)，t 為計算 時刻，z 為海水深度垂直距離，T 及 S 則為計算位置及時刻之海水溫度及鹽度， T 及S 則為計算時間段內的平均溫度及鹽度。目前本研究有關的 SSL 資料，同 樣是以月為時間單位的檔案，資料網格亦為1°×1°，時間段同為 1992 年 10 月至 2007 年 12 月。圖 3-3 為 1997 年 11 月的 steric 效應所造成的海水面高度異常圖， 圖中清楚地顯示南美洲厄瓜多爾及秘魯西部海岸外的聖嬰現象。 圖3-3：1997 年 11 月 steric 效應之海水面高度異常

  ₂₆

第 3-4 節 GRACE 資料

GRACE 的資料吾人是採用了 CSR、JPL 及 GFZ 3 個科學資料系統 (Science Data System, SDS) 的 GSM-RL04 資料，並依據 Chambers (2007) 的研究建議， 平均了 3 個 SDS 的觀測資料，資料期間自 2002 年 8 月至 2007 年 12 月，扣除資 料品質較差的月份，共計使用了 60 個月資料。 GRACE 資料除含有偶然誤差 外，尚有著明顯的系統性誤差，經多位學者研究指出，GRACE 觀測數據在處理 計算為Level-2 資料的過程中，所使用各種改正模式（例如，海潮、地潮、大氣、 水文等）所內含的模式誤差及 GRACE 衛星本身儀器誤差的影響，致使 Level-2 的重力位係數含有系統性的誤差，依據這些含有系統性誤差的重力位係數所繪製 出的全球重力異常圖或大地起伏異常圖，均可明顯可見南北方向的 ”strips” (Swenson and Wahr，2007)，此外 Swenson 亦研究指出 GRACE 的重力位係數 中亦含有係數間的相關性誤差 (correlated Errors)，因此在使用這些重力位係數 前，應先降低”strips”及相關性誤差的影響量 (Swenson and Milly， 2007； Chen et al.， 2005；Chambers et al.， 2004；Tapley et al.， 2004；Wahr et al.， 2004)。 一般而言，降低”strips” 誤差效應的方法大多使用空間平滑化理論，例如平均高 斯濾波法 (isotropic Gaussian filter ) (Wahr et al.，1998)、非平均高斯濾波法 (nonisotropic filter based on the calibrated error spectrum) ( Han et al.，2005)、最佳 濾波法 (the optimal filters based on a priori estimates of signal and measurement error variances) ( Seo and Wilson，2005)。而在消除重力位係數間相關性誤差方法 有多項式擬合法 (Polynomial Fit Menthod，PFM) (Swenson and Wahr，2007)。

在降低重力位係數間相關性誤差方法上，本文是採用Swenson 和 Wahr (2007) 的多項式擬合法，處理時是使用5 階多項式。在消除 strips 方法上，則是使用 Jekeli (1981) 的高斯平均濾波演算法 (Gaussian mean filter algorithm)，這高斯平均濾波 演算法名稱的來自于鐘型 (bell-shaped) 常態機率密度函數，而其權函數為球面 角距Ψ，相關公式的推導請參考 Jekeli (1981)。公式 (3-2) 為高斯平均濾波演算