應用衛星影像決定波向線及估算水深之初步探討

(1)

國立交通大學

土木工程學系

碩士論文

應用衛星影像決定波向線及估算水深之初

步探討

Preliminary investigation of determining wave rays and

bathymetry using satellite images

指導教授：張憲國博士

研究生：柯紳彥

中華民國一百零一年七月

(2)

(3)

國立交通大學

土木工程學系

碩士論文

應用衛星影像決定波向線及估算水深之初

步探討

Preliminary investigation of determining wave rays and

bathymetry using satellite images

指導教授：張憲國博士

研究生：柯紳彥

中華民國一百零一年七月

(4)

I

應用衛星影像決定波向線及估算水深之初

步探討

研究生：柯紳彥

指導教授：張憲國博士

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

摘要

不同於往昔利用船隻載具單束或多束測深儀來量測海底地形，本文提出以光學衛星影像來決定波向線及估算水深的新技術。本文主要運用直方圖均化及空間濾波提高影像對比度，以影像侵蝕濾除雜訊，再採用Otsu 演算法求得波紋辨識之門檻值，擷取影像圖中的波紋線。垂直於波紋線的定義決定波向線，本文在不考慮因波浪交會能量集中產生波峰變化的情況，以平滑化修正波長，並運用波浪理論中的分散關係式估算海底水深。本文利用屏東海岸實測地形評估本模式決定的波向角度及水深之精度，以本模式推算的波向角度與數模計算所得比較二者的誤差均方根約為 4.5 度，平均絕對值誤差約為 3 度，顯示本方法具有可靠的精度來決定波向角度。至於估算水深方面，本模式估算於-5m 至-7.5m 間的水深較實測水深低，但-7.5m 至-10m 的水深卻略高，而水深-10m 以下的情況，二者估算值接近。本法尚無法完全確實辨認外灘有無沙堆存在。因此，當無外灘沙堆的剖面，估算水深與實測水深的誤差均方根約為0.62m，當有外灘沙堆的剖面，估算水深與實測水深的誤差均方根升高約至0.9m。在明顯繞射區域本文所提出波向需要依繞射理論之波峰線特性修正擷取波紋線之方法。並以龜山島地形遮蔽所產生的波浪繞射現象驗證本文所提決定波浪繞射區的波向線方法是可行的。

(5)

II

Preliminary investigation of determining wave rays

and bathymetry using satellite images

Author：Shen-Yen Ke Advisor：Dr. Hsien-Kuo Chang Institute of Civil Engineering National Chiao Tung University

ABSTRACT

This paper presents a new method of determining wave rays and bathymetry from satellite images instead of the previous bathymetry survey using a vehicle with single/multi-beam echo sounder. Some techniques of image

processing are used to obtain wave rays and water depth. The techniques include: (1) histogram equalization and spatial filtering for increasing image contrast; (2) image erosion for filtering out the noise; (3) Otsu algorithm used for the threshold of extracting crest lines; (4) some algorithms of determining wave rays normal to crest lines; (5) corrected wavelength on one wave ray between two adjacent fronts; (6) neglecting wave acoustics water depth determined by the dispersion relation of the Airy wave theory with known the wave period and wave length. Finally, the whole bathymetry can be obtained at the intersection between all determined crest lines and wave rays.

The proposed method is assessed valid for estimating wave rays and water depth by comparing the results obtained by a numerical model and by bathymetry survey with those obtained by the proposed method at the Pingtung coast. The root mean square error and the mean absolute error between the estimated wave rays and those calculated by the numerical model are about 4.5 degrees and 3 degrees. The estimated water depths by the proposed methods are lower than those by the bathymetry survey for water depth from -5m to -7.5m. However, the estimated water depths by the proposed methods are slightly higher than those by the bathymetry survey for deeper depth from -7.5m to -10m. Both results are approximate when water depth exceeds -10m. The proposed

(6)

III

method cannot exactly and accurately identify the existence of offshore bars. The root mean square error between the estimated water depths and measured ones along a non-dune profile is 0.62m and is 0.9m.along a dune profile.

The proposed method of determining wave rays is corrected to the case in the diffracted zone according to the characteristics of crest lines in the diffracted zone of the theory of water wave diffraction. The corrected method of wave rays is examined good in the diffracted zone of the Guishan Island.

(7)

IV

誌謝

求學期間，承蒙恩師張憲國教授的細心教誨，以邏輯的訓練方式使得學生能夠在學術及各種方面具有獨立思考的能力，亦師亦友的風格讓我獲益良多。感謝老師，在此致上對您最高的敬意。論文口試時，承蒙王順寬教授、莊文傑博士及劉勁成博士提供寶貴意見與指導，經過改正後使本文更趨完善，特此銘謝。感謝研究室的伙伴們及中興大學的各位朋友們，在研究期間給予我莫大的幫助與鼓勵，有你們的陪伴，使得求學的生涯更加多采多姿。感謝我的父母及姐姐，沒有你們就沒有今日的我，謝謝你們的教養之恩，對於我的求學之路總是永遠的支持，使我得以順利完成學業，僅以本文獻給你們。

(8)

V

目錄

摘要 ... I ABSTRACT ... II 誌謝 ... IV 目錄 ... V 圖目錄 ... VI 表目錄 ... VI 第一章緒論 ... 1 1-1 研究動機與目的 ... 1 1-2 文獻回顧 ... 2 1-3 文章架構 ... 4 第二章研究背景及數值模擬簡介 ... 5 2-1 研究範圍介紹 ... 5 2-2 福爾摩沙衛星二號介紹及資料來源 ... 5 2-3 數值模擬簡介 ... 7 第三章波紋辨識及波向推算 ... 10 3-1 影像增強 ... 10 3-2 波紋線的辨識 ... 17 3-3 波向線的決定 ... 26 第四章波長修正及水深推算 ... 33 4-1 波長的修正 ... 33 4-2 水深的推算 ... 35 第五章結論與建議 ... 45 5-1 結論 ... 45 5-2 建議 ... 47 參考文獻 ... 48

(9)

VI

圖目錄

圖 2‐1 (A) 枋寮漁港外側海域之地理位置圖 ... 5 (B) 龜山島右側海域之地理位置圖 ... 5 圖 3‐1 直方圖均化修正圖 ...12 圖 3‐2 (A)範例原圖 (B)平滑濾波後影像圖 (C)增強濾波後影像圖 ...14 圖 3‐3 影像侵蝕範例圖 ...15 圖 3‐4 OTSU 門檻值解說圖 ...16 圖 3‐5 (A)原始 FORMOSAT‐2 衛星影像圖 (B)影像增強圖 ...17 圖 3‐6 (A) 座標轉換示意圖 (B) 影像旋轉前後圖 ...18 圖 3‐7 影像分區示意圖 ...19 圖 3‐8 有效波紋線位置的初步門檻值的示意圖 ...20 圖 3‐9 決定初始波紋線的位置的示意圖 ...21 圖 3‐10 決定同一波紋線的示意圖 ...22 圖 3‐11 合理波紋線的示意圖 ...23 圖 3‐12 初始整體波紋線的示意圖 ...24 圖 3‐13 決定合理整體波紋線的示意圖 ...25 圖 3‐14 波紋線平滑化的修正示意圖 ...26 圖 3‐16 繞射區波向線示意圖 ...27 圖 3‐17 前進區與繞射區分區示意圖 ...28 圖 3‐18 繞射區波向修正範例圖 ...28 圖 3‐19 推算衛星影像的枋寮海域的波向圖 ...30 圖 3‐20 以數模計算枋寮海域的波向圖 ...30 圖 3‐21 推算衛星影像的龜山島上半海域的波向圖 ...31 圖 3‐22 推算衛星影像的龜山島下半海域的波向圖 ...31 圖 3‐23 以數模計算龜山島海域的波向圖 ...32 圖 4‐1 相對水深與相對波長對應圖 ...33 圖 4‐2 波向方向之波長修正範例圖 ...34

(10)

VII 圖 4‐3 原推算波長與修正後等波長之比較 ...35 圖 4‐4 民國 99 年 11 月實測水深分佈圖 ...38 圖 4‐5 衛星影像估算水深分佈圖 ...39 圖 4‐6 選擇 10 個剖面的位置圖 ...39 圖 4‐7 各剖面之水深圖 ...44

(11)

VIII

