國立臺灣大學理學院海洋研究所 碩士論文
Institute of Oceanography College of Science
National Taiwan University Master Thesis
高頻雷達觀測流場補缺值問題之探討
─以台灣東北部海域為例
A Study on Data Filling from Incomplete Dataset of HF Radar Measured Ocean Currents
─A Case Study of the Flow Field Northeast of Taiwan
陳俞彣 Yu-Wen Chen
指導教授:王冑 博士、詹森 博士 Advisor: Joe Wang, Ph.D.; Sen Jan, Ph.D.
中華民國 102 年 7 月
口試委員會審定書
誌謝
首先感謝指導教授王冑老師的耐心指導,無論是在學術研究、課業或是完成 論文等方面都給予我極大的幫助,在研究態度、邏輯思考等方面也從老師的教誨 中受益良多。感謝口試委員陳慶生老師、詹森老師以及楊穎堅老師提供許多珍貴 的意見與錯誤的指正,使本文更加完善。另外還要感謝海洋所的老師們在課業上 的教導,讓我學習到許多專業知識並奠定研究的基礎,由衷感謝。
感謝實驗室蔡雅鈴學姊、郭怡君學姊在日常生活以及學業上的照顧與解惑,
也要感謝122 實驗室的各位學長及助理在平日的協助,特別謝謝郭天俠學長在許 多事情的幫助與鼓勵,尤其是時常督促我的研究進度,使本文能順利完成,非常 感謝。同時要感謝我在碩士生涯中的好伙伴王心秀與廖康淵,不管是在課業、研 究或是平日生活方面,都給予我相當多的幫助以及陪伴。還要感謝其他陪我度過 這段研究生活的人,感激不盡。
最後感謝我所有的家人以及朋友們,時常關心、鼓勵我,並陪伴、支持我完 成學業,謝謝。
摘要
高頻雷達是近代廣用的一種即時監測海流的遙測工具,但高頻雷達訊號常因 環境雜訊以及電離層等干擾影響,會造成在觀測區內出現海流觀測值缺漏的情況,
因此如何有效填補缺漏值即成為海流監測即時作業一個重要的課題。本文使用臺 灣大學海洋所、海軍軍官學校以及海軍大氣海洋局共同合作在臺灣東北部海域建 置的CODAR 雷達系統觀測資料做為研究的根據。
我們選用實向量經驗正交函數(real-vector EOF)及 Karhunen-Loève 展開法 (KLE)兩種統計方法探討流場填補缺漏值問題。這兩種方法均是根據長期觀測資料 之變量統計求出對應之特徵向量(eigenvector)作為模組(mode)基底,而觀測資料在 各基底之投影即為該模組之振幅(amplitude)。經比較分析後,我們選用前 20 組模 組(可解釋 96%以上的總變異量),分別使用最小平方法(least square)與迭代法 (iteration)以求解不完整資料情形下各模組之振幅,並以此進行不同個數缺漏值之 填補實驗,以檢驗填補效果及填補誤差,藉以推測使用不同方法之缺漏個數上限。
填補實驗的統計結果顯示,以最小平方法進行填補實驗時,在缺漏個數小於 資料點總數之 57%的情況下,其誤差隨缺漏點個數的變化並不大,若以迭代法搭 配實向量EOF 法或是 KLE 法分別進行填補,則誤差會隨缺漏點個數變多而呈線性 增長,不過各方法在空間上之分佈皆甚為類似,在離雷達站較近處誤差較小;四 種方法中以最小平方法配合實向量EOF 法填補後之誤差相對較小,效果較佳。當 缺漏個數為資料點總數之 57%時,最小平方法填補後所得之誤差開始變得明顯,
當缺漏個數達資料點總數之 71%以上時,其變化趨勢會呈指數成長;而迭代法的 實驗結果則是幾乎在所有缺漏情況下,其誤差皆大於最小平方法的誤差。在缺漏 個數為資料點總數之71%以上的情形時,這些方法均不適合用來填補缺漏值。
關鍵字:高頻雷達、CODAR、經驗正交函數、Karhunen-Loève 展開法、填補缺值
Abstract
The Coastal Ocean Dynamics Applications Radar (CODAR) is a High-Frequency (HF) radar system with compact antennas. Recently, CODAR becomes widely used, for monitoring ocean surface currents remotely in nearly real time; its ability of large coverage on ocean surface and high resolution both in time and in space, makes CODAR an ideal tool for the operational oceanography, such as now-cast of currents and data-assimilation tasks. However, environmental effects such as interferences from obstacles, and/or from ionospheric disturbances, often hampers or weakens the strength of CODAR system, might deteriorate the data quality of CODAR, inducing incomplete datasets with missing data or holes in the designated observation region. The development of appropriate methodology for filling missing data is therefore a necessity deserving for further studies, in prior to the development of operational framework.
In the context, we have analyzed a nearly two-year long dataset of CODAR, which was observed by two HF radars located at Suao and Han-Ben, respectively, and provided by the Surface Current Observations at North-East Taiwan (SCONET) project;
basic statistics of the measurements show that the SCONET dataset satisfies the normality and weakly stationary condition. Further processing, by using both modal decomposition methods of Real-vector Empirical Orthogonal Function (REOF) and the Karhuren-Loeve Expansion (KLE), respectively, reveals that almost more than 96% of total variances of currents in the whole observation region can be interpreted by the first 20 modes of both methods. Therefore the first 20 modes of both methods are used for data reconstruction and data filling later.
An independent month-long time series of CODAR measurements is adopted for the data filling experiment, by which the incomplete dataset is generated by depleting
artificially assigned grid points in the original complete dataset, the measurements are therefore treated as the ground truth for later comparison. Monte Carlo simulation is used for the tests of data filling experiment when the missing points of data exceed 3 points. We have used both the EOF and the KLE methods, in accompany with the least square and the iteration procedures for the estimation of the amplitude of each modes, for the study of data filling experiment. Results show that the EOF method in accompany with the least square procedure is the best among the four methodologies, when the percentage of occurrence of the missing data is less than 57% of the whole dataset. However, all these four methods are not adequate for filling incomplete dataset, if the percentage of occurrence of the missing data exceeds 71%.
Key words: HF radar, CODAR, EOF, Karhunen-Loève expansion, data filling
目 錄
口試委員會審定書 ... i
誌謝 ... ii
摘要 ... iii
Abstract ... iv
目 錄 ... vi
圖目錄 ... viii
表目錄 ... xii
第一章 緒論 ... 1
1.1 高頻雷達測流簡介 ... 1
1.2 臺灣東北海域 CODAR 系統簡介 ... 4
1.3 缺漏值填補相關之文獻回顧 ... 8
1.4 論文內容架構 ... 9
第二章 資料統計特性 ... 11
2.1 資料來源 ... 11
2.2 時間域基本統計 ... 18
2.3 頻率域分析以及濾潮處理 ... 25
2.4 平穩性檢定 ... 30
第三章 模分析方法與分析結果 ... 35
3.1 經驗正交函數 ... 35
3.2 Karhunen-Loève 展開法 ... 47
3.3 以模重組法填補資料缺漏 ... 58
3.3.1 最小平方法 ... 58
3.3.2 迭代法 ... 60
第四章 模擬實驗與結果討論 ... 62
4.1 實驗設計 ... 62
4.2 窮舉法實驗:任意 1~3 個資料點缺漏 ... 66
4.3 蒙地卡羅法之實驗:多個資料點出現缺漏 ... 74
4.4. 綜合討論 ... 79
第五章 結論 ... 85
參考文獻 ... 87
圖目錄
圖1-1 (a)雷達波與波浪作用示意圖;(b)上:無海流作用時的雷達回波頻譜圖,下:
有海流作用時的頻譜圖(縱軸表示回波能量密度,橫軸為頻率) (引自 Barrick et al., 1977)。 ... 2 圖1-2 美國加州利用 CODAR 雷達並配合 AVHRR 衛星在 Monterey Bay 觀測到的 海灣內及外海的海流和海表溫(引自 Paduan and Rosenfeld, 1996)。 ... 3 圖1-3 臺灣周圍海底地形變化示意圖(海底地形資料來源:國科會海洋學門資料
庫)。 ... 4 圖1-4 臺灣四周海域表層海流即時觀測網各 CODAR 雷達站站位分佈示意圖(站
位地點與名稱係參考國研院臺灣海洋科技中心海洋資料庫之CODAR 資料庫)。
... 5 圖1-5 蘇澳及漢本二雷達站(紅色方塊)及合成海流資料點位(深橘色十字)示意圖
(海底地形資料來源:國科會海洋學門資料庫)。 ... 6 圖1-6 合成海流資料點位之資料缺漏數量(組數)統計(2011 年 4 月 14 日至 2012 年
12 月 31 日,每組為一小時)。 ... 7 圖2-1 2012 年 4 月 30 日 1200UTC 蘇澳站短時間互能譜資料檔(CSS File)中三支
天線回波的自能譜(auto-spectra)。橫軸為頻率(中央為 0,左側為負,右側為正),
縱軸為與接收天線的水平距離(km),顏色越亮表示回波訊號越強。 ... 12 圖2-2 蘇澳站 2012 年 4 月 30 日 1200UTC 徑向海流分佈。徑向海流速度即是指
觀測區內某點海流流速指向該點與雷達站連線之射線上的流速分量,也就是 海流接近或遠離雷達站的流速分量。 ... 12 圖2-3 2012 年 4 月 30 日 1200UTC 的合成海流流速分佈圖,係由蘇澳與漢本二站
徑向海流資料經中控站系統結合而成。 ... 13 圖2-4 合成海流資料點位出現資料缺漏情況之統計(2011 年 4 月 14 日至 2012 年
12 月 31 日):(a)資料缺漏之百分比分佈;(b)不完整資料中缺漏點數之百分比 統計直方圖;(c)長期統計資料缺漏百分比之平均日變化。 ... 15 圖2-5 資料點位編號之相對位置圖。 ... 21 圖2-6 觀測區內之平均流場及水深 20 公尺之長期平均流場(海底地形與長期平均 流場資料來源:國科會海洋學門資料庫)。 ... 21 圖2-7 CODAR 觀測海流資料之標準差分佈:(a)東西向流速;(b)南北向流速。
... 22 圖2-8 (a)東西向流速分佈圖;(b)南北向流速分佈圖;(c)東西向流速 P-P 圖;(d)
南北向流速P-P 圖。紅線表示常態分佈曲線。 ... 24 圖2-9 中央氣象局蘇澳站(121.87°E、24.59°N) 2011/8/16~2011/11/15 的(a)潮位觀
測資料及(b)各分潮振幅。 ... 25 圖2-10 (a)第 5 點位、(b)第 69 點位之流速動能譜圖。 ... 27 圖2-11 潮流能量空間分佈圖:(a)全日潮分量;(b)半日潮分量。 ... 28 圖2-12 觀測區內第 1 個資料點海流資料(2011 年 8 月 16 日至 11 月 15 日)經調和 分析後得到的各分潮潮流橢圓長軸係數。 ... 29 圖2-13 流速速率資料之平穩性檢定結果,彩色區塊表示通過檢定,表示資料具
有穩定性,而黑色區塊表示未通過檢定:(a)平均值之檢定結果;(b)標準差之 檢定結果。 ... 32 圖3-1 實向量 EOF 分析法重建流場之特徵值分佈圖與累積特徵值分佈圖:(a)全
部140 個模組;(b)前 20 個模組。 ... 40 圖3-2 實向量 EOF 法前 20 個模組重建流場與變動流場之變異數百分比空間分佈
圖。 ... 40 圖3-3 實向量 EOF 分析法求出之第 1 模組:(a)模向量空間分佈;(b)模振幅之時
間序列。 ... 42 圖3-4 實向量 EOF 分析法求出之第 2 模組:(a)模向量空間分佈;(b)模振幅之時
間序列。 ... 43 圖3-5 實向量 EOF 分析法求出之第 3 模組:(a)模向量空間分佈;(b)模振幅之時
間序列。 ... 44 圖3-6 原始觀測流場與實向量 EOF 法前 20 個模組重建流場之偏差均方根分佈圖。
... 46 圖3-7 KLE 分析法重建流場之振幅平方和分佈圖與累積變異數分佈圖:(a)全部
140 個模組;(b)前 20 個模組。 ... 51 圖3-8 KLE 法前 20 個模組重建流場與變動流場之變異數百分比空間分佈圖。 52 圖3-9 KLE 分析法求出之第 1 模組:(a)模向量空間分佈;(b)模振幅之時間序列。
... 54 圖3-10 KLE 分析法求出之第 2 模組:(a)模向量空間分佈;(b)模振幅之時間序列。
... 55 圖3-11 KLE 分析法求出之第 3 模組:(a)模向量空間分佈;(b)模振幅之時間序列。
... 56 圖3-12 原始觀測流場與 KLE 法前 20 個模組重建流場之偏差均方根空間分佈圖。
... 57 圖4-1 蒙地卡羅法前級試算結果:假設缺漏個數為 5,以最小平方法搭配實向量
EOF 法填補,橫軸為實驗次數,縱軸為填補值與真實值之誤差均方根平均值。
... 64 圖4-2 任意 1 個資料點缺漏情形之填補實驗結果:(a)以最小平方法配合實向量
EOF 法;(b)以最小平方法配合 KLE 法;(c)以迭代法配合實向量 EOF 法;(d) 以迭代法搭配KLE 法填補之誤差均方根空間分佈。 ... 70 圖4-3 任意 1 個資料點缺漏情形之填補實驗結果(已扣除 4.1 節所述第 1 種誤差):
(a)以最小平方法配合實向量 EOF 法;(b)以最小平方法配合 KLE 法;(c)以迭 代法配合實向量EOF 法;(d)以迭代法搭配 KLE 法填補之誤差均方根空間分
佈。 ... 73 圖4-4 任意 10 個資料點缺漏情形之填補實驗結果:(a)以最小平方法配合實向量
EOF 法;(b)以最小平方法配合 KLE 法;(c)以迭代法配合實向量 EOF 法;(d) 以迭代法搭配KLE 法填補之誤差均方根空間分佈。 ... 75 圖4-5 任意 10 個資料點缺漏情形之填補實驗結果(已扣除 4.1 節所述第 1 種誤差):
(a)以最小平方法配合實向量 EOF 法;(b)以最小平方法配合 KLE 法;(c)以迭 代法配合實向量EOF 法;(d)以迭代法搭配 KLE 法填補之誤差均方根空間分 佈。 ... 76 圖4-6 任意 35 個資料點缺漏情形之填補實驗結果:(a)以最小平方法配合實向量
EOF 法;(b)以最小平方法配合 KLE 法;(c)以迭代法配合實向量 EOF 法;(d) 以迭代法搭配KLE 法填補之誤差均方根空間分佈。 ... 77 圖4-7 任意 35 個資料點缺漏情形之填補實驗結果(已扣除 4.1 節所述第 1 種誤差):
(a)以最小平方法配合實向量 EOF 法;(b)以最小平方法配合 KLE 法;(c)以迭 代法配合實向量EOF 法;(d)以迭代法搭配 KLE 法填補之誤差均方根空間分 佈。 ... 78 圖4-8 以實向量 EOF 法、KLE 法與最小平方法、迭代法進行填補不同缺漏個數
實驗結果:(a)填補值與真實值之誤差均方根平均值;(b)填補值與 20 個模組重 建流場之誤差均方根平均值分佈。 ... 80 圖4-9 以實向量 EOF 法、KLE 法與最小平方法、迭代法進行填補不同缺漏個數
實驗結果:(a)填補值與真實值之誤差均方根標準差;(b)填補值與 20 個模組重 建流場之誤差均方根標準差分佈。 ... 81 圖4-10 實向量 EOF 法與 KLE 法誤差均方根平均值之比較。 ... 82 圖4-11 F 檢測:比較(a)最小平方法;(b)迭代法之結果。 ... 84
表目錄
表2-1 合成海流資料中發生缺漏情形(共 429 組)之統計。 ... 16
表2-2 長期統計資料缺漏百分比之平均日變化。 ... 17
表2-3 各點位之東西向流速(U)與南北向流速(V)的基本統計。 ... 19
表2-4 全部資料東西向流速(U)與南北向流速(V)的基本統計資料。 ... 23
表2-5 各點位之東西向流速(U)、南北向流速(V)及速率(speed)的連個數。 ... 33
表3-1 原始觀測流場與實向量 EOF 法前 20 個模組重建流場之偏差均方根數值。 ... 46
表3-2 原始觀測流場與 KLE 法前 20 個模組重建流場之偏差均方根數值。 ... 57
表4-1 70 個資料點不同缺漏個數情況之變化種數。 ... 63
表4-2 填補任意 1 個資料點為缺漏情形之實驗結果。 ... 68
表4-3 填補任意 1 個資料點為缺漏情形之實驗結果(扣除空間上之不均勻性)。 71 表4-4 顯著水準α=0.05時與各自由度相對應之信賴區間。 ... 79
表4-5 F 分配表中,顯著水準為α及自由度為 n 時的上、下限數值。 ... 83
第一章 緒論
1.1 高頻雷達測流簡介
隨著文明以及經濟發展,人類社會對於海上活動的需求以及頻度均大幅增加。
基於安全以及經濟效益等因素考量,海上活動諸如航運、遊憩、救難搜索或是汙 染追蹤等都需要參考海洋的資訊,特別是海面流場變化更是其中重要的一環。海 表面流場除了受背景海流場之作用外,還會受到波浪、風場、潮汐等外在因素的 影響,其變化十分複雜。海流觀測傳統上是利用漂流浮標(drifter),或是錨碇(mooring) 海流儀(如都卜勒流速儀 Acoustic Doppler Current Profiler, ADCP,或其他型式海流 儀)等方法,缺點是在時間和空間尺度上難以同時兼顧。近代發展出的海流遙測技 術,例如岸基高頻雷達(shore-based high-frequency radars),則具有觀測資料空間涵 蓋範圍大(尺度從數百公尺到數百公里)和時間解析度高(可以小至數分鐘)的優點,
目前雖然其量測精確度尚不能與傳統方法相比,但對作業化海流觀測而言這已是 一種可行的觀測工具(Barrick et al., 1977; Paduan and Rosenfeld, 1996)。
國外運用高頻(high-frequency, HF)雷達觀測海流之研究已發展了近 40 年之久 (Teague et al., 1997),其遙測原理如圖 1-1 所示:由雷達天線發射出的電磁波在傳 播過程中接觸到海面粗糙構造時會產生後向散射(backscattering),後向散射的回波 其性質與會影響到海面粗糙度特性和分佈的物理因子(如波浪、海流等)之間有密切 的關係,所以利用雷達回波的訊號可以反推出這些物理因子的分佈特性(Barrick et al., 1977; Paduan and Rosenfeld, 1996)。後向回波訊號主要是由布拉格散射(Bragg scattering)所造成,也就是當波長為 λ 的雷達波遇到海面上波長為
(
1,2,3,)
2 n= nλ
的水波時,雷達後向散射回波會發生建設性的疊加效應,因此回波信號就會較強(圖 1-1a),其中尤以n=1時為最強。當海面無流時,雷達回波頻譜圖上會在此
2 λ波長
水波所對應之頻率處出現能量尖峰(如圖 1-1b 上圖);但當海面有流時,由於海流
會對水波產生都卜勒頻移(Doppler shift),因此在雷達回波頻譜圖上能量尖峰頻率 亦會發生偏移(圖 1-1b 下圖),故從雷達回波的頻移可以反算出表面海流的大小。
不過單一雷達只能測到雷達天線徑向方向上的海流速度分量,需要二正交之分量 才能結合成為海流速度向量,是以至少要結合兩座以上的雷達才能得到共同涵蓋 範圍內的海流向量流場資訊(Barrick et al., 1977; Paduan and Rosenfeld, 1996)。
(a)
(b)
圖1-1 (a)雷達波與波浪作用示意圖;(b)上:無海流作用時的雷達回波頻譜圖,下:
有海流作用時的頻譜圖(縱軸表示回波能量密度,橫軸為頻率) (引自 Barrick et al., 1977)。
Emery et al. (2004)指出目前使用的高頻雷達測流系統大致可分成兩類:其一為 波束形成(beam forming)雷達,例如 Wellen radar (WERA; Gurgel et al., 1999);另一 種則是測向(direction finding)雷達,如 Coastal Ocean Dynamics Application Radar (CODAR)。其中 CODAR 系統因為體積小便於攜帶(Barrick et al., 1977),且較早商
業化(1984 年即已開始;參見 CODAR Ocean Sensors 網站,http://www.codar.com/),
目前使用較廣泛(Fang et al., 2011)。Paduan and Rosenfeld (1996)曾在美國加州的 Monterey Bay 地區利用 1994 年的 CODAR 雷達資料與 AVHRR 衛星資料所觀測到 海灣及外海的海表溫相互比對(如圖 1-2),其結果顯示 CODAR 測得的平均流場(黑 色箭頭)在灣口處呈現出一股往南的海流,而灣內的海表溫也同樣呈現往南方擴展 的狀況,這種分佈態勢相當符合當地歷史水文資料統計特性,因此他們認為高頻 雷達測流系統在長期觀測以及即時監測海流等作業方面具有很大的發展潛力。
圖1-2 美國加州利用 CODAR 雷達並配合 AVHRR 衛星在 Monterey Bay 觀測到的 海灣內及外海的海流和海表溫(引自 Paduan and Rosenfeld, 1996)。
1.2 臺灣東北海域 CODAR 系統簡介
臺灣位於東亞海運樞紐,週遭有幾道重要海流,如黑潮、黑潮支流、大陸沿 岸流、南海暖流等通過,另外鄰近海域海底地形崎嶇(圖 1-3),風場多變,更增加 了流況的複雜性(Jan et al., 2002)。由於運用 CODAR 雷達即時監測海流可以獲得空 間以及時間解析度均甚高的流況資料,我國海洋學界多年前即倡議使用此系統來 建立臺灣四周海域表層海流即時觀測網,所獲資料除可強化海洋資料庫外並可結 合數值模式發展預報技術,以提供各界在研究、防災、救難及海上休憩活動等方 面上的需求,目前觀測網建置作業已初步完成(環台各 CODAR 雷達站分佈如圖 1-4,
參見國研院臺灣海洋科技中心網站,http://www.tori.org.tw/)。
圖 1-3 臺灣周圍海底地形變化示意圖(海底地形資料來源:國科會海洋學門資料 庫)。
圖 1-4 臺灣四周海域表層海流即時觀測網各 CODAR 雷達站站位分佈示意圖(站 位地點與名稱係參考國研院臺灣海洋科技中心海洋資料庫之CODAR 資料庫)。
在圖1-4 所示之環臺 CODAR 觀測站中,蘇澳與漢本二站係由國立臺灣大學海 洋所、海軍軍官學校以及海軍大氣海洋局共同合作建置並維持,詳細情形可參見
「臺灣東北海域海流觀測」(Surface Currents Observation North East of Taiwan, SCONET)網頁(http://www-codar.oc.ntu.edu.tw)。二站量測區涵蓋範圍包括臺灣東北 部宜蘭海脊以及和平海盆交界處一帶海域,當地海底地形變化相當劇烈(圖 1-3),
亦是黑潮流經之處。二站中,蘇澳站(Suao)位於北方澳海邊的海岬上,經度 121.8825oE,緯度 24.6002oN,漢本站(Han-Ben)則位於和平溪北岸海邊,經度 121.7695oE,緯度 24.3301oN;二雷達站均使用長距型 CODAR 系統(Long-Range SeaSonde),各具有一支直立的發射天線以及三支相互垂直的接收天線,所用頻率 均為4.3~5.4 MHz,空間解析度則取決於所用之頻寬,目前蘇澳與漢本二站均為 8.16 km (Fang et al., 2011)。
雖然長距型 CODAR 系統在最佳情況下理論上的量測距離可達 220 km,但 Fang et al. (2011)發現蘇澳、漢本二站均會受到與那國島(Yonaguni,與二站距離均 為 110 km 左右)以及電離層干擾影響,尤其是後者在本地時間晚上 6 點到凌晨 2 點之間其效應最為嚴重。此外,雷達天線場型也會受到周遭環境的影響,若是雷 達接收天線設置的位置不佳,例如附近有較高的樹木或是建築物等之類的高大物 體,都會對天線場型產生不良效應(Kohut and Glenn, 2003; Ramp et al., 2008),有可 能導致系統在估算流速時發生誤差(Fang et al., 2011)。Fang et al. (2011)經由海流比 對、校正,並依資料統計結果以及考慮兩站徑向交角大小等效應後,將蘇澳和漢 本二站CODAR 系統合成海流之產出結果限定在如圖 1-5 所示之 70 個站點內。在 此觀測範圍內約有 26 個觀測點位於宜蘭海脊上,另外 44 個觀測點則位於較深處 的和平海盆。為了簡化,以下將蘇澳、漢本二站 CODAR 系統之觀測海域簡稱為 SCONET 海域。
圖 1-5 蘇澳及漢本二雷達站(紅色方塊)及合成海流資料點位(深橘色十字)示意圖 (海底地形資料來源:國科會海洋學門資料庫)。
蘇澳、漢本二站自2011 年 4 月正式開始連線觀測以來,迄 2012 年 12 月所收 到的海流資料中,不免會因雷達站斷電、儀器故障,或是電離層、環境雜訊干擾 以及雷達波被遮蔽等因素影響,致使回波訊號太弱或接收不到雷達回波訊號,因 而造成此70 個觀測點合成海流資料中不時會發生有些點位無資料的情形,此期間 內各觀測點之資料缺漏數量統計如圖1-6 所示。整體而言,此期間內每個觀測點缺 漏資料之組數約介於 600~900 組之間,其數量與觀測點和雷達站之距離似乎呈正 比關係,距離雷達站越遠處之資料缺漏情形則越嚴重(圖 1-6)。如何填補觀測資料 缺漏值是發展高頻雷達海流現報作業化系統過程中一個實際的問題,因此也構成 了一個有意義的研究課題。
圖1-6 合成海流資料點位之資料缺漏數量(組數)統計(2011 年 4 月 14 日至 2012 年 12 月 31 日,每組為一小時)。
1.3 缺漏值填補相關之文獻回顧
「如何填補缺漏資料(missing data)」是不同領域的資料分析,特別是對於遙測 資料或是一些歷史水文資料而言,這是一個很常見的問題(Beckers and Rixen, 2003;
Kondrashov and Ghil, 2006; Oliveira and Gomes, 2010)。造成觀測記錄中出現缺漏值 的可能原因包括有儀器故障、有遮蔽物無法獲得資料(例如被雲覆蓋)、受到雜訊干 擾等,特別是分析時空變化資料時,有缺漏值往往使傳統方法無法順利處理因而 形成問題(Beckers and Rixen, 2003; Kondrashov and Ghil, 2006)。Han and Kamber (2006)曾提出資料清理(data cleansing)的概念,也就是在分析數據時要先進行包括 填補缺漏值、消除雜訊及校正數據等前級處理過程,他們同時也提出了缺值填補 方法的大方向,例如可以刪去有缺漏值的整筆資料、以平均值來填補,或者是使 用以最有效的分析方法所得到之可能性最大的值來填補,後者亦是當下常用的補 值方法。
填補缺值可以用參數相依模型(parameter-dependent models)作為模型基底 (model-based)來推估缺漏值,例如使用最佳內差法(optimal interpolation)、經驗正交 函數(empirical orthogonal function, EOF)、卡爾曼濾波(Kalman filtering) (Reynolds and Smith, 1994; Smith et al., 1996; Kaplan et al., 1997; Kondrashov and Ghil, 2006),
或是使用客觀分析法(objective analysis) (Bretherton et al., 1976)等方法。在這些方法 中,經驗正交函數在物理海洋學以及大氣方面應用甚廣,例如Boyd et al. (1994)以 EOF 法分析 Sargasso Sea 的 CTD 溫度資料,他們分別測試在不同深度有缺值時,
使用前5 個模組(mode)所合成的溫度資料與實際 CTD 量測之溫度資料間的誤差情 形,結果發現結合水溫氣候值(climatological values)和 EOF 法合成的資料與 CTD 直接量測值間的誤差會比僅用水溫氣候值來推估之誤差要小;又如 Beckers and Rixen (2003) 選擇義大利半島和巴爾幹半島間的亞得里亞海(Adriatic Sea)中一塊 區域做資料填補研究,此區域的AVHRR 影像常受雲層及雜訊的影響,資料品質並 不好,而利用EOF 法則能夠對原始資料做過濾(filter)和內差補值(interpolation)處理;
另外Alvera-Azcárate et al. (2005)也是使用亞得里亞海的 AVHRR 資料進行研究,運 用在不同的雲覆蓋量,也就是有不同數量的資料缺值情況下,以EOF 法填補之結 果與原始資料相比較,發現填補結果之誤差雖會隨雲量變多而增加,但並不至於 增加太多。
此外, Karhunen-Loève 展開法(Karhunen-Loève expansion, KLE) (Preisendorfer, 1988)是一種本質上與 EOF 法非常類似的方法,在流體力學以及工程領域亦有學者 將其應用於填補缺漏值方面:例如Bui-Thanh et al. (2004)使用在機翼表面量測到的 壓力資料,經由KLE 法填補缺漏資料之後,發現 KLE 法之重建資料與計算流體力 學模式算出的結果相當吻合;而Snihir et al. (2006)則是使用 KLE 法預測由不同電 流將電池充電後所獲得的開路電壓(open-circuit voltage)曲線,其實驗結果認為此方 法在電池管理系統(battery management systems)的降維技術(dimension reduction techniques)上有很大的潛力。
1.4 論文內容架構
由前節的文獻回顧可知,利用參數相依模型為基底之統計方法可以有效過濾、
填補缺漏資料,這類方法亦有學者應用於高頻雷達資料上(例如 Emery et al., 2006)。
將 CODAR 即時觀測資料裡的缺漏值適當填補後有助於海流資料分析,特別是資 料同化(data assimilation)作業。Kohut et al. (1999)在十餘年前即倡議將 CODAR 的 即時海流觀測資料進行資料同化後加入數值模式,他們認為如此當可有效提升海 洋環境變化預報的正確性。由於經驗正交函數是物理海洋研究上常被使用的方法 (Beckers and Rixen, 2003),此法能有效地將資料特性用簡潔(compact)的形式呈現 (Kaihatu et al., 1998),故有不少學者將其應用在填補缺漏值方面(Boyd et al., 1994;
Beckers and Rixen, 2003; Alvera-Azcárate et al., 2005);另方面,與 EOF 法類似之 Karhunen-Loève 展開法則在工程計算或地球物理領域上使用甚多(Tiampo et al.,
2004; Dong, et al., 2006; Snihir et al., 2006),但未見用於海洋資料處理方面,有鑑於 此,本文將使用此方法處理CODAR 觀測海流資料,並進一步與 EOF 法之結果相 互比較。在以下的章節內,我們將分別探討以EOF 法以及 KLE 法對 CODAR 海流 資料進行模分解(modal decomposition)之效果,並比較以二法重組資料之優缺點,
另外亦在不同缺漏值個數的情況之下,檢驗上述二法之填補結果成效與誤差大小,
並推測使用此二法時缺漏資料點個數的上限。
本文架構如下,第一章簡介 CODAR 系統之運作原理以及臺灣東北海域 CODAR 觀測作業現況,另外也對前人如何應用統計方法填補缺漏值進行簡要的文 獻回顧;第二章介紹 CODAR 雷達觀測海流資料之處理方法、即時資料傳輸流程 以及資料統計特性;第三章敘述 CODAR 資料之模分解分析方法與分析結果,比 較原始觀測流場與經過模分解後再重建的流場資料之差異度;第四章以真實觀測 數據為根據,假設資料出現不同數量的缺漏情況,統計各不同填補方法所得結果 與真實值之差異情形;第五章則為結論,比較各方法之優劣與限制後總結全文。
第二章 資料統計特性
在這一章內我們將簡介SCONET 海域 CODAR 系統處理海流觀測資料的方法 與傳輸流程,並說明資料的統計特性,包含平均流場、流速分量機率分佈情形等 基本統計。另外,亦對流速資料之常態性以及平穩性進行了統計檢定。
2.1 資料來源
如前章所述,蘇澳與漢本二站均是使用長距型的CODAR高頻雷達測流系統,
根據Fang et al. (2011),二站CODAR系統之作業方式為在接收到回波信號後即先將 其快速傅立葉轉換(Fast Fourier Transform, FFT)轉為能譜(spectrum)形式暫存,接著 約經每15分鐘累積後自動產生一組互能譜資料檔(Cross Spectrum File, CSQ File),
以及約每30分鐘系統再將前後之CSQ File進行平均,產生一組稱為短時間互能譜的 資料檔(Short Time Cross Spectrum File, CSS File) (如圖2-1)。為了有效提升訊噪比,
CODAR系統又會以經驗時間長度(一般多使用3小時,但蘇澳與漢本二站係設定使 用4小時,Fang et al. 2011)對所收集的CSS File再進行一次平均(Ramp et al., 2008),
並以此作為估算海流的基礎,二站CODAR系統每小時產出的流速資料即是根據這 種平均後的能譜所估算,會以雷達站為中心按極座標方式輸出每間隔1度方位角以 及徑向上每8.16 km一點之徑向海流資料檔(Radial Velocity File, RUV File) (如圖 2-2)。最後再將不同雷達站的同步徑向海流觀測資料回傳到中控站系統並由後者做 整合,產出合成表面海流資料檔(Total Velocity File, TUV File) (如圖2-3),徑向海流 資料檔與合成表面海流資料檔皆為每小時輸出一次。
圖2-1 2012 年 4 月 30 日 1200UTC 蘇澳站短時間互能譜資料檔(CSS File)中三支 天線回波的自能譜(auto-spectra)。橫軸為頻率(中央為 0,左側為負,右側為正),
縱軸為與接收天線的水平距離(km),顏色越亮表示回波訊號越強。
圖2-2 蘇澳站 2012 年 4 月 30 日 1200UTC 徑向海流分佈。徑向海流速度即是指 觀測區內某點海流流速指向該點與雷達站連線之射線上的流速分量,也就是海流 接近或遠離雷達站的流速分量。
圖2-3 2012 年 4 月 30 日 1200UTC 的合成海流流速分佈圖,係由蘇澳與漢本二站 徑向海流資料經中控站系統結合而成。
Fang et al. (2011)將國科會海洋學門資料庫累積之歷史海流資料與蘇澳、漢本 二 CODAR 長期觀測資料相比對,並參考二站互能譜資料中所記錄之資料品質參 數統計,最後設定了資料品質相對較好且穩定的 70 個點(如圖 1-5 所示)作為中控 站系統產生合成海流場之輸出範圍;我們在以下各章節內所使用的資料均為這類 產品。本文分析所用資料是從2011 年 4 月 14 日 2300UTC 開始,到 2012 年 12 月 31 日 2200UTC 為止,每小時一組(如圖 2-3 即為其中一例),共計 15048 組。其中 前20 個月(即由 2011 年 4 月 14 日 2300UTC 至 2012 年 11 月 30 日 2300UTC 止)共 14305 組資料將用於 EOF 以及 KLE 模分析(詳見第三章),而最後一個月(2012 年 12 月 1 日 0000UTC 至 2012 年 12 月 31 日 2200UTC 止)共 743 組資料則當作獨立 樣本,作為填補缺值實驗之根據。
如1.2 節所述,蘇澳、漢本二站有時會因各種因素影響會造成合成海流資料中 出現某些點位無資料的情形。在 15048 組觀測資料中,各觀測點出現缺漏值之數 量統計如圖 1-6,換算成百分比出現率則如圖 2-4(a)所示,可看出各點資料之最大 缺漏率幾乎都在5%以下。此外,資料中完全沒有任何缺漏值的組數為 14014,約 佔全部資料的93.13%,有缺值情況的組數為 1034,共佔 6.87%;後者中,因儀器 故障完全沒有任何觀測資料的組數則為605,佔全部資料的 4.02%;而其餘的 429 組,資料缺漏量由1 至 69 點,其百分比出現率(以總樣本數 429 計)分佈則如圖 2-4(b) 所示,其中以每組缺漏一點的情形數量最多,共 126 次,佔有缺漏值資料組的 29.37%,缺二點佔 11.42%,缺三點則為 9.79%,顯示觀測值發生缺漏情形時,約 有一半左右的情況都是在三點以內,至於其他缺漏數量情況詳見表 2-1。另方面,
將15048 組資料中所有有缺漏值的組數(共 1034 組)按其每日出現之小時數另行統 計,則每小時總共約發生30~60 次缺漏,詳細數據見表 2-2。圖 2-4(c)為每日資料 缺漏量百分比之逐時變化分佈,缺漏量之每小時出現率(以樣本數 1034 計)都在 6%
以下,但在每日UTC 時間 10 時至 20 時(約本地時間晚上 6 點到凌晨 4 點)則比較 容易發生資料缺漏情形,此時段內缺漏事件之出現率約為4~6%,與 Fang et al. (2011) 所述SCONET 海域 CODAR 系統受電離層干擾影響時間大致吻合。
(b)
(c)
圖2-4 合成海流資料點位出現資料缺漏情況之統計(2011 年 4 月 14 日至 2012 年 12 月 31 日):(a)資料缺漏之百分比分佈;(b)不完整資料中缺漏點數之百分比統計 直方圖;(c)長期統計資料缺漏百分比之平均日變化。
(a)
表2-1 合成海流資料中發生缺漏情形(共 429 組)之統計。
缺漏值個數 組 % 缺漏值個數 組 % 缺漏值個數 組 %
1 126 29.37 24 1 0.23 47 1 0.23 2 49 11.42 25 2 0.47 48 2 0.47 3 42 9.79 26 1 0.23 49 3 0.70 4 21 4.90 27 1 0.23 50 1 0.23 5 10 2.33 28 3 0.70 51 0 0.00 6 10 2.33 29 0 0.00 52 2 0.47 7 10 2.33 30 1 0.23 53 2 0.47 8 8 1.86 31 0 0.00 54 3 0.70 9 12 2.80 32 3 0.70 55 4 0.93 10 9 2.10 33 0 0.00 56 8 1.86 11 7 1.63 34 0 0.00 57 6 1.40 12 6 1.40 35 2 0.47 58 1 0.23 13 5 1.17 36 1 0.23 59 4 0.93 14 8 1.86 37 1 0.23 60 1 0.23 15 5 1.17 38 1 0.23 61 0 0.00 16 4 0.93 39 1 0.23 62 2 0.47 17 4 0.93 40 2 0.47 63 6 1.40 18 3 0.70 41 3 0.70 64 0 0.00 19 4 0.93 42 2 0.47 65 1 0.23 20 3 0.70 43 0 0.00 66 0 0.00 21 1 0.23 44 2 0.47 67 0 0.00 22 2 0.47 45 0 0.00 68 1 0.23 23 4 0.93 46 1 0.23 69 0 0.00
表2-2 長期統計資料缺漏百分比之平均日變化。
UTC Hour No Missing Missing Missing (%)
0 592 35 3.38
1 587 40 3.87
2 589 38 3.68
3 591 36 3.48
4 590 37 3.58
5 594 33 3.19
6 591 36 3.48
7 592 35 3.38
8 591 36 3.48
9 588 39 3.77
10 582 45 4.35 11 579 48 4.64 12 584 43 4.16 13 578 49 4.74 14 573 54 5.22 15 575 52 5.03 16 571 56 5.42 17 570 57 5.51 18 575 52 5.03 19 576 51 4.93 20 578 49 4.74 21 588 39 3.77 22 591 36 3.48 23 589 38 3.68 Total 14014 1034 100.00
2.2 時間域基本統計
為了便於分辨與說明起見,首先我們將圖1-5 所示之 SCONET 海域合成海流 資料點位按照由下而上、由左而右之順序編碼命名(由 1 至 70),其分佈如圖 2-5。
接著將2011 年 4 月 14 日至 2012 年 12 月 31 日期間此 70 個資料點之東西向流速(U) 與南北向流速(V)分別求其平均值(mean)、標準差(standard deviation, SD),以及中 位數(median),統計結果詳如表 2-3。
SCONET 海域內平均之觀測海流分佈情形如圖 2-6 所示,圖中深橘色箭頭即 為平均海流流速向量,各點中最大流速為84.87 cm/s,最小則為 28.79 cm/s,流向 為北至東北,和我們所認知的黑潮流況相符。觀測區內平均海流流速較小處有二,
其一是位於觀測區東南側,大約在122.25oE~122.5oE 之間以及 24.25oN 以南的海域,
另一處則是在宜蘭海脊上接近蘇澳的海域。圖2-6 中以黑色箭頭表示的是海洋資料 庫發佈的船載ADCP 水深 20 公尺長期平均流場(1991 年~2008 年),後者之海流流 速範圍在 4.52~112.41 cm/s 之間,流向為東北,海流流速較小處亦位於 SCONET 海域的東南部與西北部。比較兩種不同來源的長期平均流場,可看出二者之間的 差異在觀測區中部二者流速向量的偏差較小,且西部又較東部偏差小。二者之差 異可能和船載ADCP 與 CODAR 觀測資料之水深不同(前者水深約 20 公尺,後者 約1 公尺)有關。
計算出各資料點的平均海流流速後,再根據平均值計算各點位之海流分量標 準差,可反映海流流速變動之幅度大小。東西向流速標準差(圖 2-7a)介於 17.53~
38.62 cm/s 之間,此值在 SCONET 海域西北方較小、西南方較大,東部則大多介 於20~25 cm/s 之間;南北向流速的標準差分佈(圖 2-7b)則在 21.59~39.63 cm/s 之間,
較小值出現在SCONET 海域中南部,而北部的標準差大致上都較大。標準差與海 流流速變動幅度成正比,故根據標準差的分佈可看出在SCONET 海域西南部東西 向流速變動較大,而在SCONET 海域東北部則南北向流速變動較大。
表2-3 各點位之東西向流速(U)與南北向流速(V)的基本統計。
U Mean V Mean U SD V SD U Median V Median 1 12.00 79.81 33.75 31.53 13.30 81.67 2 13.94 75.82 27.32 31.39 14.75 78.24 3 19.57 68.34 21.70 31.41 20.03 71.68 4 17.56 50.32 19.45 30.57 17.26 53.41 5 20.50 36.27 18.22 31.36 19.99 37.48 6 29.03 29.53 17.57 30.98 29.74 28.84 7 33.86 24.26 19.16 29.53 36.26 23.70 8 16.62 78.57 38.62 29.76 18.05 80.30 9 12.10 75.96 30.20 29.65 12.78 77.00 10 14.78 72.80 25.95 29.19 14.82 73.65 11 18.61 65.24 24.23 28.47 18.91 66.65 12 23.19 53.89 20.31 30.09 23.32 55.75 13 26.51 36.21 17.53 30.84 26.35 36.22 14 35.21 31.23 17.84 30.53 35.82 30.08 15 37.21 27.80 18.99 29.93 40.00 27.23 16 34.50 24.48 20.44 28.59 38.06 26.17 17 2.50 64.09 34.73 26.26 3.77 64.99 18 3.47 63.15 26.50 25.78 3.44 63.86 19 9.45 61.66 23.65 25.36 8.76 61.96 20 16.82 58.02 22.70 24.72 15.68 57.99 21 25.11 51.55 21.76 27.58 25.12 50.48 22 34.53 42.75 18.85 28.79 34.09 41.30 23 41.67 34.90 17.96 30.18 42.41 32.60 24 41.25 32.84 18.74 29.63 43.77 32.44 25 40.49 27.51 20.80 31.67 43.54 28.79 26 -3.81 50.82 27.13 22.42 -3.62 50.93 27 0.22 50.74 21.57 22.61 -0.43 51.20 28 6.90 51.04 21.33 23.31 5.75 51.34 29 16.97 53.47 20.91 23.42 15.38 53.40 30 26.67 50.94 21.09 26.43 25.49 49.83 31 39.80 51.53 19.92 29.79 39.06 49.78 32 45.55 45.96 18.35 28.78 45.56 44.72 33 48.39 40.98 19.72 30.78 49.99 40.24 34 51.85 21.23 23.20 32.62 53.88 20.66
(cm/s) D.P
35 -13.10 32.50 23.02 21.59 -13.70 32.27 36 -3.52 37.79 20.19 22.75 -4.83 37.69 37 5.30 41.27 21.45 24.47 3.72 41.00 38 18.06 49.81 21.43 24.91 16.55 48.80 39 30.05 56.82 21.58 26.58 28.58 55.99 40 41.62 60.81 21.50 29.26 40.63 59.79 41 51.19 63.73 21.06 29.83 50.58 62.71 42 52.86 37.83 21.11 31.79 52.54 36.68 43 59.17 28.17 24.19 35.89 59.51 28.52 44 -11.01 28.63 21.14 22.69 -12.31 28.22 45 0.33 35.57 21.03 24.28 -1.46 35.02 46 9.09 42.28 22.39 26.82 6.95 40.50 47 21.49 54.90 22.38 27.46 19.68 53.25 48 32.33 63.51 22.48 28.07 30.83 62.43 49 41.98 66.65 22.83 30.69 40.73 65.79 50 52.02 67.06 23.07 30.89 50.87 65.34 51 57.22 43.65 23.31 33.48 55.87 42.95 52 -5.40 28.28 21.94 24.22 -6.86 26.42 53 6.97 40.82 23.09 27.23 5.14 38.58 54 15.16 50.25 23.40 28.38 13.24 48.20 55 24.71 60.51 23.05 28.94 23.19 59.22 56 32.24 63.92 22.83 30.36 30.67 62.44 57 39.48 58.70 23.17 31.38 37.89 56.84 58 48.39 55.69 23.80 33.10 46.65 53.82 59 56.75 53.37 24.63 34.57 54.72 52.51 60 13.29 49.49 25.03 30.38 11.33 47.09 61 19.09 56.12 24.57 30.93 17.15 53.77 62 25.53 61.20 23.62 31.47 24.13 59.64 63 30.49 56.46 22.78 30.66 29.15 55.45 64 36.45 47.42 23.25 32.78 34.77 45.84 65 44.64 43.01 23.65 36.52 42.56 41.26 66 20.42 56.33 25.20 33.73 18.50 53.68 67 25.36 55.48 23.55 31.61 24.07 53.80 68 28.31 44.89 22.45 31.12 27.26 44.49 69 33.55 35.81 22.99 35.47 32.24 35.30 70 41.01 34.54 23.10 39.63 38.93 32.29
圖2-5 資料點位編號之相對位置圖。
圖2-6 觀測區內之平均流場及水深 20 公尺之長期平均流場(海底地形與長期平均 流場資料來源:國科會海洋學門資料庫)。
(a)
(b)
圖2-7 CODAR 觀測海流資料之標準差分佈:(a)東西向流速;(b)南北向流速。
為了觀察CODAR 流速觀測資料之機率分配情形,我們將 2011 年 4 月 14 日 至2012 年 12 月 31 日所有東西向流速與南北向流速資料統計分析,分別計算出各 自的平均值、標準誤(standard error, SE)及中位數(如表 2-4)。另外,再根據表中之 平均值以及標準誤,分別將東西向流速與南北向流速標準化(normalization),以Z表 示( s
Z = X −μ ,X 為流速資料,μ 、 s 則分別是資料的平均值與標準誤),並繪出 機率分佈直方圖(圖 2-8a 與 b),圖中直方圖是將標準化後的流速資料Z 等分成 50 個區間後再繪出,每個區間大小(ΔZ)分別為 0.26(圖 2-8a)與 0.24(圖 2-8b);紅色曲 線為流速資料的平均值與標準誤所構成的常態分配曲線,也就是說,圖2-8(a)與圖 2-8(b)的紅色曲線實際上分別是平均值為 25.65 cm/s、49.23 cm/s 以及標準誤為 28.86 cm/s、32.80 cm/s 的常態分配曲線。將此曲線與直方圖互相對照,可以從圖中看出 東西向與南北向流速的分佈都相當近似於常態分配。此外,圖2-8(c)與(d)則分別是 東西向和南北向流速的P-P 圖(probability- probability plot),藍圈是標準化後流速資 料Z與流速資料之累積分佈函數(cumulative distribution function, CDF)的分佈圖,
若資料為常態分配,則藍圈的分佈情形應與紅色對角線吻合,由 P-P 圖中二者幾 乎完全重疊可證明資料屬於常態分配。此外,將各點位之東西向、南北向流速同 樣分別經標準化後以Z表示並觀察其分佈情形,70 個點位的資料皆為常態分配(因 篇幅關係,此部份未列入論文)。
表2-4 全部資料東西向流速(U)與南北向流速(V)的基本統計資料。
Mean (cm/s) SE (cm/s) Median (cm/s) U 25.65 28.86 26.39 V 49.23 32.80 48.74
(a) (b)
(c) (d)
圖 2-8 (a)東西向流速分佈圖;(b)南北向流速分佈圖;(c)東西向流速 P-P 圖;(d) 南北向流速P-P 圖。紅線表示常態分佈曲線。
2.3 頻率域分析以及濾潮處理
潮汐與潮流是近岸海域裡重要的物理海洋現象;在淺水域裡潮流強度往往比 背景海流還要強勁。由劉 (1999)知,位於臺灣東北方的 SCONET 海域潮汐主要是 受混合潮型所主控,圖2-9(a)為 2011 年 8 月 16 日至 11 月 15 日蘇澳港的逐時潮位 曲線,可明顯看出本地潮汐14 天週期的大小潮變化以及小潮時潮型以全日潮為主 的特性。將這段潮位資料以T-tide (Pawlowicz et al., 2002)進行調和分析,圖 2-9(b) 即為主要一些分潮振幅分佈情形,根據Defant (1961)的潮型指標(Form Number),
2 2
1 1
S M
O F K
+
= + ,可算出蘇澳港的潮型指標F =0.52,因此可知鄰近蘇澳的SCONET
海域潮汐應當係屬以半日潮為主的混合潮型。
(a)
(b)
圖2-9 中央氣象局蘇澳站(121.87°E、24.59°N) 2011/8/16~2011/11/15 的(a)潮位觀 測資料及(b)各分潮振幅。
接下來進一步以能譜分析觀察各點位流速資料在頻率域之特性。首先將每個 點位的15048 筆東西向流速以及南北向流速資料皆分成 4 段(每段 4096 點),分別 扣除各段平均流速分量後,經快速傅立葉轉換(FFT)得出四組東西向和南北向流速 之能譜密度(power spectral density),然後將四組平均並在頻率域中每三個頻率相加 平均及重新取樣,如此即得到自由度(degree of freedom)為 24 的動能譜(kinetic energy spectrum),圖 2-10(a)為第 5 點位(相關位置參照圖 2-5)之動能譜,可以看出 頻率 f 越低,能量越大,此外在週期約為12 小時及 24 小時(即 2 cpd 和 1 cpd)處有 明顯的峰值,分別代表全日潮與半日潮的成份,另外在 0.05~1 cpd 之間,動能譜 有正比於 f−1.1的趨勢,不過在所有70 個點位中,這種特性之能譜僅出現在觀測區 西部的資料點上,在觀測區東部(以圖 2-10b 為例)及南部,海流能譜中全日潮與半 日潮分量之能量尖峰較不明顯,而在0.05~1 cpd 之間其能譜能階則約正比於 f−0.9。 此外,將各資料點動能譜中與全日潮、半日潮頻率相對應的能量挑出並繪成水平 分佈圖(圖 2-11),可見全日潮潮流在觀測區東北部、西北部的能量較強(圖 2-11a),
半日潮潮流則是在北部的能量較大(圖 2-11b),另方面半日潮潮流的能量在有些資 料點上會比全日潮潮流還大,最大的大約可達全日潮潮流能量的兩倍以上,但在 另外一些資料點上則又呈現相反的情形,這種能量分佈不均勻的特性可能和當地 的內潮(internal tide)有關(Jan et al., 2012)。由動能譜可知,CODAR 觀測資料中低頻 波動的能量較潮流為大。
若將各個資料點所觀測到的海流資料進行調和分析,可得出各資料點的分潮 潮流橢圓長軸係數,圖 2-12 即為一例,該圖顯示第 1 個資料點位由 2011 年 8 月 16 日至 11 月 15 日資料所算出之潮流橢圓長軸係數,在全日潮及半日潮(O1、K1
與M2、S2)以及低頻成份處潮流振幅均較大,全日潮及半日潮潮流橢圓長軸係數約 在 10 cm/s 左右,至於其他點位之潮流橢圓長軸係數則多小於 10 cm/s。與表 2-3 相比可知本地潮流流速頗小於平均流速,當然其中也有部分原因可能係緣於2.1 節
所述,為了提升訊雜比(Ramp et al., 2008),我們所用之徑向海流資料檔是經由 4 小 時平均處理後的CSS File 計算得出,因此合成海流資料中週期低於 4 小時的訊號 會被相當程度的抑制,而潮流成份也會比實際情形為弱(Fang et al., 2011)。為了去 除潮流的影響,我們先將前20 個月的原始觀測流場利用調和分析法求出各資料點 位之潮流調和常數,再追算合成出該時間內各點位之潮流,然後將原始觀測資料 扣除潮流,如此處理後即得出經過濾潮的海流資料,後者將使用於下一節之後的 各項分析。
(a)
(b)
圖2-10 (a)第 5 點位、(b)第 69 點位之流速動能譜圖。
全日潮
半日潮
全日潮 半日潮
(a)
(b)
圖2-11 潮流能量空間分佈圖:(a)全日潮分量;(b)半日潮分量。
圖2-12 觀測區內第 1 個資料點海流資料(2011 年 8 月 16 日至 11 月 15 日)經調和 分析後得到的各分潮潮流橢圓長軸係數。
2.4 平穩性檢定
對我們在第四章將要探討的填補資料缺漏實驗而言,由於填補方式是藉助於 實向量EOF法或KLE法(詳見第三章)得到的基底模組,但這些基底模組是根據「過 去」歷史資料之統計求出,因此如果統計結果是非穩定的(nonstationary),那麼後 續的補缺過程在邏輯上便會喪失立論基礎,是以確定資料是否具備統計穩定性便 是一項首要的工作。根據Bendat and Piersol (1971),如果一個時序資料其資料總體 平均(ensemble mean)以及共變異數(covariance)不隨起始時間的改變而變,而僅與各 資料點之時間差有關,那麼這組時序資料即滿足弱平穩(weakly stationary)的假定。
Bendat and Piersol (1971)指出連檢定(run test)是檢定資料平穩性(stationarity)的 一種作法。連檢定是一種無母數統計分析法(non-parametric statistics),或稱不受分 配限制的統計法(distribution-free),按其作法我們的檢定步驟如下:(1)首先定義虛 無假設(null hypothesis) H0為:資料是平穩的,而其對立假設(alternative hypothesis) H1則為:資料不是平穩的;(2)將前述經過濾除潮流後的14305組(約20個月)的 CODAR海流觀測資料等分為16段;(3)求出每一段之平均值與標準差,如此就得出 16個平均值與標準差的時序;(4)目測觀察此16個平均值與標準差之時序,如果沒 有明顯的趨勢(trend)變化,那麼H0就很有可能成立;(5)求出二序列數據的中位數,
以及二中位數與二序列資料之交點,根據交點可將序列分成R個區間,也就是R個 連,再由連分配(run distribution)表(Bendat and Piersol, 1971)檢定當設定之顯著水準 (significance level)為α時,如果H0成立,那麼連的個數R應介於rn;1−α2與rn;α2之間,
反之若R落在rn;1−α2與rn;α2的範圍外,則H0不成立;上述中rn;1−α2以及rn;α2為連分 配表給定之上、下限數值。
如果H0成立,根據α =5%之顯著水準以及n= N 2=16 2=8查表後可得出連 個數R應當介於r8;1−α2 =r8;1−0.052 =r8;0.975 =4和r8;α2 =r8;0.052 =r8;0.025 =13之間。表2-5 即為針對SCONET海域內70個資料點之連檢定分析結果,表中數字分別為各資料點
東西向流速(U)、南北向流速(V)以及流速速率(speed)之平均值與標準差的連個數,
表中第一行(D.P.)係指資料點位之編號(對應位置參見圖2-5),第二至第四行是以平 均值所計算出的連個數,最後三行則是用標準差所求出的連個數。圖2-13(a)與(b) 則分別是以各資料點之流速速率平均值和標準差進行連檢定的結果,圖中彩色區 域表示通過檢定(即各資料點的連個數介於4~13之間),而黑色區域則表示不通過檢 定,從二圖可知全部資料點的資料都通過檢定,即所有資料皆滿足弱平穩之假定,
應當可以向外延伸用於觀測時段之外。
(a)
(b)
圖 2-13 流速速率資料之平穩性檢定結果,彩色區塊表示通過檢定,表示資料具 有穩定性,而黑色區塊表示未通過檢定:(a)平均值之檢定結果;(b)標準差之檢定 結果。
表2-5 各點位之東西向流速(U)、南北向流速(V)及速率(speed)的連個數。
D.P. Umean Vmean Speedmean USD VSD SpeedSD
1 9 8 8 8 8 5 2 9 10 9 8 8 5 3 9 10 9 5 8 5 4 7 7 9 6 9 9 5 9 7 8 4 7 8 6 11 11 8 8 7 10 7 9 9 8 8 10 11 8 5 10 8 8 8 5 9 9 10 10 4 8 7 10 9 10 9 7 9 7 11 7 9 9 7 9 7 12 9 9 8 6 9 5 13 11 7 8 4 9 8 14 11 7 8 6 11 8 15 13 11 8 9 10 12 16 13 11 4 7 10 8 17 4 10 6 6 8 5 18 11 10 6 5 5 7 19 9 10 7 7 9 7 20 7 9 7 4 7 9 21 7 7 9 6 9 7 22 7 7 8 6 9 7 23 13 7 8 7 9 7 24 13 11 10 9 11 7 25 13 9 6 7 9 9 26 6 12 8 8 11 9 27 9 10 5 5 11 9 28 9 11 5 5 11 9 29 9 11 5 4 9 9 30 7 7 7 4 9 7 31 9 7 9 6 9 9 32 10 5 9 6 9 7 33 10 9 4 7 9 9 34 11 7 4 7 9 9
35 7 9 7 5 9 9 36 9 11 9 5 11 8 37 9 11 5 5 9 10 38 7 11 5 5 9 9 39 7 9 5 4 9 9 40 9 9 8 6 9 9 41 10 9 8 6 9 9 42 10 7 6 6 4 7 43 10 7 8 7 6 6 44 9 8 5 5 9 9 45 10 8 9 5 7 6 46 12 8 9 5 9 6 47 10 10 9 5 9 9 48 10 10 9 5 9 9 49 10 11 9 4 9 9 50 10 9 7 6 9 9 51 10 7 6 6 4 8 52 10 10 5 5 7 9 53 10 6 9 9 9 9 54 10 10 10 7 9 9 55 10 10 9 7 9 11 56 10 12 9 7 9 11 57 10 11 9 4 9 9 58 10 11 9 4 7 9 59 10 9 6 6 8 8 60 10 6 11 9 10 10 61 10 10 10 7 11 9 62 10 10 11 7 11 11 63 10 10 9 7 11 11 64 10 10 9 4 11 9 65 10 11 12 6 12 9 66 10 10 10 7 11 9 67 10 10 10 7 11 9 68 10 10 9 7 9 11 69 10 10 12 7 12 9 70 10 9 12 6 12 10
第三章 模分析方法與分析結果
本章討論經驗正交函數與 Karhunen-Loève 展開法二種模分析法,簡要說明二 種方法所分解出之模組成份、構成以及主要模組之物理意義。另外,亦分析以主 要模組疊加重組之資料與原始資料之誤差及重組還原度。最後討論以模組重組法 填補觀測缺漏值時之實際作法。
3.1 經驗正交函數
經驗正交函數(EOF)是物理海洋學常用的一種統計方法,係將資料分解成空間 和時間上的許多成份,藉以探討影響資料行為的主要因素。這種方法在不同的科 學領域裡有不同的名稱,例如海洋學和氣象學稱為經驗正交函數,統計學稱為主 成份分析(Principal Component Analysis, PCA),心理學或經濟學領域則使用因素分 析(factor analysis)的名稱(Beckers and Rixen, 2003; Kerschen et al., 2005)。Kundu et al.
(1975)提到 EOF 法是個不需要任何動力假設,只需要做資料統計的分析方法。近 幾十年來,氣象學和海洋學方面更是廣泛地運用 EOF 法來分析資料(Beckers and Rixen, 2003)。由於 EOF 法可以從資料裡提取出主要的結構,能有效地將資料用簡 潔的形式,也就是少量的主分量表現出來(Kaihatu et al., 1998),因此成為簡化或過 濾大量資料時常用的有用工具。EOF 法可以用於純量,例如海表溫(Kutzbach, 1967),
或向量資料,如風(Hardy, 1977)或海流(Kundu and Allen, 1976)等。
Kaihatu et al. (1998)指出,對於向量型態的資料,EOF 法有兩種不同的處理方 式,其一如 Kundu and Allen (1976)係使用複數型態的 EOF 法,也稱為 CEOF (Complex Empirical Orthogonal Functions)法,是先將二維的水平流速向量以複數型 態表示,接著再根據EOF 法的步驟求解複數矩陣之特徵向量(各成份可能為複數);
另一種則是稱為實向量EOF (real-vector EOF)法,此方法是使用向量原型直接進行 分析,所求出之特徵向量各成份均為實數。由於 CEOF 法算出的各組特徵向量其 方向並非固定值,但實向量EOF 法卻不會有此種情況,且後者並不侷限於二維,
可以延展至N 維向量(Kaihatu et al., 1998),因此本文採用實向量 EOF 法來進行分 析。
在SCONET 海域內共有 70 個資料點(圖 1-5),每個資料點所測到的都是一組 二維的水平流速向量,在此我們採用的是濾潮後之流速
(
U ,i Vi)
, i=1 ,2, ,70。使用實向量EOF 法計算時,需先將(U ,V )減去各資料點的海流平均值(U ,V ),以二 者 之 差 值(u,v ) 進 行 後 續 分 析 , 為 簡 化 計 下 文 中 u 和 v 將 簡 稱 為 變 動 流 速 (fluctuation)。
( x
it
k) ( U x
it
k) ( ) U x
iu , = , −
( x
it
k) ( V x
it
k) ( ) V x
iv , = , −
(1)在(1)式中, x 表 示 第i i 個 資 料 點 (i=1 ,2, ,70)的水平位置,t 表 示 第k k 小 時 (k =1,2,,K,K =13280)。由(1)式,共70個觀測點而每點又有兩個流速分量,因此 可以組成一個140×140的共變異數矩陣(covariance matrix),以R表示:
[ ]
, , 1,2, , ,1 1
1 1
1 1
1 1
1 1 1 1
1 1 1 1
N q
p
v v u v
v u u u v
v u v
v u u u
v v u v
v u u u v
v u v
v u u u
R R
M M M M
M M M M M
M
M M
M M
M M
pq
=
=
= (2)
式中,
( ) ( )
=
K= kk j k i j
i
u x t u x t
u K u
1
, 1 ,
( ) ( )
==
Kk
k j k i j
i
u x t v x t
v K u
1
, 1 ,
( ) ( )
==
Kk
k j k i j
i
v x t u x t
u K v
1
, 1 ,
. 2 70 ,
140 = =
= N
M N
( x t ) v ( ) x t i j M
K v v
v
K
k
k j k i j
i
1 , , , , 1 , 2 , ,
1
=
=
=
K 為資料筆數。得出R後即可計算出N個特徵值(eigenvalue)λn,n=1,2,,N ,以 及相對應的N組互相正交(orthogonal)的特徵向量(eigenvector)φn,n=1,2,,N;即
N q
N n
R
n nqN
p
p n
pq
, 1 , 2 , , ; 1 , 2 , ,
1
=
=
=
=
φ λ
φ
(3)式中,
φn又可稱為第 n 模之模組(mode)或基底向量(base vector),T表示轉置矩陣。得出 模 組 向 量 後 , 再 計 算 變 動 流 速 對 各 模 組 向 量 之 投 影 , 令
[ ] w
p[ u ( ) ( ) ( ) ( ) x v x u x v x u ( ) ( ) x
Mv x
M]
w = =
1,
1,
2,
2, , ,
(p=1 ,2, ,N),即:( ) ( )
=
= +
=
= + +
=
+
=
+
=
=
N
p
p n k p N
p
p n k p
N
p
p v n k p N
p
p u n k p N
p
p n k p k
n
t x v t
x u
x t
x v x
t x u t w t
E
, 6 , 4 ,
2 2
, 5 , 3 ,
1 2
1
, 6 , 4 ,
2 2 2
, 5 , 3 ,
1 2
1 2
1 1
, ,
, ,
φ φ
φ φ
φ
式中En即為t 時間第k n 模組之振幅。由所有模組振幅之時序資料亦可再行合成變
動流速,即:
( ) ( ) ( ) ( )
[ ]
( ) ( )
( ) ( )
⋅
=
M v n
M u n v n u n
M M
x x x x
x v x u x
v x u
φ φ φ φ
1 1
1
1
, , , ,
[ ] [ ]
( ) ( ) ( ) ( ) ( ) ( )
[
nu nv nu nv nu M nv M]
TN T n n
n T
N p
p n n
x x
x x
x
x φ φ φ φ φ
φ
φ φ
φ φ
φ
, ,
, ,
, ,
, , ,
2 2
1 1
2 1 ,
, 2 , 1
=
=
=
=(4)