表目錄

表2-1 商用光學衛星空間解析度 ... 6

表2-2 FORMOSAT-2 衛星資料的基本特性 ... 6

表2-3 選擇影像背景資料 ... 7

表4-1 選擇四個剖面計算出沙丘影響波長的變化 ... 36

表4-2 各剖面的實測水深與推算水深的 RMSE(m) ... 38

(12)

1

第一章

緒論

1-1 研究動機與目的

海岸向來是濱海國家人文發展、經濟活動及生態聚集之地，為一敏感且脆弱的地帶，任何海洋特性的變化皆可能對海岸造成影響，進而破壞人文社會之發展。自古以來，台灣海岸在自然營力（如海浪、海流及風等）及人為力量（如興建水庫、砂石濫採及超抽地下水等）等作用下，各處海岸之侵淤變化多端，有些區域尚屬於穩定情況，但部份海岸已受到嚴重侵蝕，如台南二仁溪北岸灣裡岸一帶，在 1993 年至 2004 年之間，海岸線後退最大可達 100 公尺左右，且侵蝕情況持續惡化中。故海岸侵淤之問題，對於沿岸居民的安全及經濟影響，不容小覷。對於海岸侵蝕問題的處理方式，以人工岬角、突堤、離岸堤或潛堤等海岸保護工法皆能達到減緩海岸侵蝕的作用，但長時間下來，仍會造成沉積物的堆積而掩蓋某些海岸岩台。因此若能正確、有效的設計有利於海岸保護及生態保育之工程規劃，則能減少對海岸造成的衝擊與破壞。工程的事前規劃，需要有海域地形、海象資料、沿岸海流等資訊來設計分析，如波向資料可用來獲得折、繞射現象及海流方向並得知沿岸漂沙之資訊，而水深資料則可用來了解海域地形之變化等。目前國內外對於海洋特性監測資料的獲得主要由海洋研究站、海洋監測浮標、商用光學衛星及海洋調查船等組成，若為實現大範圍海洋監測，則須投入大量海洋監測浮標或增加調查船航行次數，但對於環境惡劣的海域，浮標的放置與回收及船隻量測人員的安全都是非常棘手的問題。以台灣而言，往昔的水深資料多為船隻載具單束或多束測深儀來量測，不僅費時、費力，而且花費金額亦高，亦無法做長期且密集的量測，船隻因海域環境的限制條件，對於許多近岸或礁石過多的海域皆無法取得詳細水深資料，為改善其作業上的缺點，雖在船隻上裝置GPS 定位系統，但對於實際

(13)

2 量測路線規劃上仍然無法精確掌握。近年來高解析度的商用光學衛星紛紛發射，從1999 年 Landsat-7 ETM+ 的15 公尺空間解析度，至 2003 年 SPOT-5 衛星的 10 公尺空間解析度，到目前我國於 2004 年所發射的福爾摩沙衛星二號衛星(FORMOSAT-2)的 2 公尺空間解析度。隨著太空遙測技術的日趨完善，衛星影像具有全面性、機動性及高空間解析度等許多優點，利用衛星影像來辨識波紋線進而推算波浪特性的資料是一種新嘗試的解決方法。以衛星影像圖來取得波向線的方式經濟且省時，經過後續分析處理，能獲得更多的波浪資料來提供海岸施工的初步運用。因此本文希望藉由高空間解析度的FORMOSAT-2 影像，配合影像處理技術來提升影像的對比度，以提高波紋線辨識的精確度，並利用辨識所得的波紋線來推算波向線及計算波長，最後運用波浪理論來求得更多波浪特性資料，以提供港灣建設初步估算使用。本文選擇屏東海岸及宜蘭龜山島海岸二處作為研究基地，並驗證此模式的可行性及準確性，確定模式整體邏輯的可行性後，再對於模式中各小細節做深入的探討與修正，以提高精確度，並運用至台灣其他海域部份。

1-2 文獻回顧

以衛星影像來探討海洋特性的新技術在國內外逐漸被提出，如 Vogelzang 等人(1992) 、Durand 等人(1998)及 Pleskachevsky 等人(2011)，其結果皆顯示高空間解析度的衛星影像在地理資訊學以及河口與海岸調查上有所幫助，並在未來具有很大的發展性。

衛星影像運用在灘線擷取上的探討，Ryu 等人(2002)指出利用多頻譜影像中的熱紅外光(Thermal-Infrared Rays)、近紅外光(Near-Infrared Rays)及短波紅外光(Short-Wave Infrared Rays)的組合，可在退潮時提升衛星影像擷取灘線的成效。吳等人(2003)利用許多不同時間的遙測資料分析台灣西海岸的變遷，顯示港南海岸是台灣西岸侵蝕較嚴重的區域。呂等人(2004)以涵蓋範

(14)

3

圍大但空間解析度差的 SPOT 衛星影像利用區域成長趨近法進行澎湖及東沙島地區的海岸線與面積的遙測分析，在忽略潮位影響後的遙感偵測結果其平均誤差達10％以下。楊(2007)利用 SPOT-5 衛星中之近紅外光波段影像辨識水線。

Chen and Chang(2009)及陳(2010)提出灘線平移修正法(one-line shifting method, OSM)，並透過潮位資料融合法，使 Nao.99b 潮位模式的 RMSE 在 10m 以下，其誤差改善率在西部海岸為 64％，西南海岸為 73％，東部海岸為17％。吳(2010)提出 HIS 影像融合提高影像圖解析度，並以 NOVI-Green 設定灘線門檻值，擷取影像圖之灘線，探討台南的黃金海岸灘線變化，結果顯示黃金海岸為一個夏侵冬淤的海岸型態。張等人(2011)探討光學衛星影像的雲層覆蓋率，並以 GMM 配合群集分析法並同時考慮時間及空間差異對台灣海岸進行分類，提供海岸地區光學衛星影像可用率的初估。衛星影像運用在探討水深的新技術也逐漸被提出，Vogelzang 等人(1992) 以合成孔徑雷達(SAR)影像圖來推算近岸水深可得較高精確度，因 SAR 影像圖具有高空間解析度，而目前 SAR 影像圖已達 1m 以下的空間解析度。 Durand 等人(1998)提出使用衛星遙測數據來推算湖泊或海洋水深的可行性，而Kuo 等人(1999)驗證了從 ERS-1 SAR 的功率譜中找到主波方向和波峰值頻率的可行性。Tripathi 等人(2002)從 IRS-1D LISS-III 獲得 Kakinada Bay 地區的遙測數據，以最小二乘回歸分析海域水深與反射光譜之間的回歸式，再經過濁度校正後，選用Chi-square 來做驗證，其於頻段 0.52μm-0.58μm 的結果最佳。Leu 等人(2005)利用 Wave spectrum bathymetric(WSB)方法從遙測光學影像中獲得淺水域的水深，其檢測水深範圍約略低於12m。

Pleskachevsky and Lehner(2011)提出位於 0m 至 20m 水深可由衛星光學遙測數據獲得，而10m 至 100m 水深處則需利用 SAR 影像來獲得，介於 10m 至20m 的水深則可從以上兩個來源互相應用來取得合理值。SAR 影像與其他衛星光學遙測影像有不同天氣條件的限制，而融合二者的數據在水深

(15)

4 20m 至 60m 之間，其精度可達±15％。由上述研究指出，衛星影像的空間解析度對於推估結果的精確度具有高度的影響性，而不同衛星遙測的優缺點如何互補也是一大課題。若要將衛星影像探討海洋特性資料的技術用於工程上，勢必需要提升更高的空間解析度，以增加其可靠性。

1-3 文章架構

本文以衛星影像圖進行波紋辨識及波向推算後，進而平滑化修正經計算所得的波長，並運用波浪理論中的分散關係式來估算水深。整體架構上依照五個章節來探討研究流程及方法。第一章是前言論述，主要說明本文之研究動機及目的、文獻回顧及文章架構。第二章是簡介資料來源的商用光學衛星，且說明在衛星影像圖選擇上所需要注意的事項，並依序列出本文所選用的研究基地之衛星影像圖時間、衛星種類及衛星的空間解析度，最後再簡介數值模擬MIKE 21 軟體。第三章是介紹影像前置處理的方法，以增強影像圖的對比度，並利用 Otsu 演算法求得波紋線門檻值。接著，本文提出一系列的波紋線辨識流程，再運用辨識所得之波紋線上各座標，依照法線方向求得波向線，並針對繞射區提出波紋線修正的方法。第四章是介紹以波紋線上之波向點計算所得的波長，在本文僅考慮水深為平緩變化下，須做平滑化修正，修正後的波長配合海洋浮標測得的週期，經由波浪理論中的分散關係式來估算水深。第五章是針對本文分析所得的結果，進行探討並給予建議，以利後續之研究探討。

(16)

5

第二章

研究背景及數值模擬簡介

2-1 研究範圍介紹

本文選擇枋寮海岸及龜山島二處作為研究基地。屏東縣枋寮鄉枋寮漁港海域，長度約達2.5 公里，海岸線平直無曲折，水深變化平緩。宜蘭縣外海的龜山島，海岸線多曲折變化，二者的地理位置分別如圖2-1(a)及 2-1(b) 所示。圖 2-1 (a) 枋寮漁港外側海域之地理位置圖 (b) 龜山島右側海域之地理位置圖

2-2 福爾摩沙衛星二號介紹及資料來源

目前商用光學衛星影像逐漸提高空間解析度，如表 2-1 所示。從 1992 年Landsat-7 ETM+的 15 公尺空間解析度、2003 年 SPOT-5 衛星的 2.5 公尺空間解析度，到現在福爾摩沙衛星二號(FORMOSAT-2)以及美國的 IKONOS 與Quickbird 衛星約 1 公尺的空間解析度。然而，每種衛星各有其優缺點，本研究考慮高空間解析度及資料取得方便性，以福爾摩沙衛星二號為主，以下簡述此衛星之特性。

(17)

6 福爾摩沙衛星二號於 2004 年 5 月 21 日成功發射，為我國第一個自主性遙測科學衛星，屬於太陽同步衛星，每日繞地球 14 圈，具左右各 45°傾斜拍攝能力，其掃描像幅寬為24 公里，攝影模式為衛星本體旋轉同步取樣。應用領域包含土地利用與變遷、農林規劃、環境監控、災害評估及科學研究與教育等方面。福爾摩沙衛星二號影像圖的座標系統為TWD67 系統，而實測地形資料為 TWD97 系統，為利後續研究分析，需將影像的 TWD67 座標轉換成 TWD97 系統，以統一所有座標系統。表2-1 商用光學衛星空間解析度空間解析度 15 公尺 2.5 公尺 2 公尺以內衛星名稱 Landsat-7 ETM+ SPOT-5 FORMOSAT-2 IKONOS Quickbird 本文向國立中央大學太空遙測中心申購的福爾摩沙衛星二號影像圖，已經過系統改正及精密幾何改正，其影像空間解析度與感測器光譜模式如表2-2 所示。表2-2 FORMOSAT-2 衛星資料的基本特性感測器軌道光譜模式波段數光譜解析度(μm) 空間解析度 RSI 太陽同步全色態(PAN) 1 0.52~0.82 2m 多頻譜態(XS) 4 0.45~0.52(藍) 0.52~0.60(綠) 0.63~0.69(紅) 0.76~0.90(近紅外) 8m

(18)

7 本文主要辨識衛星影像的波紋線並推算出其波向及波長，為提高辨識的準確度，須採用較高品質的影像圖。若衛星影像有雲霧覆蓋於圖標區上方，波紋線則無法達到有效辨識，而在颱風前後有較大波浪，波紋線較為明顯，因此本文選定無雲、矇氣低及波紋線明顯之影像圖來分析。基於上述條件，影像圖可選取數量少，經過比對選擇後，使用的影像圖資料如表 2-3 所示。表 2-3 選擇影像背景資料研究基地日期時間衛星種類座標系統對應颱風枋寮 2007/08/05 02:08:56 FORMOSAT-2 TWD67 PABUK (08/06~08/08) 龜山島 2011/05/29 02:00:45 FORMOSAT-2 TWD67 SONGDA (05/27~05/28)

2-3 數值模擬簡介

本文採用丹麥水力研究所(DHI：Danish Hydraulic Institute)研發完成之 MIKE 21 軟體進行數值模擬，DHI 係一個獨立國際諮詢及科學研究機構，其所研發之MIKE 21 是模擬水動力、水質、泥沙、波浪之專業工程軟體，主要應用於港口、河流、湖泊、河口海岸及海洋，具有先進前後處理功能及用戶介面。本文主要利用MIKE 21 SW 波浪模式進行推算，以下簡述此模式。 MIKE 21 SW 主要用於模擬離岸及近岸區，風浪與湧浪之成長、衰減及傳播變形。模式以有限體積法(Finite Volume Method)將控制方程式在空間上進行離散，所採用之網格系統為非結構性三角形網格，每一個三角形都代表一個元素，每一時間步均可解出各元素中心之N 值。MIKE 21 SW 包含

(19)

8

1.全譜公式

全譜公式基於 Komen 等人(1994)與 Young (1999)之波動守恆方程式。全譜模式包含下列物理現象：風浪成長、波與波非線性交互作用、白沫消散、底床摩擦消散、水深引起之波浪破碎、水深變化引起之折射與淺水變形、波流交互作用、水深隨時間變化產生之影響、波浪場冰層覆蓋產生之影響。該模式係模擬波浪在各種外力作用下之波能變化，其求解之方程式為波浪作用力守恆方程式(Wave action conservation equation)，此方程式為二維時變域之偏微分方程式，表示如下：

 



S N v t N ___ _   (2-1)

其中N



x,



,



,t



為波力頻譜密度函數(wave action density function)，與波能頻譜密度函數(wave energy density function)之關係為 N=E/σ；x

 

x,y 為卡氏坐標；t 為時間；σ為角頻率；而θ則為波向角。v 



cx,cy,c,c



為一個 4 維空間向量，



cx,cy



dx/dt cg U，代表群波在空間上之傳遞速度，其中c_g為群波速度，而

U



則為海流之流速，模式可考量波浪在海流流動狀況下之傳遞情形。c_ 



/dt，代表角頻率σ之變化速度。c_ 



/dt，代表波向角θ之變化速度。



     



 





  / x, / y, / , / 為一個 4 維空間之向量運算子，S 則為源項(Source term)，代表波浪傳遞過程時可能發生之波浪能量成長、消散及非線性交互作用等物理現象，表示如下： surf bot ds nl in S S S S S S      (2-2) 其中Sin為風浪之成長項，Snl為波與波非線性交互作用項，Sds為白沫消散項， Sbot為底床摩擦項，而Ssurf則為碎波消耗項。

(20)

9

2.方向參數化解耦公式

方向參數化解耦公式基於波動守恆方程之參數化形式，按照 Holthuuijsen (1989)理論，參數化在頻率空間內進行，引入波譜之零階與一階矩作為決定變量，可用於以波浪碎波後而產生之沿岸流之泥沙傳輸計算。參數化之方程式表示如下：

 



 







0 0 0 0 0 C m _T y m C x m C t m gx gy _           



 _(2-3)

 



 







1 1 1 1 1 C m _T y m C x m C t m gx gy _           



 _(2-4)

其中m0(x,y,



)為波譜 N(x,y,ω,



)之零次矩，m1(x,y,



)為波譜 N(x,y,ω,



)之一次

矩，Cgx與 Cgy分別為群波速度 x、y 方向之分量，C為方向之波浪行進速度，



為角頻率，T0與 T1為 Source Terms，包含風、底床摩擦、碎波等因素。 n 次矩 mn(



)定義為：







_





0 , , , , ,y w N x y w dw x m_n



n



(2-5)

(21)

10

第三章

波紋辨識及波向推算

本研究項目主要包含影像增強、波紋線辨識及波向線推算等三部份。各研究項目之目的皆為提升衛星影像圖於波向推算之正確性。本文同時收集枋寮海岸地形資料(民國 99 年度屏東海岸基本資料監測調查計畫報告)，此地形資料可提供數值模擬出波向，並驗證影像圖推算波向及水深的正確性。龜山島雖無實測地形資料以利驗證波向及水深，但因其曲折海岸線會有明顯波浪折、繞射現象，可以說明影像圖推算波向之能力。因衛星影像圖的對比度不強，導致波紋線不易辨識，因此本文採用一系列影像增強方法，提高影像圖對比度，另外，波紋辨識需像素門檻值來區分波峰與波谷，所以本章節提出Otsu 演算法

(Otsu’s method)

選定門檻值。最後本文提出決定波紋及波向線的方法。

3-1 影像增強

為提升波紋辨識的精確度，需先進行影像前置處理。其中包括提高對比度的直方圖均化

(histogram equalization)

及空間濾波，濾除雜訊值的影像形態學之侵蝕法

(erosion)

及判定門檻值的 Otsu 演算法，以下分別簡述各種影像前置處理方法。

1.直方圖均化

直方圖在此表示數位影像圖中，統計灰階度之像素個數的直方分布圖，由於灰階度代表影像的亮度，因此調整直方圖整體分布可達到影像強化之目的。直方圖均化為使原影像圖中灰階度大小從[fmin , fmax]區間，能均勻分佈於[gmin , gmax]區間，以提升影像對比度。其轉換方程式表示如下：

 

f T g  (3-1)

(22)

11 此轉換須滿足下列二個條件： (a)T（f）在 fmin≦f≦fmax 的區間上為單調遞增。 (b)相對於 fmin≦f≦fmax，gmin≦T（f）≦gmax。 其中f 為原始灰階度，g 為輸出灰階度。輸入機率密度函數形式和輸出機率 密度均勻分佈形式，分別如下(3-2)及(3-3)式所示。

 

, f_min f f_max N n f p_f  f   (3-2)

 

_min _max min max

,

1 g

g

p

_g







(3-3) 其中 N 為所有像素個數，nf為灰階度 f 出現的像素個數，pf（f）和 pg（f）

分別為f、g 的機率密度(probability density function)。

根據機率密度理論得知其累積機率密度(cumulative density function)相等，Pg（f）= Pf（f），如式(3-4)。將(3-2)及(3-3)式代入(3-4)式得式(3-5)。圖 3-1 為直方圖均化的範例說明，原圖灰階範圍集中在 70 至 150 之間，經過處理後擴展至0 至 255 之間。

 

_

 



_gg p_g d  _ff p_f d min min     (3-4)



_{ }

      f f i i _g _f _f _f N n g g g min max min min min max , (3-5)

(23)

12 0 1 2 3 4 5 6 7 8x 10 4 gray scale th e nu m be r of pi xe ls 0 50 100 150 200 250 圖 3-1 直方圖均化修正圖

2.空間濾波

空間濾波為遮罩(mask)逐點方式對影像做處理，其遮罩為矩陣形式，陣列大小依照需求而定，常用遮罩矩陣如下(3-6)、(3-7)、(3-8)、(3-9)式所示。空間濾波可模糊影像，降低雜訊值，亦可增強影像，提高對比度。而濾波器類型又分平滑濾波與增強濾波，以下分別做簡述。 (1)平滑濾波低通濾波器為用來模糊影像及降低雜訊，分別有平均濾波與中值濾波。平均濾波是計算遮罩內的像素值平均值後，作為濾波器輸出，中值濾波是將遮罩內的像素值由小至大排列，選取中間值作為濾波器的輸出。平均濾波的效果不只是去除雜訊，還模糊影像的邊緣，而造成影像失真。如果目的只是去除雜訊，中值濾波有較佳效果且不失真。 (2)增強濾波高通濾波器為用來突顯影像的細節或邊緣，達到影像增強目的，其種類有基本高通空間濾波、高頻加強濾波及差分型濾波器三種類型。基本高通空間濾波的中心位置具有較大灰階度的像素，經過處理後，旁邊像素灰 0 1 2 3 4 5 6 7 8x 10 4 gray scale th e numb er of pixe ls 0 50 100 150 200 250

(24)

13 階度之間的差異會被放大。高頻加強濾波是將原始影像乘上大於1 的倍率，再減去影像經低通濾波後的結果，以提高邊緣灰階度的差異性。差分型濾波器是將影像邊緣灰階度放大，達到凸顯物體邊緣輪廓的效果。說明的範例原圖、平滑濾波圖及增強濾波圖分別於圖3-2(a)、3-2 (b)及 3-2 (c)所示。圖(3-2)可看出經平滑濾波後的影像圖，紅色圓圈內的雜訊點相較於原圖同一區域已經降低，而經增強濾波後的影像圖，邊緣區域明顯，黑白對比提高。 3×3 平均濾波器遮罩：           1 1 1 1 1 1 1 1 1 9 1 (3-6) 3×3 高通空間濾波器遮罩：                   1 1 1 1 8 1 1 1 1 9 1 (3-7) 3×3 高通增強濾波器遮罩： , 1 1 1 1 1 1 9 1 1 1 1 9 1 _                      (3-8) 3×3 差分型濾波器行、列遮罩：                           1 0 1 1 0 1 1 0 1 3 1 1 1 1 0 0 0 1 1 1 3 1 (3-9)

(25)

14 圖3-2 (a)範例原圖 (b)平滑濾波後影像圖 (c)增強濾波後影像圖

3.影像形態學-侵蝕

侵蝕為影像形態學中，常見的方法之一，能夠消除影像邊緣的細小突起及區塊間細長的連接，讓影像區域四周往內縮，其運算如(3-10)式所示。

 



w

Y

X



Y

X

_new







_w



(3-10) 其中w 為二值影像點位，X 為處理資料，Y 為結構元素，Θ 為侵蝕運算子， X new為處理資料。說明範例如下圖3-3 所示。 a b c

(26)

15

圖3-3 影像侵蝕範例圖

4.門檻值判定的 Otsu 演算法

Otsu 演算法是由 Otsu 在 1979 年提出的一種自動二值化方法(automatic binary method)，藉由最大化群間變異數(between class variance)或最小化群內變異數(with-in class variance)達到二值化目的。其運算式如下(3-11)式。









2 1 1 2 0 0 2 2 1 1 2 0 0 2

































w b

F

Minimize

F

Maximize

(3-11)

其中_b2為群間變異數(between class variance)，_w2為群內變異數(with-in class variance) ，₀ 為第一群機率值，₀為第一群灰階平均值，₀2為第一群變異數；₁為第二群機率值，₁為第二群灰階平均值，₁2為第二群變異數， 為整體灰階平均值，2為整體變異數。數位影像中，灰階度 1 至 254 依序分別假定為門檻值，依照門檻值將所有像素分為第一群與第二群，再使用(3-11)式計算與門檻值相對應的 F 值。各灰階度皆對應一個 F 值，從所有 F 值中，選擇最大或最小所對應的灰階度作為Otsu 門檻值。所有像素再依照 Otsu 門檻值分成兩群，灰階度大於門檻值的一群給定255(表示白色)，另一群給定 0(表示黑色)，如下圖 3-4 所示。

(27)

16 圖3-4 Otsu 門檻值解說圖影像前置處理為將影像圖陸域部份先行遮蔽，再以直方圖均化，改變灰階度的分布形式，達到影像強化效果。經強化後影像圖，使用平滑濾波柔化影像圖，再利用增強濾波提升對比度。經上述方法後，影像圖雖已提升對比，但相對也產生影響辨識的雜訊值，因此本研究使用影像型態學中的影像侵蝕概念濾除雜訊。為提升辨識精確度，採用門檻值判定的Otsu 演算法，推算Otsu 門檻值，並劃分波紋及非波紋群集，再取波紋值群集做後續分析。本研究以 2007 年 8 月 5 日的枋寮漁港外側海域作為範例，原始 FORMOSAT-2 衛星影像圖與經過處理後影像圖互做比對，二者分別為圖 3-5(a)及 3-5(b)。圖 3-5(a)整體對比度低，波紋線位置模稜兩可，不易判別，相較下圖3-5(b)經影像前置處理後，波紋點位明顯，提高波紋辨識的精確度。因此本研究判定高對比度可提高影像中決定波紋的精確度。

(28)

17 圖3-5 (a)原始 FORMOSAT-2 衛星影像圖 (b)影像增強圖

3-2 波紋線的辨識

經上節前置處理後的影像圖可清楚顯示波紋點的位置，接著需將波紋點位連接成波紋線，再找出曲線之法線方向。處理過程包含影像旋轉、波紋分區辨識、合併波紋、波紋線修正及影像點位復原座標等。各方法及其步驟簡述如下：

1.影像的旋轉

本文需將原影像圖座標（E，N）轉成平行垂直海岸座標（X，Y），以利後續分析。轉換座標示意圖及旋轉影像前後圖分別如圖3-6(a)及圖 3-6(b) 所示，其中X 和 Y 為影像旋轉後座標軸。座標軸轉換方程表示如下：                     N E Y X



cos sin sin cos (3-12) 其中θ 為旋轉角度，逆時鐘為正值，順時鐘為負值。 a b

(29)

18 圖 3-6 (a) 座標轉換示意圖 (b) 影像旋轉前後圖

2. 波紋分區辨識

大範圍的影像圖中因為內含波紋點位數多，可能無法選擇出連續合理的波紋線。因此本節首先將原影像圖分割，進行分區辨識，期望能夠獲得波紋的細節變化，提高辨識的準確度。本文提出三個步驟來處理分區辨識，並簡述如下：

(30)

19 步驟一：分區切割本步驟的目的為適當切割的部份影像來判斷波紋的合理位置。影像分區的示意圖如圖 3-7 所示，以 50×50 像素（pixel，P）大小切割原影像圖，且於同一列各區間有50%重疊，其標註 11 為第 1 列之 1 區，標註 12 為第 1 列之2 區，標註 21 為第 2 列之 1 區，標註 22 為第 2 列之 2 區，重疊 50% 部份劃以斜線表示，X 和 Y 為旋轉影像圖後座標軸，xi和 yi為第 i 分區座標軸。圖3-7 影像分區示意圖步驟二：初始波紋位置的決定擷取第i 區每個波紋位置的 y 座標，由小至大排列後，計算某 y 值上波峰個數（N(y)）及個數平均值（N_y）與母體標準偏差（ ）。本文選擇平_y 均值加n 倍母體標準偏差（N_y  n _y）當為判斷可能波紋線位置的初步門

(31)

20 檻值，其中n 為 1、1.5 及 2 測試後，以 1.5 倍結果佳。若 N(y)大於N_y _y 可視為有效資料。圖3-8 為選擇一個分區範例來說明此步驟的示意結果，圖中圓點表示影像擷取波紋位置，個數平均值N_y=(13+25+…+25+14)/7=18.1，母體標準偏差 ={[(13-18.1)_y 2+…+(25-18.1)2+(14-18.1)2]/7}0.5 =6.4，門檻值 y y N  =24.5。即在圖中，紅色框框為選擇有效資料 y=8、32 及 33 的個數都大於門檻值。圖3-8 有效波紋線位置的初步門檻值的示意圖從上述案例的說明可知，選出三組有效的波紋線的可能位置後，計算各組的 y 值差量（y），即 =32-8=24 及y₁  =33-32=1，並計算差量平均y₂ 值（y），即y=(24+1)/2=12.5。本文再選擇y為判斷初始波紋線平均位置之門檻值。若y大於y，即視有效波紋線的平均位置為可辨識的波紋線；若y小於y則視有效波紋線的平均位置為相近似之同一條波紋線，在此狀況則再取平均值當為此處波紋線的位置。如圖 3-8 中，有效波紋線 N (32) 0 10 20 30 40 50 0 10 20 30 40 x y N(5)=13 N(8)=25 N(10)=12 N(25)=13 N(32)=25 N(33)=25 N(35)=14

(32)

21 及N (33)即可視為相同之波紋線，而 N (8)及 N (32)為二條不同之波紋線。圖3-9 為範例說明，判斷的門檻值y=12.5，決定初始波紋線的y 值為8 及 32.5。圖3-9 決定初始波紋線的位置的示意圖步驟三：合理波紋線位置的決定步驟二選出二組初始波紋線y=8 及 32.5 後，計算各組y 值差量（y），即y=32.5-8=24.5，再計算此差量的平均值（y），即y=24.5。因本範例只找出二組，因此y=y。若可找出多於二組，則y則會不同於y。本文再選擇y的0.15 倍當為搜尋波紋範圍之門檻值，即0.15y=3.68。計算影像圖中所有波紋點的y 座標與各組y值的絕對差值（ y ），本文認定y 小於 y  15 . 0 的波紋點位，視為同一波紋線。圖 3-10 為選取同一條合理波紋線上的位置說明圖。二條實線為初始波紋線，紅色虛線為各組y 值加減0.15y，即32.5±3.68 及 8±3.68，於紅色虛線內圓點皆為同一波紋線。

(33)

22 圖3-10 決定同一波紋線的示意圖如同步驟二，再計算已選擇在同一波紋線可能範圍內的波紋點的個數（N）及其平均值（N_y）。本文再選擇N_y/2為判斷同一條合理波紋線的門檻值。若N 大於N_y /2，即為同一條合理波紋線的位置。再次計算同一波紋線的平均y 值（y），定為同一條合理波紋線的y 座標。如圖 3-10 中，同一波紋線內波紋個數 N(8)=50 及 N (32.5)=64 ，波紋個數平均值 N_y =(50+64)/2=57，門檻值N_y /2=57/2=28.5，而 N(8)> 28.5 且 N (32.5)>28.5，因此兩組為同一條合理波紋線。進而計算兩組的y，即y 1=7.8 及 y2=33.4，此二值則為決定二條合理波紋線的平均y 座標，如圖 3-11 所示。 0 10 20 30 40 50 0 10 20 30 40 x y y=4.32 y=8 y=11.68 y=28.82 y=32.5 y=36.18

(34)

23 圖3-11 合理波紋線的示意圖

3.合併辨識

經上述處理後，本文已取得各分區合理波紋線的 y 座標，再將各分區的波紋線合併成整體波紋線。處理合併的過程有二個步驟，並簡述如下：合併步驟一：初始整體波紋線的決定各分區波紋線的y 座標還原至原 Y 座標，再將分區同一行的 Y 值排列一起，計算行內各組 Y 值間差量（Y ），接著選取每行Y 計算差量平均值（Y ），本文再選擇Y 作為判斷初始整體波紋線平均位置的門檻值。此時將所有分區合理波紋線的Y 值由小至大排列在一起，以最小 Y 值為第一波紋線，若下一個波紋線Y 值與第一波紋線差量大於Y，則視為第二波紋線，若差量小於Y 則忽略此波紋線，而逐次判斷。以此步驟獲得初始整體波紋線的平均Y 座標。如圖 3-12 說明此步驟，其門檻值Y =13，以 Y=10 為第一波紋線，下一波紋線 Y=20 和第一波紋線的差量為 20-10=10，10 小 0 10 20 30 40 50 0 10 20 30 40 x y y=7.8 y=33.4

(35)

24 於Y ，因此繼續往下一波紋線 Y=25，計算其差量為 25-10=15，15 大於Y，則Y=25 視為第二點，以此類推，紅色框框即為選出的初始整體波紋線平均 Y 座標位置，其中 Y=10 及 Y=20 二個波紋線間距太短，不可能代表是二個有波長距離的波紋線。圖3-12 初始整體波紋線的示意圖選定影像圖最接近初始整體波紋線的平均 Y 值的座標為波紋線起始座標（X，Y）。各波紋線從起始座標開始平行 X 軸方向搜尋，限制條件為初始整體波紋線各組起始 Y 座標的差量 0.5 倍（

0 .

5 

Y

）及相同 X 座標下的 Y 值與各起始 Y 座標最小的差量，則可取得各波紋線內所有點位。合併步驟二：整體波紋線的決定如同前述分區波紋線位置的決定步驟，計算整體波紋線。同一波紋線的Y 座標平均值（Y ）及Y 的差量（



Y

）與其平均值（Y ），本文選定Y 0 10 20 30 40 50 0 20 40 60 80 100 X Y Y=10 Y=20 Y=25 Y=43 Y=62 Y=78 Y=90

(36)

25 為判斷合理整體波紋線的平均位置為門檻值。若



Y

大於Y 則為合理整體波紋線的平均位置，若



Y

小於Y 則合併兩初始波紋線的平均位置為合理整體波紋線的平均位置。圖3-13 為一部份的決定合理波紋線範例，其標號 Yi代表各波紋線段座標平均值，Y2和Y3因差量小於門檻值，故合併成虛線部份。圖3-13 決定合理整體波紋線的示意圖

4.波紋線凹凸段修正

本文取得合理整體波紋線後，波紋線部份線段可能會上下變化太大，波紋線呈現不平滑現象。本文不針對波浪交會而產生局部能量集中，波峰會變形的特殊情況，僅處理平滑波紋線的情況。因此當擷取波紋線有凹凸變化太大時則進行平滑處理。在同一波紋線上的點位斜率變化量應為平緩，如圖3-14 為有凹凸情況的波紋線示意圖，實線代表原波紋線，而在a、b 點之間原波紋線的斜率從正值轉為負值，再轉為正值。因此當斜率變化太大，則代表波紋線有凹凸現象，本文選擇斜率變化量為0.2 內作為門檻值。首先選擇波紋線上最接近平均 Y 座標的點位為起點，X 軸自左至右計算波紋線點位的斜率，不符合門檻值者，進行三次多項（cubic spline）內插平滑波紋線。此內插方法精確度高，但計算時間較線性內插多。

(37)

26 圖3-14 波紋線平滑化的修正示意圖

5.影像座標復原

經上述四個步驟處理後，已獲得最後整體波紋線的座標。再利用式(3-12) 將X、Y 座標轉換回 E、N 座標，即完成波紋線的辨識。

3-3 波向線的決定

1.法線方向的決定

經上節獲得的同一波紋線座標後，從最外海的波紋線為決定波向的起始波紋線。對波紋線上每一座標（Xi，Yi），先進行計算鄰近二點之斜率（tan  ）。如式 3-13 所示，在二端點位使用鄰近二點的割線（secand）當為斜率，而在中點點位則用旁邊二點的割線當為斜率。法線向量(normal vector)則為tan



i  j。此結果可利用matlab 中，梯度運算（gradient 指令）獲得法線向量，而法線方程式(normal equation)為Y Y_i [(X  X_i)/tan]。

(38)

27 ) /( ) ( tan 1 ,..., 2 ) /( ) ( tan 1 ) /( ) ( tan 1 1 1 1 1 1 1 1 1                        i i i i i i i i i i i i X X Y Y n i if n i X X Y Y n i if X X Y Y i if



(3-13) 通過（Xi，Yi）沿著法線方向之直線與下一條波紋線相交的X 座標和 Y 座標，即為波向前進至新的波紋線上的點位。若依上述方法再計算此處法線向量，此值與原來（Xi，Yi）平均後的法線方向，再修正下一個波紋線之交點，如此重覆步驟後，則可求得所有波向線與波紋線的交叉點位。但為平順波向線，本文最後用三次多項式內插法來平滑波向線。

2.繞射區（diffracted zone）波向修正

影像圖中，繞射明顯的區域波向的轉彎比一般折射明顯，如圖3-16 所示。因此處理方法不同於上述波向前進的步驟。本研究提出另外一種決定波向線的方法。首先決定一繞射點，並利用決定的波向線計算平均波長的2 倍（2L）為基準，將繞射點向外移動，往內、外分為繞射區及前進區，如圖 3-17 為前進區與繞射區的示意圖，紅點為繞射點，綠點為移動2L後位置。圖 3-16繞射區波向線示意圖

(39)

28 圖3-17 前進區與繞射區分區示意圖繞射區內，波紋線為弧形，因此本文選用各段平均波長（L_i）為半徑，繞射點移動2L後的位置為圓心，分別畫弧形線，如圖3-18 所示。L₁、L₂ 及 3 L 分別代表各段平均波長，藍線為修正後波紋線。最後將繞射區內修正的波紋線取代原先的波紋線，再重新推算波向線。圖3-18 繞射區波向修正範例圖

(40)

29

為評估推算波向的精確度，本文選用數模計算的波向角度與推估波向角度的均分根誤差(root mean square error，RMSE)。RMSE的定義如式(3-14)

所示。





N d d RMSE N i



 m c   1 2 (3-14) 其中d 為推算波向角度，_c d 為數模計算波向角度，_m N為資料數。此值越大表示二者偏差越多，即推估誤差越大，越不吻合本文推算96年8月5日屏東縣枋寮漁港外側海域的波向圖示如圖3-19，而數模以99 年7月地形計算所得的波向圖示如圖 3-20。此地區地形平緩，灘線筆直，故波向方向變化不大。外海波浪的波向角度約為 226 度，本文所提供的方法推算出的波向角度介於 211 度至 235 度之間，而數模計算所得的波向角度介於 214 度至 232 度之間，二者的平均絕對值誤差量為 3.14 度，RMSE為4.13度。本文劃分龜山島的研究範例圖為上半部與下半部分別探討。推算龜山島上下半部的波向圖分別如圖3-21與圖3-22所示，而數模計算所得的波向圖示如圖 3-23。外海波浪的波向角度為 58 度，圖 3-21 中龜山島上半部因海岸形狀遮蔽，波浪能量側向傳播而產生波浪繞射，本文所提出的方法可有效推算出此區的繞射現象，其推算的波向角度介於 25 度至 60 度之間，而數模計算所得的波向角度介於 32 度至 58 度之間，二者的平均絕對值誤差量為2.32度，RMSE為3.75度。圖3-22中龜山島下半部海岸也因海岸形狀的遮蔽關係，波浪也會有繞射現象，本文亦推得此區的繞射現象，其推算的波向角度介於53 度至 86 度之間，而數模計算所得的波向角度介於 55 度至85度之間，二者的平均絕對值誤差量為3.11 度，RMSE為4.64度。

(41)

30

圖 3-19 推算衛星影像的枋寮海域的波向圖

(42)

31

圖3-21 推算衛星影像的龜山島上半海域的波向圖

(43)

32

(44)

33

第四章

波長修正及水深推算

4-1 波長的修正

經過波紋辨識與波向推算後，獲得波向線點位，再由各波紋線上的波向點位可計算出各段波長。本文僅考慮水深變化為平緩情況，波長會隨水深變淺而逐漸遞減，示如圖4-1。因此假設水深為平滑變化，計算的波長需加以修正，其中包含波向方向修正及波紋方向修正二部份，並簡述如下：圖4-1 相對水深與相對波長對應圖

1.波向方向修正

同一波向線的波長變化量應為平緩，因此對波長進行平滑修正。首先計算同一波向線的波長平均值（L）及母體標準偏差（_L），再從波向線的波長中選擇與L差距最小的波長當為起始點，往上或下的波長總變化量為 1.5倍母體標準偏差作為判斷合理波長的門檻值（1.5_L /n_up、1.5_L /n_down）， 0 0.1 0.2 0.3 0.4 0.5 0.6 0 0.5 1 1.5 2 2.5 h/L0 L/ L0

(45)

34 其中n_up及n_down 分別為起始點往上下的波長個數。若波長變化小於門檻值則視為合理波長，反之波長變化大於門檻值則修正波長。圖4-2為選取某條波向線來當修正波長範例，右邊L值為原波長，位置以紅色圈表示，左邊L 值為修正後波長，位置以綠色方形表示。波長 L=13、20、30 及 35，而平均波長L =24.5，母體標準偏差1.5_L=14.8，L=20為最接近L的波長，故視為起始點，其上下門檻值1.5_L /n_up=14.8/1=14.8及1.5_L /n_down=14.8/2=7.4。起始點向上判斷，L=13 時，20-13=7 小於 14.8，視其為合理波長。起始點向下判斷，L=30 時，30-20=10 大於 7.4，則利用起始點與另一波長作為修正，L=（20+35）/2=27.5，再次判斷此值是否合理，27.5-20=7.5 大於 7.4，則利用起始點加上門檻值再修正，L=20+7.4=27.4，並選定L=27.4 為新起始點，L=35時，35-27.4=7.6大於7.4，但因L=35為最後一波長值，因此修正方法為起始點加上門檻值，L=27.4+7.4=34.8，完成波長修正。圖 4-2 波向方向之波長修正範例圖 0 2 4 6 8 10 10 20 30 40 profile w ave le ng th (m ) L=13 L=20 L=30 L=35 L=13 L=20 L=27.4 L=34.8

the original wave length the new wave length

(46)

35

2.波紋方向修正

同一波紋線的波長值應相近，因此對波長進行平滑修正。首先計算同一波紋線的波長平均值（L）及母體標準偏差（_L），若波紋線上L與L差值小於_L，視為合理波長，若波紋線上 L與L差值大於_L，則 L修正為L，重覆修正直至波長的母體標準偏差變化量小於0.2，完成波長修正。經上述二種波長修正處理前後，如圖4-3 所示。從圖4-3(a)中發現，未經修正前其波長於波向方向及波紋方向的變化量大，而圖4-3(b)修正後的等波長圖顯示波長變化量小。 (a) 原推算波長圖 (b)修正後等波長圖圖4-3 原推算波長與修正後等波長之比較

4-2 水深的推算

經波長修正後，本文再引用波浪理論中的分散關係式，配合中央氣象局的海洋浮標實測週期資料，估算水深值，如(4-1)式。本文僅做水深初步探討，故以單一週期來做計算。 kh k g c2  tanh (4-1)

(47)

36

其中c=L/T為波速(wave celerity)，g為重力加速度(gravitational acceleration)，

k=2

π

/L為波數(wave number)，h為水深(water depth)，L為波長，T為週期。

本文獲得屏東縣枋寮漁港外側海域的實測水深和推算水深示於圖 4-4 及圖 4-5。圖 4-4 中南側近岸標記橢圓處為有沙丘地形區域，但在圖 4-5 相對區域並未判斷出此處有沙丘地形。為探討原因，本文選取四個經過沙丘的剖面出來探討，示如圖4-4，而以微小振幅波理論計算沙丘影響波長變化如表4-1。表4-1之第二及三行分別為沙丘平均水深的波長與沙丘底部水深的波長，第四行為二者之差值（L），分別為 2.1m、4.6m、9.5m及7.3m，而此區的平均波長（L）約為 100m，沙丘寬度（

W

）分別為 50m、70m、 90m及70m，由表 4-1可知波長在沙丘地形影響下波速變慢。考慮沙丘寬度影響波長的比例，依LW /L 計算波長因為沙丘之改變量估算值為1.05m、 3.22m、8.55m及 5.11m。本文選用的衛星影像圖之解析度為2m，對於1.05m 的波長改變量辨識上不易，而3.22m、8.55m及5.11m 則可辨識出來，但因平滑波長處理的過程下，原本辨識所得的沙丘現象卻被修正去除。表4-1 選擇四個剖面計算出沙丘影響波長的變化沙丘平均水深的波長沙丘底部水深的波長  (m)L

W

(m) L(m) L W L /  (m) 剖面1 78.3 80.4 2.1 50 100 1.05 剖面2 81.8 86.4 4.6 70 100 3.22 剖面3 85.9 95.4 9.5 90 100 8.55 剖面4 83.6 90.9 7.3 70 100 5.11 本文再沿海岸切10個水深剖面，探討各剖面的 RMSE來說明本文所提方法推估水深的精確性。選擇剖面的位置如圖 4-6 所示。各水深剖面圖及 RMSE值表分別如圖4-7及表4-2所示。因本文方法無法推估至岸邊的波長，而無-5m以上水深的數據，若使用外插法去推估-5m以上的水深，其結果不佳，因此本文僅列出-5m以下水深的部份。從10個水深剖面圖中一致發現，估算水深於-5m至-7.5m之間較實測水深低，-7.5m至-10m 之間略高於實測水深，最後-10m以下的水深值二者較吻合。

(48)

37 兩次實測地形的時間分別為99年7月及99 年11月，時間在間隔三個月內，地形高程的RMSE在 0.25m至0.66m間。以下依特性分段說明。 1. 剖面1至剖面3為無沙丘情況，波長未受沙丘影響，平均RMSE 約為 0.27m。 2. 剖面4至剖面6為有沙丘情況，而沙丘的位置變化平均為50m，高度變化平均為0.5m，平均RMSE約為 0.35m。 3. 剖面 7 及剖面 8 為有沙丘情況，而沙丘的位置變化平均為 100m，高度變化平均為2m，平均RMSE約為 0.55m。 4. 剖面9為有沙丘情況，沙丘的位置無變化，其高度變化為2.5m，RMSE 約為0.5m。 5. 剖面10為有沙丘情況，沙丘的位置變化為100m，其高度變化為 1m， RMSE約為 0.66m。本文選用的衛星影像圖時間為96年8月，與實測地形資料的時間相隔三年，地形高程的RMSE 可達 0.4m到 1.1m之間的變化。以下依特性分段說明： 1. 剖面 1 至剖面 3，實測水深與估算水深都無沙丘影響，二者的平均 RMSE約為 0.62m。 2. 剖面4、剖面 5及剖面10，實測水深有沙丘現象，其沙丘平均高度約為 2m，而於估算水深的同一區域並無沙丘現象，二者的平均 RMSE 約為0.75m。 3. 剖面 6至剖面9，實測水深有沙丘現象，其沙丘平均高度約為 3.5m，而於估算水深的同一區域並無沙丘現象，二者的平均 RMSE 約為 0.9m。

(49)

38 表4-2 各剖面的實測水深與推算水深的RMSE(m) 與第一次實測比較與第二次實測比較二次實測比較沙丘有無剖面1 0.643 0.584 0.283 No 剖面2 0.609 0.486 0.307 No 剖面3 0.558 0.405 0.254 No 剖面4 0.785 0.653 0.368 Yes 剖面5 0.745 0.733 0.336 Yes 剖面6 0.879 0.914 0.352 Yes 剖面7 1.052 1.023 0.579 Yes 剖面8 0.929 0.973 0.444 Yes 剖面9 1.182 0.772 0.506 Yes 剖面 10 0.753 0.820 0.663 Yes 圖4-4 民國 99年11月實測水深分佈圖 P1 P2 P3 P4

(50)

39 圖4-5 衛星影像估算水深分佈圖圖4-6 選擇10 個剖面的位置圖 X Y water depth 200 400 600 800 1000 1200 1400 0 200 400 600 800 1000 1200 -20 -18 -16 -14 -12 -10 -8 -6 -4 -2 0 p1 p2 p3 p4 p5 p6 p7 p8 p9 p10

(51)

40 200 400 600 800 1000 1200 -20 -15 -10 -5 0 distance(m) w a te r de pt h( m ) PROFILE1 estimation (96/08) measurement (99/07) measurement (99/11) 200 400 600 800 1000 1200 -20 -15 -10 -5 0 distance(m) w a te r de pt h( m ) PROFILE2 estimation (96/08) measurement (99/07) measurement (99/11)

(52)

(53)

(54)

(55)

44 圖4-7 各剖面之水深圖 200 400 600 800 1000 1200 -20 -15 -10 -5 0 distance(m) w a te r de pt h( m ) PROFILE9 estimation (96/08) measurement (99/07) measurement (99/11) 200 400 600 800 1000 1200 -20 -15 -10 -5 0 distance(m) w a te r de pt h( m ) PROFILE10 estimation (96/08) measurement (99/07) measurement (99/11)

(56)

45

第五章

結論與建議

5-1 結論

本文提出以光學衛星影像來獲得波浪資料及海底地形資料的新技術。因衛星影像具有對比低的缺點，本文運用直方圖均化及空間濾波來提高影像對比度，以影像形態學中的影像侵蝕濾除雜訊，再採用Otsu演算法來求得波紋辨識的門檻值，以擷取合理波峰點。本文將原影像座標（E，N）轉換成垂直平行海岸座標（X，Y），為了防止大範圍辨識下，小細節變化易忽略的缺點，以 50×50 像素的矩陣大小分區辨識。分區辨識上，本文以不同資料特性所定義之門檻值判斷方式來分析，獲得波紋線合理的平均 y 值後，再將波紋線的平均 y 值還原至原 Y 值。接著，本文計算各Y值間差量（



Y

）及差量平均值（Y ），以Y 來判斷合理整體波紋線，並以垂直於波紋線的定義決定波向線。本文在不考慮波浪交會產生能量集中所造成的波峰變化之情況，以平滑化修正波長，最後將影像座標（X，Y）轉回原圖座標（E，N）。本文以求得的波向線點位計算各段波長，假設水深為平緩變化，根據波浪理論中相對水深與相對波長對應關係，波長會隨水深改變而平緩的變化，分別以同一波向線的波長平緩變化及同一波紋線的波長之變異量小二種方式修正波長，二者皆是利用波長的平均值（L）及母體標準偏差（L）來做修正。本文再以修正後的波長，根據波浪理論中的分散關係式及海洋浮標獲得的週期來估算，即可獲得水深資料。本文利用屏東海岸的實測地形評估本模式決定的波向角度及水深之精度，以本模式推算的波向角度與數模計算所得比較二者的誤差均方根約為 4.5度，平均絕對值誤差約為3度。顯示本方法具有可靠的精度來決定波向角度。至於估算水深方面，從實測地形得知枋寮海域南側外灘有沙堆存在。本文針對沙堆寬度對波長變化可能造成的影響性，以微小振幅波理論進行

(57)

46 評估，結果顯示沙堆造成的波長變化量約大於2m，相對於衛星影像圖的空間解析度而言，是可辨識出沙堆現象。但本法尚無法完全確認未修正的波長於某區域所產生的變化量是否合理存在，故假設地形為平緩變化下，針對波長進行平滑修正，導致估算的海底地形於同區域並無沙堆現象存在。本文選取10個水深剖面進行精確度評估，因本文無法推估至岸邊的波長，故無-5m以上水深的數據，在此僅列出-5m以下水深部份來探討。本文估算於-5m至-7.5m間的水深較實測水深低，但-7.5m至-10m的水深卻略高，而水深-10m以下的情況，二者估算值接近。二次實測地形的時間間隔三個月，當無外灘沙堆現象的剖面，二次實測的誤差均方根約為0.27m，當有外灘沙堆現象的剖面，隨著沙堆的高度與位置之變化不同，二次實測的誤差均方根約為 0.35m 至 0.66m 之間。而本文選取的影像時間與實測資料時間間隔約三年，當無外灘沙堆現象的剖面，估算水深與實測水深的誤差均方根約為0.62m，當有外灘沙堆現象的剖面，隨著沙堆的高度變化，估算水深與實測水深的誤差均方根約提高為 0.75m 至1.0m之間。由上述得知，本文估算的水深與實測水深資料，整體誤差均方根約為1m，顯示本方法估算水深精確度的可靠性佳。在明顯繞射區域本文所提出波向需要依繞射理論之波峰線特性修正擷取波紋線的方法。繞射區域的波向轉彎比一般折射區域明顯，本文以各段平均波長為半徑，繞射點外移2L為圓心畫弧，重新擷取波紋線，並以龜山島地形遮蔽所產生的波浪繞射現象驗證本文所提決定波浪繞射區的波向線方法是可行的。本文將龜山島劃分上半部與下半部分別探討，上半部以本模式推算的波向角度與數模計算所得比較二者的誤差均方根約為 3.75 度，平均絕對值誤差約為 2.32 度，下半部以本模式推算的波向角度與數模計算所得比較二者的誤差均方根約為4.64度，平均絕對值誤差約為3.11 度，顯示本方法具有可靠的精度來決定繞射區的波向角度。

(58)

47

5-2 建議

本文透過遙測技術及影像處理的整合，進而考慮波浪特性來推算波向線並獲得水深，以期望能在海岸工程規劃上提供便利的監測技術。由於涉及的範圍廣泛，整體研究流程中可能會發生的誤差難以估算。從目前研究成果來看，本文針對可改進及修正的部份，提出探討與建議，以利後續研究方向，簡述如下。 1. 影像處理的方法，需能更準確的擷取出波紋點，降低不合理波紋點之數量，以提升波紋線辨識的精確度。 2. 波向線的推算，本文採用各點斜率逐步推算，隨著推算的長度越長，誤差累積越大，若能採用其他方式來推估，則可提高波向線的精確度。 3. 繞射點的給定方式，應採用自動判斷的方式找出合理位置，以增加模式的客觀性。 4. 波長修正部份，可採用近岸區和海域區來做不同修正方式，對於分界的離岸長度，建議可用沙丘存在的範圍來做限制。 5. 計算水深僅用單一週期配合分散關係式做初估，可深入探討計算水深之不同公式的準確性。 6. 同一區域選取不同時間的多張衛星影像圖去估算水深，並整合及決定出此區域合理之水深。

(59)

48

參考文獻

1. 郭一羽(2001)，「海岸工程學」，文山書局。 2. 吳哲榮、吳啟南(2003)，「遙測技術應用於臺灣西海岸五十年來變遷分析」，航測及遙測學刊，第八卷，第三期，第95-110頁。 3. 呂黎光(2004)，「海岸線與面積之遙感測繪應用研究」，海洋工程學刊，第四卷，第一期，第71-88 頁。 4. 朱子豪(2004)，「數值遙測判釋」，國立台灣大學地理環境資源學系課程講義。 5. 楊勤儀(2007)，「利用衛星影像萃取近岸地形-以台灣北部為例」，國立中央大學地球物理研究所，碩士論文。 6. 陳蔚瑋(2009)，「衛星影像的灘線辨識及其應用至灘線變遷之研究」，國立交通大學土木工程研究所，博士論文。 7. 吳政杰(2010)，「衛星影像灘線辨識之精確度評估研究」，國立交通大學土木工程研究所，碩士論文。 8. 張憲國、陳蔚瑋、劉勁成(2011)，「台灣海岸地區雲量分布與衛星影像可用率分析」，第三十三屆海洋工程研討會論文集，高雄，第 699-704 頁。

9. Alpers, W. and Hennings, I. (1984) A theory of the imaging mechanism of underwater bottom topography by real and synthetic aperture radar, Journal

of Geophysical Research, 89 (C6), 10529–10546.

10. Chen, W. W. and Chang, H. K. (2009) Estimation of shoreline position and change from satellite images considering tidal variation, Estuarine, Coastal

and Shelf Science, 84(20), 54-60.

11. Durand, D. and Cauneau, F. (1998) Towards a new method for

shallow-water monitoring using remote sensing, Future Trends in Remote

(60)

49

12. Hongxing, L. (2003) Derivation of surface topography and terrain

parameters from single satellite image using shape-from-shading technique,

Computers & Geosciences, 29, 1229-1239.

13. Ingo, H. and Margitta, M. (2002) The influence of quasi resonant internal waves on the radar imaging mechanism of shallow sea bottom topography,

Oceanologica Acta, 25, 87-99.

14. Kuo, Y. Y. and Leu, L. G. and Kao, Y. L. (1999) Directional spectrum analysis and statistics obtained from ERS-1 SAR wave images, Ocean

Engineering, 26, 1125-1144.

15. Leu, L. G. and Chang, H. W. (2005) Remotely sensing in detecting the water depths and bed load of shallow waters and their changes, Ocean

Engineering, 32, 1174-1198.

16. Otsu, N. (1987) A threshold selection method from gray-scale histogram,

IEEE Trans. Systems, Man, and Cybernetics, 8, 62-66.

17. Pleskachevsky, A. and Lehner, S. (2011) Estimation of underwater

topography using satellite high resolution synthetic aperture radar data,

Proceedings of 4th TerraSAR-X Meeting, Oberpfaffenhofen, Germany.

18. Ryu, J. H., and Won, J. S. and Min, K. D. (2002) Waterline extraction from Landsat TM data in a tidal flat: A case study in Gomso Bay, Korea, Remote

Sensing of Environment, 83, 442-456.

19. Tripathi, N. K. and Rao, A. M. (2002) Bathymetric mapping in Kakinada bay, India, using IRS-ID LISS-III data, International Journal of Remote

Sensing, 23 (6), 1013–1025.

20. Tay, L. T. and Daya Sagar, B. S. and Chuah, H. T. (2005) Derivation of terrain roughness indicators via granulometries, International Journal of

(61)

50

21. Tuncay, K. and Abdulaziz, G. and Fevzi, K. and Mustafa, D. (2011)

Automatic detection of shoreline change on coastal Ramsar wetlands of Turkey, Ocean Engineering, 38, 1141-1149.

22. Vogelzang, J. (1992) Sea bottom topography with X-band SALR: the

relation between radar imagery and bathymetry, International Journal of

Remote Sensing, 13 (10), 1942-1958.

23. Volker, L. and Christian, H. and Randolph, K. (2006) Derivation of

planetary topography using multi-image shape-from-shading, Planetary

and Space Science, 54, 661-674.

24. Walter, H. F. Smith and David, T. S. (1994) Bathymetric prediction from dense satellite altimetry and sparse shipboard bathymetry, Journal of

Geophysical Research, 99 (B11), 803–824.

25. Yijun, H. and Hui, S. and William, P. (2006) Remote sensing of ocean waves by polarimetric SAR, Journal of Atmospheric and Oceanic

應用衛星影像決定波向線及估算水深之初步探討

國立交通大學

土木工程學系

碩士論文