國
立
交
通
大
學
土木工程學系
碩士論文
加權基本不震盪法結合類神經網路與遺傳
演算法應用於一維淺水波方程式之求解
Applying Weighted Essentially Non-oscillatory
Schemes combined with Artificial Neural
Network and Genetic Algorithm to Solving 1-D
Shallow Water Equations
學生：李勁頤
指導教授：葉克家 博士
加權基本不震盪法結合類神經網路與遺傳
演算法應用於一維淺水波方程式之求解
學生：李勁頤 指導教授：葉克家
國立交通大學土木工程學系
摘要
本研究以原始高階 WENO5 算則(Jiang and Shu 1996)為基礎所提 出在強震波區域的改進方法：修正平滑指示器 WENO5 算則(Zhang and Shu 2007)，將之應用於求解一維淺水波方程式，發現在某些情況 下，會產生較不好的模擬結果。為改善此模擬結果，本研究以修正平 滑指示器 WENO5 算則的收斂理念為基礎，結合類神經網路(ANN) 與遺傳演算法(GA)的原理，透過學習嘗詴發展出一套以類神經網路來 判斷收斂程度，間接決定 WENO 算則權重的新算則。最後以一簡單 之平滑曲線的微分案例與各種一維明渠流案例進行模擬，比較新算則 與各式 WENO5 算則所得之結果，以評估新算則模擬之收斂性與精確 性。 關鍵詞： 加權基本不震盪法、類神經網路、遺傳演算法、一維淺水波方程式
Applying Weighted Essentially Non-oscillatory Schemes
combined with Artificial Neural Network and Genetic
Algorithm to Solving 1-D Shallow Water Equations
Student: Chin-I Li Advisor: Keh-Chia Yeh Institute of Civil Engineering
National Chiao Tung University
Abstract
The WENO5 scheme with the modified smoothness indicators (Zhang and Shu 2007) derived from the original WENO5 scheme (Jiang and Shu 1996), is a correction method used in the region where shock waves happen. By applying this scheme to solving 1-D shallow water equations, we can find that the simulations are not precise enough under some situations. For the more precise result, this study uses the convergent theory based on the WENO5 scheme with the modified smoothness indicators, and combines the principles of artificial neural network (ANN) and genetic algorithm (GA). Through the learning process, we try to develop a new scheme using ANN for determining convergent degrees, and deciding weight distributions of the WENO5 scheme at second hand. Finally, by testing a simple differential case of a smooth curve and some cases of 1-D open-channel flow, we can compare the simulation results of using the new scheme with those of using various WENO5 schemes, as a way to evaluate the convergence and accuracy of the new scheme applied to the simulations.
謝誌
首先要感謝我的父母能養育我至今，並支持我走完這段求學生涯； 而在研究所的這段日子裡，無論是教授、學長、同學、學弟妹與朋友 們，我非常幸運地能遇到一群有趣又心地善良的人們，在此就祝福各 位好人都能有個幸福快樂的未來囉。最後就容小弟我說一聲：謝謝啦！
目錄
摘要 ... I Abstract ... II 謝誌 ... III 目錄 ... IV 表目錄 ... VI 圖目錄 ... VII 符號說明 ... IX 第一章 緒論 ... 1 1-1 研究動機與目的 ... 1 1-2 文獻回顧 ... 3 1-2-1 類神經網路與遺傳演算法 ... 3 1-2-2 高解析算則的發展 ... 4 1-2-3 守恆定律的探討 ... 8 1-3 研究方法與流程 ... 9 第二章 理論基礎 ... 12 2-1 ENO 算則 ... 12 2-1-1 一維數值通量的建構與近似模擬 ... 12 2-1-2 各計算元的數值通量 ... 16 2-1-3 推展與通量分離 ... 18 2-2 WENO 算則 ... 20 2-2-1 上風中央算則與理想權重 ... 20 2-2-2 平滑指示器 ... 22 2-2-3 WENO5 算則 ... 23 2-3 差分式 WENO 算則 ... 25 2-4 Mapped WENO 算則 ... 25 2-5 修正平滑指示器 WENO 算則 ... 27 第三章 WENO 運用於一維淺水波方程式 ... 31 3-1 統御方程式 ... 31 3-2 數值方法 ... 32 3-2-1 時間離散 ... 32 3-2-2 空間離散 ... 33 3-2-3 守恆形式差分式 WENO 算則 ... 34 3-3 乾溼交界處理 ... 40 3-4 稀疏波 ... 41 第四章 模式建置 ... 444-2 模式架構 ... 52 4-2-1 水理模式 ... 52 4-2-2 類神經網路模式 ... 52 4-3 學習案例 ... 55 4-3-1 案例選擇 ... 55 4-3-2 參數設定 ... 57 第五章 模式驗證與結果分析 ... 61 5-1 學習結果 ... 61 5-2 無劇烈通量變化之詴驗案例 ... 67 5-2-1 簡單平滑曲線微分案例 ... 67 5-2-2 穩態流流經障礙物案例 ... 69 5-2-3 潮汐波案例 ... 72 5-3 水面分離案例 ... 74 5-4 潰壩案例 ... 80 第六章 結論與建議 ... 84 6-1 結論 ... 84 6-2 建議 ... 86 參考文獻 ... 88 附錄 A. 類神經網路 ... 97 附錄 B. 遺傳演算法 ... 103
表目錄
表 2- 1 重建常數表(k=1-4) ... 18 表 4- 1 JSWENO 與 ZSWENO 算則結果誤差表 ... 49 表 4- 2 方案 A ~ C 選用常數與代表誤差表 ... 58 表 5- 1 類神經網路學習結果誤差表 ... 61 表 5- 2 100 格網數誤差表 ... 67 表 5- 3 500 格網數誤差表 ... 67 表 5- 4 1000 格網數誤差表 ... 68 表 5- 5 2000 格網數誤差表 ... 68 表 5- 6 案例 A 與案例 B 之誤差表 ... 70 表 5- 7 潮汐波案例之誤差表 ... 73圖目錄
圖 1- 1 研究流程圖 ... 11 圖 2- 1 各計算元 i+1/2 處三階精度通量示意圖 ... 17 圖 2- 2 上風中央算則之近似 i+1/2 通量示意圖 ... 21 圖 2- 3 原始權重經 Mapped 函數轉換值圖 ... 26 圖 2- 4 階梯函數圖 ... 29 圖 2- 5 階梯函數的一次微分與二次微分項圖 ... 29 圖 2- 6 ㄧ維尤拉方程式駐波案例圖 ... 30 圖 2- 7 駐波案例之 、 與合值分佈圖 ... 30 圖 3- 1 計算格網示意圖 ... 32 圖 3- 2 有限體積法之數值通量示意圖 ... 33 圖 3- 3 數值通量示意圖。 ... 37 圖 3- 4 源項通量示意圖 ... 37 圖 3- 5 乾溼交界示意圖 ... 41 圖 3- 6 一維潰壩案例示意圖 ... 43 圖 4- 1 初始條件示意圖 ... 45 圖 4- 2 JSWENO 模擬水位結果圖 ... 46 圖 4- 3 ZSWENO 模擬水位結果圖... 47 圖 4- 4 初始條件設定圖 ... 48 圖 4- 5 JSWENO 模擬水位與解析解之比較圖 ... 49 圖 4- 6 JSWENO 模擬 q 與解析解之比較圖 ... 50 圖 4- 7 ZSWENO 模擬水位與解析解之比較圖... 50 圖 4- 8 ZSWENO 模擬 q 與解析解之比較圖... 51 圖 4- 9 網路概念圖 ... 54 圖 4- 10 穩態時 分佈圖 ... 55 圖 4- 11 穩態時 分佈圖 ... 56 圖 4- 12 穩態時 分佈圖 ... 56 圖 4- 13 穩態時 分佈圖 ... 56 圖 4- 14 類神經網路架構圖 ... 57 圖 5- 1 方案 A (1)模擬水位與解析解之比較圖 ... 62 圖 5- 2 方案 A (1) 模擬 q 與解析解之比較圖 ... 62 圖 5- 3 方案 A (2)模擬水位與解析解之比較圖 ... 63 圖 5- 4 方案 A (2) 模擬 q 與解析解之比較圖 ... 63 圖 5- 5 方案 B (1)模擬水位與解析解之比較圖 ... 63 圖 5- 6 方案 B (1) 模擬 q 與解析解之比較圖 ... 64 圖 5- 7 方案 B (2)模擬水位與解析解之比較圖 ... 64 圖 5- 8 方案 B (2) 模擬 q 與解析解之比較圖 ... 64圖 5- 9 方案 C (1)模擬水位與解析解之比較圖 ... 65 圖 5- 10 方案 C(1) 模擬 q 與解析解之比較圖 ... 65 圖 5- 11 方案 C (2)模擬水位與解析解之比較圖 ... 65 圖 5- 12 方案 C(2) 模擬 q 與解析解之比較圖 ... 66 圖 5- 13 代表誤差收斂圖 ... 66 圖 5- 14 簡單平滑曲線 圖 ... 69 圖 5- 15 案例 A 水位模擬結果圖 ... 70 圖 5- 16 案例 A 單寬流量模擬結果圖 ... 71 圖 5- 17 案例 B 水位模擬結果圖... 71 圖 5- 18 案例 B 單寬流量模擬結果圖... 72 圖 5- 19 潮汐波案例水深模擬結果圖 ... 73 圖 5- 20 潮汐波案例速度模擬結果圖 ... 74 圖 5- 21 AWENO 模擬水位結果圖 ... 76 圖 5- 22 水面分離處模擬水位結果圖 ... 76 圖 5- 23 底床高起處模擬水位結果圖 ... 77 圖 5- 24 水面分離案例模擬水位結果圖(t = 25 sec) ... 77 圖 5- 25 水面分離案例模擬水位結果圖(t = 45 sec) ... 78 圖 5- 26 平滑曲線處模擬水位結果圖(t = 25 sec) ... 78 圖 5- 27 上游邊界處模擬水位結果圖(t = 45 sec) ... 79 圖 5- 28 5000 格網數模擬水位結果圖(t = 15sec) ... 81 圖 5- 29 各式 WENO 算則模擬水位結果圖(t = 15sec) ... 81 圖 5- 30 水位變化處模擬水位圖(t = 15sec) ... 82 圖 5- 31 5000 格網數模擬水位結果圖(t = 60sec) ... 82 圖 5- 32 各式 WENO 算則模擬水位結果圖(t = 60sec) ... 83 圖 5- 33 水位變化處模擬水位圖(t = 60sec) ... 83
符號說明
：重建常數 Eh：水深總絕對值誤差 Eq：單寬流量總絕對值誤差 ：代表誤差 f：數值通量 ：構成平滑指示器五點絕對值之最大值 ：適合度值 g：重力加速度 ：mapped 函數 h：水深 ：最小水深容忍值 k：精度 L：基因長度 ：左特徵向量 ：空間運算子 Ma_Rate：交配機率 Ma_N：每一學習循環交配數 Mu_Rate：突變機率 max_eq：單寬流量最大絕對值誤差 max_eh：水深最大絕對值誤差 n：曼寧糙度 ：平滑指示器含二次微分項部份：平滑指示器含一次、三次微分項部份 q：單寬流量 ：右特徵向量 S：源項通量 t：時間 v：流速 W：連結加權值 ：上限值 ：下限值 z：底床高程 ：空間間距 ：修正係數 ：平滑指示器 ：一微小數值 θ：門閥值 ：特徵速度 ：權重 ：理想權
第一章 緒論
1-1 研究動機與目的
淺水波方程式(shallow-water equations)，亦稱為聖凡南方程式 (Saint-Venant equations)，在數值模擬上已被廣泛地應用於河川水流變 化、河口潮汐流況、波浪傳遞等流體運動的模擬；而淺水波方程式屬 非線性雙曲線型方程式，在震波的傳遞時可能有除峭且劇烈的變化， 其解亦可能有不連續面的產生。加權基本不震盪法(WENO - weighted essentially non-oscillatory) 在求解雙曲型方程式，對於震波等不連續面之傳遞分辨率高，且易於 推廣至高維度情況和帶有複雜源項之問題，故已成為在求解非線性雙 曲線型方程式經常選用的高解析算則之一。後續陸續發展出一些改善 原始 WENO 算則之收斂性與解析度的修正法，如守恆 WENO 法 (Vukovic and Sopta 2002)、結合 exact C (conservation)性質(Xing and Shu 2005)、Mapped WENO 法(Henrick et al. 2005)、修正平滑指示器 (MSI - modified smoothness indicator)WENO 法(Zhang and Shu 2007)、 逆擴散通量修正(anti-diffusive flux corrections)WENO 法(Xu and Shu 2005)等算則。
滑指示器，改進在強震波區域的收斂性，因此吾人將此算則應用於水 面分離與穩態流經障礙物形成水躍。經模擬結果發現在水面分離的案 例中，對於乾濕點交界上的處理有很好的效果，但對於另一穩態流經 障礙物案例有較差的模擬結果，依據文獻(Zhang and Shu 2007)中的數 值案例測詴，修正平滑指示器 WENO 算則應用在受邊界條件影響較 大的案例上會有較不好的結果，推測是由於受其上游邊界的影響，造 成在不連續面上過度收斂而影響了附近水面的連接與流量的守恆 性。 為改善此結果與增加修正平滑指示器的適用性，吾人提出一個理 論想法，以修正平滑指示器 WENO 算則的理論為基礎，設計一套與 組合平滑指示器五點值有其相關性的修正係數(此部份詳細介紹請見 4-1 節)藉以改進平滑指示器，希望在保持其收斂性的同時，亦能盡量 保持數值解的精確度與不受邊界之影響。另假設修正係數與五點值的 相關性能由一未知的函數式來代表，故本文應用類神經網路類神經網 路(ANN - artificial neural network)的黑盒分析概念，結合遺傳演算法 (GA - genetic algorithm)的優選特性，以穩態流經障礙物案例為學習案 例，建構類神經網路以擬合此假設的函數式，之後應用建構完成之類 神經網路於詴驗案例，作為驗證結果並與各 WENO 算則模擬之結果
1-2 文獻回顧
1-2-1 類神經網路與遺傳演算法
類神經網路是一種以電腦的邏輯運算模擬人類腦神經細胞網路 的學問。McCulloch and Pitts (1943)提出神經元數學模型，利用數學模 式對人類神經系統中的神經元運作建立數學模型，把生物神經單元視 同一個邏輯決策單元；Werbos (1974)在其博士論文中提出了隱藏層的 學習演算法，這是首次有人提出倒傳遞網路(back-propagation network) 的概念。之後 Parker (1985)再次提出倒傳遞網路，同年 Rumelheart et al. (1985)發表了一篇倒傳遞網路的文章，倒傳遞類神經網路才使其廣為 人知，而倒傳遞類神經網路模式是目前類神經網路學習模式中最具代 表性，應用最普遍的模式。經過近數十年的努力，類神經網路的理論 基礎逐漸完善，現已被視為極有效的非線性模型建構工具，因其理論 單純、結構單純，容易建立所需的特定模式，因此大量被應用於各領 域上。 遺傳演算法的理論基礎可回溯自 1859 年達爾文(Charles Darwin) 的「物種演化」 （On the Origin of Species by Means of Nature Selection） 書中的「物競天擇，適者生存」的演化及淘汰觀念。將這種自然界的 選 擇 方 法系 統化並 發 展 為可 用之模 式 ， 最早 是由密 西 根 大學 的
Holland (1975)發展出遺傳演算法搜尋技術的基本架構，並且由其學生 Goldberg (1989)成功地運用在工程問題上。而利用此法來搜索演化類 神經網路的架構，此即為混合遺傳演算法 (HGA - hybrid genetic algorithm)，在許多研究上，相對於倒傳遞類神經網路有不易陷入局 部最小值的優點，在預測上亦能有較佳的穩定性與誤差。
類神經網路與遺傳演算法在水利方面的應用上，國內外的相關研 究自然也不少。諸如應用於降雨與逕流，或結合地形條件做預測模擬 (French 1992, Sajikumar and Thandaveswara 1992, Lorrai and Sechi
1995, Loke et al. 1997 ,Zhang and Rao 2000)，(陳昶憲等 1996、孫建平
1997)、地下水位預報(Shigidi and Grcia 2003, Daliakopoulos et al. 2005,)、水質分析(郭異銘 1999、范正成等 2006)、土石流災害預警(李 心平和張斐章 1995、張東炯 2000、曾國源 2003)與水庫系統操作(Jain et al. 1999, Chaves and Chang 2008)等，應用的範圍可說是相當地廣 泛。
1-2-2 高解析算則的發展
由於雙曲線方程式存在著不連續解，這些不連續解在使用許多傳 統數值方法模擬時，易產生 spurious oscillations 的震盪問題，為了處 理震盪的問題，Harten (1983) 以 Van Leer (1974, 1979) 提出的通量修
函數(limiter)來限制守恆量或通量的梯度變化，來減低數值傳遞產生 的震盪，而後許多高階(指超過一階精度) TVD 算則被陸續發展出來， 廣泛的介紹與文獻回顧可參考 Roe (1986)、Yee (1989)及 LeVeque (1990)，高階精度的 TVD 算則發展提升了雙曲線型方程式數值解的 品質，也正式開啟了高解析算則的時代。
Harten et al. (1987)提出 ENO (essentially non-oscillatory)算則，該 算則的特點是依據所定義的平滑指示器(indicator of smoothness)，選 取出最為平滑代表的計算元(stencil)進行差分，得到均勻高階精確度 且具有基本不震盪的特性，計算元是指求某一點差分值時，其所包含 之相關點的集合，r 階精度 ENO 算則採用 r 組計算元做為差分選擇的 候選，例如計算 i+1/2 介面的數值通量正的部份，i 為自身的點位， 以三階精度為例，使用三點差分而有三組計算元可選擇，分別為： [i-2, i-1, i] ﹑[i-1, i, i+1] ﹑[i, i+1, i+2]。此算則是屬於守恆變數型態 (variables version)，其所使用的重建(reconstruction)計算是擴充 Harten and Osher (1987)所發展的二階重建模式，此外 Shu and Osher (1988, 1989)進一步發展逐點式(pointwise)的 ENO 算則來取代原先單元平均 的計算，在二維或是更高維度的計算上能獲得較佳的效率。ENO 算 則特別適用於包含震波及複雜的平滑流場結構，如紊流與震波交互作 用、渦漩流與震波交互作用等流場。ENO 算則已被應用及擴展在許
多方面，如在選取計算元的過程加入偏量(bias)以加強穩定及精度 (Fatemi et al. 1991, Shu 1990)、利用 sub-cell resolution 使接觸面更為 除峭(Harten 1989) 、多維的三角網格分析 (Abgrall 1994)等。
Liu et al. (1994)針對 ENO 算則的缺點，提出了加權型基本不振 盪(WENO)算則，取代 ENO 算則僅選取一組計算元來近似數值通量 的作法，以定義的平滑指示器來決定各計算元的代表平滑程度，依此 分配每組計算元一個加權係數，將所有採用的計算元全部組合起來。 而有關配置的原則，在平滑區儘量使其最佳化藉以達到高階準確度的 要求，而在不連續面附近則令含有不連續面之計算元的加權係數盡量 減小，最後再以外凸組合(convex combination)的方式來近似數值通量， 在單調平滑的區域，r 階精度 ENO 算則經由組合後在能提高為(r + 1) 階精度的 WENO 算則，然而無法達到最理想(2r - 1)的精度。之後 Jiang & Shu (1996)以上述之 WENO 算則為基礎，重新定義了新的平滑指示 器，發展了 WENO5 算則，使用了三階精度的 ENO 算則組合之 WENO 算則，在特定條件下能保持五階的精度，震波的傳遞計算也依然能維 持基本不震盪的原則。自 Jiang and Shu 在 1996 年發表新式 WENO 算則以來，此篇文章在 SCI 資料庫中廣被引用，其應用包括：算則發 展有關(Gerritsen and Olsson 1998, Liu and Osher 1998, Choi and Liu 1998, Suresh and Huynh 1997, Jiang and Yu 1998, Jiang and
Tadmor 1998) ， 應 用 在 二 維 非 結 構 性 網 格 (Friedrich 1998, Ollivier-Gooch 1997) ，應用在不可壓縮流(Minion and Brown 1997, Yang et al. 1998, Chen et al. 1999)，應用在可壓縮流(Yang et al. 1998, Yang et al. 1992, Hsieh et al. 2008)，與 Shallow-water 有關(Kuo and Polvani 1999, Liska and Wendroff 1999)等。
WENO 算則由於加權係數的值是由流場的物理量運算後所求得， 不需如 ENO 算則用到邏輯判斷平滑度來選擇使用之計算元，因此電 腦程式的計算方式可以向量化處理，而數值通量不再因為計算元的選 取而有跳動現象，相對地也改善了收斂的效率，WENO 算則除了保持 ENO 算則的均勻高階的優點，同時改善了 ENO 算則的收斂性。 近年來有一些研究針對原始 WENO 算則，提高其收斂性與解析 的精度。Despres and Lagoutiere (2001)首先提出一階限制函數下風算 則(limited downwind scheme)，可有效的避免流場中不連續接觸面的 糢糊化(smearing)並同時保有非線性的穩定性，之後 Bouchut (2004) 修正此算則使其滿足熵條件且提出一簡化的顯式逆擴散通量公式，在 前兩項研究的基礎下，Xu and Shu (2005)提出了高階有限差分逆擴散 通量修正之 WENO 算則，此一算則不僅在流場平滑區可維持高階精 度外，在流場不連續區域也可保持不震盪的特性與不連續接觸面除峭 的高解析度；Henrick et al. (2005)提出了 Mapped WENO 算則，此一
算則改善了原始 WENO 算則在臨界點處的精確度；Zhang and Shu (2007)針對高階 WENO 算則提出了修正平滑指示器的觀念，依據數 值詴驗的測詴結果，顯示在邊界無震波穿越的影響下，一維與二維的 高 強 度 震波 問題使 用 此 一修 正算則 可 改 進其 收斂性 至 機 械零 點 (machine zero)。
1-2-3 守恆定律的探討
長久以來，許多文章探討數學模型之守恆律(conservation law)， 即流體在某一控制體積之物理量需守恆。Berm dez and V´azquez (1994)於推求淺水波方程數學差分過程中，使用 Q-scheme 與源項採 上風法(upwinding)，考慮質量及動量守恆之觀念來計算通量梯度與源 項，故當求解問題為穩態靜止時，所得解皆能符合守恆律，此稱之為 exact C (conservation)性質，Hubbard and Garcia-Navarro (2000)則延伸 其數值方法發展出二階 MUSCL 有限體積法。其它另有針對源項及數 值通量作守恆性質探討 ，如：Greenberg and LeRoux (1996)提出 well-balanced 算則、LeVeque (1998)發展出準穩態(quasi-steady)波傳的 運算法則、Jenny and Muller (1998)針對源項採用 Rankine-Hugoniot 黎 曼解來做修正、Chinnayya and LeRoux (1999)發展出一套新的黎曼求 解方式，運用於一維淺水波方程式、Botchorishvilli et al. (2000)結合或 Roe 型態的分離與上風法，建構出一個補抓穩定型態的方法、Zhou et al. (2001)在源項的處理上，有別於一般所提倡的上風法，而引用 surface gradient method (SGM)的概念等。
Gosse (2001) 則 進 一 步 提 出 加 權 基 本 不 震 盪 法 之 守 恆 法 則 ， Vukovic and Sopta (2002)則利用 Goss 所得之法則，求解一維淺水波方 程，Crnjaric-Zic et al. (2004)推導求解守恆型式一維輸砂方程式，兩者 運算法則之概念皆為加權基本不震盪法，輔以高階數值通量、有限體 積法為基礎，以減少數值通量進出控制體積之誤差，以達質量及動量 守恆之要求。
1-3 研究方法與流程
本研究在數模計算上主要分水理計算與類神經網路兩大部份，在 水理計算上引用守恆 WENO 法(Vukovic & Sopta 2002)，搭配有限體 積 法 求 解一 維淺水 波 方 程式 ， 類神 經 網 路架 構採取 前 饋 式網 路 (forward networks)。圖 1- 1 為本研究流程圖，首先將類神經網路各單 元之網路加權值與門限值，以二進位作為基因組的編碼，並產生初始 基因組群集。在學習過程中，各基因組透過類神經網路與運算，可以 決定在 1-1 節中所定義之修正係數，將之代入平滑指示器的運算以決 定 WENO 算則的權重，並回歸水理的計算中。將計算後各基因組所得之解與解析解做誤差比較，應用遺傳演算法的優選概念，誤差小的 基因組給予較佳的適應度(在複製時擁有較大的機率被選取)，誤差大 則對應較差的適應度，而後經過複製選取、基因交配與基因突變等機 制以產生下一子代的基因組，以上整個流程稱為一次學習循環；在經 過數千或數萬次學習循環後，理論可得一組基因組以有效架構類神經 網路，使計算結果能更接近學習案例之解析解，之後使用架構之類神 經網路運用於其他案例，作為驗證結果並與各式 WENO 算則所模擬 之結果做比較分析。
基因編碼 產生初始基因組群集 水理計算 基因複製、交配與突變 Yes 是否到達最後 模擬時間 No 是否基因組 群集全運算 完畢 選用基因組群集中 未模擬過的一組基因組 Ｎo 是否已達到最大 學習循環次數 NO 產生新一代的基因組群集 學習結束 Yes 計算誤差與適合度 Yes 透過類神經網路 決定修正係數 代入平滑指示器 決定WENO之權重 類神經網路計算 圖 1- 1 研究流程圖
第二章 理論基礎
有關類神經網路與遺傳演算法的介紹請見本文附錄 A 與附錄 B。 本章針對前述基本不震盪相關算則，逐一介紹 ENO 算則、WENO 算 則、差分式 WENO 算則，並以差分式 WENO 算則為基礎分別說明 Mapped WENO 法、修正平滑指示器 WENO 法等算則。
2-1 ENO 算則
Harten et al. (1983)提出了 ENO (essentially non-oscillatory)算則， ENO 算則是均勻高階，且基本不震盪的解析包含不連續的流場。ENO 算則設計的特點在尋求最平滑的計算元進行高階差分，故可避免在極 值點附近降階以及在不連續點附近產生非物理性的大震盪現象。
2-1-1 一維數值通量的建構與近似模擬
假設一維格點分佈如下： (2-1)定義單元(cell)，單元中心(cell centers)，單元大小(cell size)如下：
， ， ，
， (2-2)
給定單元內的平均值函數為 f(x)：
求解一多項式 ，其最多為 k-1 階，使其滿足 ， ， (2-4) 求解上述問題首先定義原函數(primitive function) F(x) 如下： ∞ (2-5) 上式積分下限依需要可改成某個固定值，利用(2-3)式的單元平均 值的定義，可計算如下： ∞ ∞ ∞ (2-6) 如此，若已知單元平均值 ，則可得到原函數在每個單元邊界 點(cell boundary)的值；給定單元 的位置及精度 k，我們選定計算 元(stencil) S(i) 包括：左邊 r 個單元，右邊 s 個單元及本身 ，其中 r+s+1=k： (2-7) 若定義 P(x)滿足 ， (2-8) P(x)的微分表示為 p(x) (2-9) 而在計算元內每一個單元平均值皆可表示為
∞ ∞ ， (2-10) 上式顯示 p(x) 就是所尋找的函數。 一般而言，前述重建問題大多用來計算單元邊界點的重建量，再 以此重建量計算單元邊界位置的通量，此計算如下： ， ， (2-11) 又由於前述由單元平均量 映射至單元邊界點的重建量 及 的過程是線性的，故存在常係數 及 使得 ， (2-12) 其 中 及 與計算元內包含的左邊單元數 r、內差階數 k 與 單元大小 有關，與函數值 f 無關。上式在相同位置 得到 ± 兩個不同的重建量，分別來自於單元 與 。前述兩重建量相同， 依據對稱的概念得到下列關係式： (2-13) 總結以上的推導，另一組常數 重新構造單元邊界上 的
(2-14) 以一維簡單雙曲線方程式為例 (2-15) 假設為均勻格網，則該方程式在空間離散為 (2-16) 其中 與 為數值通量。若 (2-17) 若存在函數 h(x) 使得 (2-18) 則 可寫作 (2-19) 故若 (2-20) 為 k 階精度算則。
2-1-2 各計算元的數值通量
參照(2-20)式，各計算元之 以三階精度近似函數 h(x)，r 為計算元內，以本身 為基準，所採用的左邊單元數目，以多項式 表示為 ℎ (2-21) 將上式代入(2-18)式，可推得各計算元之 、 、 ；依據採用 不同的計算元，所得到不同的三階精度多項式表示如下： (2-22) 上述各多項式近似位於 i+1/2 的通量 ，如圖 2- 1 所示可得 (2-23)圖 2- 1 各計算元 i+1/2 處三階精度通量示意圖 求通量 僅需將所有計算元內單元皆向左移動一格計算，將 (3-23)式與 做泰勒展開可表示為 (2-24) 通式為 ， (2-25) 此三階的誤差 在之後 WENO 算則中，定義平滑指示器與決 定權重的部份佔了極為重要的角色。此外空間離散項亦會由於 與 選擇計算元的不同，而有不同的精度：
， ， (2-26) (2-23)式為(2-14)式之精度 k = 3，各計算元的展開式子，而依據上 述的觀念，同理可得各精度與各選取計算元之重建常數 ，若網格 點為均勻分佈 Shu (1997) ， 如表 2- 1： 表 2- 1 重建常數 表(k=1-4) k r j=0 j=1 j=2 j=3 1 -1 1 0 1 2 -1 3/2 -1/2 0 1/2 1/2 1 -1/2 3/2 3 -1 11/6 -7/6 1/3 0 1/3 5/6 -1/6 1 -1/6 5/6 1/3 2 1/3 -7/6 11/6 4 -1 25/12 -23/12 13/12 -1/4 0 1/4 13/12 -5/12 1/12 1 -1/12 7/12 7/12 -1/12 2 1/12 -5/12 13/12 1/4 3 -1/4 13/12 -23/12 25/12
2-1-3 推展與通量分離
傳統的二階算則，如 Lax-Wendroff 算則及 Beam-Warming 算則， 在不連續點附近產生非物理性的震盪，主要是採用固定的計算元進行差分，ENO 算則設計則選取最平滑計算元的進行差分，其計算元的 選取是隨著流場結構改變，這些改變主要是使數值差分時避開不連續 點。
Shu and Osher (1988, 1989)對原始之 ENO 算則做了些修正，在此 新的 ENO 算則中，其建構數值通量的方式避免採用單元平均值來重 建函數，而改以直接對通量做類似之基本不振盪差值來求取數值通量。 另外他們也捨棄了原 ENO 算則中的 Lax-Wendroff 形式時間離散法， 改 用 時 間 與 空 間 分 開 考 慮 之 半 離 散 化 系 統 ， 另 行 發 展 了 TVD Runge-Kutta 形式之時間離散法，如此改善了原 ENO 算則不易寫成 程式及不易推展到多維問題之缺點，也加強了 ENO 算則的效率。 在實際應用上，除前述的 ENO 差分方式外，一般會加入通量分 離(flux splitting)的概念。物理通量正、負的分離方式，有許多分離可 供選擇，數值通量則定義如下： (2-27)
K 為一分離函數，有 Godunov flux 及 Engquist-Osher flux、Lax
-Friedrich 等方法可供選擇，在此介紹總體或局部 Lax-Friedrich 之通 量分離法來求解流場，以(2-15)式為例，總體 Lax-Friedrich 之通量分 離法如下:
(2-28) 其中 ，代表所有計算網格內的最大特徵值。 局部 Lax-Friedrich 之通量分離法為： (2-29) 其中 ，代表局部計算網格內的最大特徵值。
2-2 WENO 算則
WENO 算則的基本理念：取代 ENO 算則僅選取一組計算元來 近似數值通量的作法，而將所有需要研判平滑度的計算元全部組合起 來，並賦予每組計算元一個加權係數，最後再以外凸組合(convex combination)的方式來近似數值通量： (2-30) 為了穩定性與一致性，權重 頇滿足 ， (2-31)2-2-1 上風中央算則與理想權重
與 2-1-2 節之多項式的推導相同，上風中央算則(Upstream central scheme)源自於實際通量 h(x)五階精度的多項式近似值：ℎ 將上式代入(2-18)式，可展開積分式為 ) (2-32) 利用計算元內之五個單元 ， 、 、 、 、 ，可得五組方程式聯立求解 ，再近似 位於 i+1/2 的通量 ，如圖 2- 2 所示可得 (2-33) 圖 2- 2 上風中央算則之近似 i+1/2 通量示意圖 同理，在 處，通量 為
(2-34) 因此將(2-33)式與(2-34)式代入(2-17)式，可得五階精度空間離散 的上風中央算則： (2-35) 而在平滑區域中，數值通量近似解 給定一組權重以符合上 風中央算則所得的線性組合，藉以達到高階準確度的要求。因此將式 (2-23)以式(2-30)的形式組合成式(2-33)，整合得下式(2-36)： (2-36) 因此可得到一組常數： 、 、 ，此又被 稱為理想權重(ideal weights)，在平滑區能滿足五階的精度： (2-37)
2-2-2 平滑指示器
如果區域內有不連續面的存在，在組合數值通量中，將降低含 有不連續面的計算元的權重值，以減少數值震盪的發生與影響；而 如何去衡量計算元內的平滑程度，來做權重值的調整與修正，在此利用定義之平滑指示器(Jiang and Shu 1996)： (2-38) 將(2-21)式代入，可得
(2-39) 再代入(2-22)式之 、 ，可得各計算元代表之平滑指示器為 (2-40) 而各計算元權重的修正則為 ， 其中 (2-41) 為一微小數值，以限制 值有其最大上限值，而不會因 時， 因為無限大值超過程式所能定義數值的最大邊界， 而造成程式的運算中止。
2-2-3 WENO5 算則
與 ENO 算則相同，在此導入通量分離的概念將物理通量分離成 正、負兩部分，定義出以三組計算元所組合而成的 WENO5 算則 (Jiang and Shu 1996)。以下為 WENO5-LLF-S 算則的說明：(2-42) 其中權重 計算如下 ， ， ， ， (2-43) 平滑指示器 定義為 (2-44) 數值通量負的部分 (2-45) 其中權重 計算如下 ， ， ， ， (2-46) 平滑指示器 定義為 (2-47)
2-3 差分式 WENO 算則
差分式算則將數值通量分為一階及高階通量，各種滿足熵條件 的一階算則多可直接應用。差分式算則與 Shu 等所發展的算則主要 差異是以通量差值進行分離，原型式則是以通量進行分離。差分式 算則如下： (2-48) 其中 為一階數值通量， 為高階數值通量。有許多一階算 則可供選擇為一階數值通量，如前節所介紹之 Godunov flux 及 Engquist-Osher flux、Roe flux、Lax-Friedrich flux 等算則。之後在第 四章會以一維淺水波方程式為例，對差分式 WENO 算則的計算流程， 做更詳細的介紹。2-4 Mapped WENO 算則
Henrick et al. (2005)觀察到了原始 WENO 算則中計算非線性權重 所使用的平滑指示器在某些特定的臨界點時可能會失去精確性，因此 他們提出了 Mapped WENO 算則，此一算則改善了原始 WENO 算則 在臨界點處的精確度，Mapped WENO 算則主要以 Mapped 函數來計 算原始 WENO 算則中非線性權重：
其中 ，為原始 WENO 算則計算所得的權重， 則為 上述之理想權重值。此一 Mapped 函數為單調漸增函數，具有限斜率 及下列特性： ， ， ， ， (2-50) 調整後的權重可由下式獲得：
，
(2-51) 而經 Mapped 函數轉換所得之
值
與原始值的相關位置比較如 圖 2- 3 所示： 圖 2- 3 原始權重經 Mapped 函數轉換值圖 (摘自 Henrick et al. 2005)根據 Henrick et al. (2005)分析的結果，此一算則不但可以改善平 滑區域的精度，同時受到微小參數值 ε的影響也大為減少，收斂性 則改善許多。
2-5 修正平滑指示器 WENO 算則
非線性權重與平滑指示器在震波附近之解析度上扮演很重要的 角色，同時兩者均直接影響震波前後之震盪與收斂性。Zhang and Shu (2007) 針 對 高 階 WENO5 算 則 提 出 了 修 正 平 滑 指 示 器 (modified smoothness indicator)的觀念，依據他們測詴的結果顯示，在不受邊界 條件之影響下的一維與二維強震波無黏性穩態流場問題中，此一修正 後的算則可以改進其收斂性。 將(3-40)式採泰勒級數展開： (2-52) 可將平滑指示器分為兩部份： ， 含有一次微分與三次微分項 ， 含二次微分項 (2-53) 在極劇的高強度震波傳遞中，以圖 2- 4 的階梯函數為例，圖 2- 5
為階梯函數的一次微分與二次微分項的大小分佈，可發現越靠近震波 的中心點 C 處，一次微分項越大而二次微分項則越小趨近於零，因此 若加入二次微分項做為平滑度的判別，對於震波中心兩鄰近點 B、D 會增大其平滑度，而影響震波傳遞在計算上的收歛性。
在原始平滑指示器中， 項的大小主決於一次微分項，而 項
則為二次微分項，綜合以上的結論。Zhang and Shu (2007)建議在穩定
的震波區域捨棄 項使用修正平滑指示器。原始 WENO5 算則之平 滑指示器(2-44)式與(2-47)式，經修正後表示如下： ， ， (2-54) ， 同理，修正平滑指示器亦適用於差分式 WENO5 算則。
圖 2- 6 為 Zhang and Shu (2007)中之ㄧ維駐波案例，所使用統御 方程式為一維尤拉方程式，數值求解使用 WENO 算則，其流程與計 算類似於一維淺水波方程式，此案例相關的設定參數與起始條件詳
細請見以上文獻。圖 2- 7 為此案例震波中心與鄰近處(2-53)式之 、
與 值，可明顯發現在捨去 項後，計算元之間的相對
圖 2- 4 階梯函數圖
圖 2- 6 ㄧ維尤拉方程式駐波案例圖
圖 2- 7 駐波案例之 、 與合值分佈圖
第三章 WENO 運用於一維淺水波方程式
本章說明 WENO 算則應用於一維淺水波方程式，數值方法時間 離散採用三階時間準度之 TVD Runge-Kutta 積分法，空間離散則採 用差分式的守恆形式 WENO 算則(Vukovic and Sopta 2002)，搭配有限 體積法求解，並延續 2-3 節對差分式數值通量作更進一步的說明。3-1 統御方程式
淺 水 波 方 程 式 ， 源 自 於 納 維 - 斯 托 克 斯 方 程 式 (navier-stokes equations)採水深積分的形式，並定義以下幾點以作簡化：(1)水平方 向的尺度遠大於垂直方向的尺度、(2)不可壓縮流、(3)忽略風力與柯 氏力的影響、(4)垂直壓力梯度採靜水壓力分佈、(5)在水平方向的底 床坡度變化極小、(6)在垂直方向的速度分佈為定值。 因此一維淺水波方程式守恆形式(conservation form )表示如下： (3-1) ， ， (3-2) 其中 h 為水深(m)， q 為單寬流量(cms)，v = q/h 為流速(m/s)，z 為 底床高程(m)，n 為曼寧糙度，g 為重力加速度(m/s2 )。另定義 Jacobian 矩陣 A： (3-3)上式所對應之特徵值有二： 、 (3-4) Jacobian 矩陣 A 的左、右特徵向量依上下式分為 、 ，上 標 p 表示為運算矩陣中第 p 個位置，下標(i+1/2)意指格網 i 與i+1 之 交界面位置，而交界面位置之特徵速度及特徵向量以算數平均求解， 圖 3- 1 為計算格網示意圖。 圖 3- 1 計算格網示意圖 將(3-1)式中對流項移至等號右側，再以空間運算子 代表 對流項與源項兩項之和，表示如下： (3-5) (3-6)
3-2 數值方法
3-2-1 時間離散
利用三階 TVD-Runge-Kutta 積分法，考慮半離散化(semi-descrete) 系統，將(3-5)式表為時間項之離散式後，其積分形式如下： i i+1 i+1/2(3-7)
3-2-2 空間離散
基於每個網格單元的局部通量平衡，本文採用空間離散與時間離 散分開處理的半離散有限體積法(finite volume method)。(3-5)式之空 間運算子部分 ，以有限體積法數值通量之格式表示： (3-8) 式中 為源項通量， 、 為數值通量，計算採用五階的差分 式守恆形式 WENO 算則。(3-8)式說明格網左右交界面淨通量之和， 加上該格網之源項通量即為控制體積內物理量之變化，圖 3- 2 為數值 通量示意圖。 圖 3- 2 有限體積法之數值通量示意圖 Si f -i-1/2 f +i-1/2 f -i+1/2 f +i+1/2 χi-1/2 χ χi+1/2 i χi+1 χi-1 Fi-1/2 Fi+1/2
3-2-3 守恆形式差分式 WENO 算則
引用 Vukovic and Sopta (2002)之文獻，說明差分式守恆形式 WENO 算則之理論。依據(2-48)式可知差分型式包含一階與高階數值 通量的計算，以下分別介紹 LLF 算則與 RF 算則。 a. Lax-Friedrichs flux(LLF)算則： 一階通量的部份採 (3-9) ，為全計算域最大特徵速度之絕對值 高階通量則定義為 (3-10) 組合式(3-9)與式(3-10)可得完整之 LLF flux 算則： (3-11) b. Roe flux(RF)算則： 一階通量的部份採
， ， (3-12) 高階通量則為： ， ， (3-13) 組合式(3-12)與式(3-13) 得完整之 LLF flux 算則： ， ， (3-14) 在數值模擬時因流場特徵速度傳遞方向不同，而選擇不同的數 值通量分離方法，分別為 LLF flux 及 RF flux formulations，以下整 理各種情況為：
(a) (LLF flux formulations)
(3-15) (b) (RF flux formulations) (3-16) 上兩式中， 、 兩項代表高階之正負數值通量，計算以 WENO 算則格式表示如下： (3-17)
上式 則為：
(a) (LLF flux formulations）
(3-18) (b) (RF flux formulations) (3-19) 其中，下標 。 (3-15)式說明兩格網間特徵速度之傳遞方向互為反方向時，代表 數值通量在 i+1/2 位置有進有出，所以高階數值通量應同時考慮 、 兩項，如圖 3- 3(a)所示；反之 (3-16)式則說明若兩格網 間特徵速度之傳遞方向皆為同方向時，代表數值通量在 i+1/2 位置 只 往 單 一 方 向 傳 遞 ， 因 此 高 階 數 值 通 量 只 需 依 傳 遞 方 向 選 擇 或 ，示意圖如圖 3- 3(b)及圖 3- 3(c)。 (a) χi+1
Fi+1/2=0.5(fi+1+αui+1)li+1/2+P-i+1/2 Fi+1/2=0.5(fi-αui)li+1/2+P+i+1/2
χi+1/2 λi+1 (p) λi (p) χi LLF flux
(b) (c) 圖 3- 3 數值通量示意圖 (a) ， (b) ， (c) 源項部分，同理使用通量分離的概念，將源項影響分為左、右兩 部分如所示。 (3-20) 圖 3- 4 源項通量示意圖 χi+1 Fi+1/2=fi·li+1/2+P+i+1/2
χi+1/2 λi+1 (p) λi(p) χi χi+1 Fi+1/2=fi+1 ·li+1/2+P-i+1/2
χi+1/2 λi+1(p) λi(p) χi χi+1 Si-1/2,R χi+1/2 χi-1 χi-1/2 χi Si+1/2,L RF flux RF flux
同理依據特徵速度傳遞方向不同，分別選擇 LLF flux 及 RF flux
的數值通量分離方法，將 、 分離表示如下：
(a) (LLF flux formulations)
(3-22) (b) (RF flux formulations) (3-23) 以上 G 為某定義函數， 及 計算方法可參考 (3-17) ~ (3-19)式，以 WENO 算則格式表示如下： (3-24) 上式 則為：
(a) (LLF flux formulations)
(3-25) (b) (RF flux formulations) (3-26)
其中 ，G、 Z 函數則依守恆律定義為： (3-27) (3-28) (3-27)式與(3-28)式中， 、 與 、 為對應其點位 1 與 2 之水深與底床高。 在原始的差分式 WENO 算則中，對流項(3-17)式中的權重依據 (3-18)與(3-19)式的 通量計算決定，為 形式；源項(3-24) 式的權重，則是以(3-25)式與(3-26)式的 通量為依據，為 形式。
而差分式守恆形式 WENO 算則依據結合 exact C (conservation) 性質，水體靜止時，則必滿足： (3-29) 再結合了對流項與源項所對應的通量，權重皆以( )作為 新輸入通量計算決定，為 ( )的形式，因此重新定義(3-17) 式與(3-24)式： (3-30) (3-31)
3-3 乾溼交界處理
淺水波方程式於數值計算中常含有源項，如應用於實際不規則地 形之淹水模擬時，對乾、濕交界之處理不可避免。當水流遇到高出水 面之地形時，水深變為零，在計算特徵速度時，分母有一水深 h 項， 計算到水深為零處，會產生發散問題。 本文參考 George (2004)之推導，其概念為先求解濕鋒之特徵速度， 再利用積分曲線推求乾鋒之特徵速度，如此可避免乾鋒水深為零處產 生之數值發散問題。 此外定義一最小水深容忍值 ， 之值可依模擬案例之尺度加 以選擇，由前一時刻計算所得之水深 h 與 比較，以判斷計算點之 乾溼狀況。若水深 h 大於 則該點屬於濕床情況，無頇任何處理即 可進行數值計算；當計算遇到水深為零或小於 時，則表示該點為 乾床狀態，需推求乾鋒之特徵速度。 當溼鋒由左至右傳遞至乾床，如圖 3- 5(a)所示，位置 i+1/2 之特 徵速度為： (3-32) 反之溼鋒由右至左傳遞至乾床，如圖 3- 5(b)所示時，位置 i+1/2 之特徵速度為：(3-33) 圖 3- 5 乾溼交界示意圖 a 為溼鋒由左至右傳遞，a1 為水面線，a2 為特徵速度示意圖。 b 為溼鋒由右至左傳遞，b1 為水面線，b2 為特徵速度示意圖。
3-4 稀疏波
當控制體積內各點之特徵曲線兩兩相互遠離，且無交會之可能， 該現象稱稀疏波(rarefaction wave)。 以潰壩案件為例，當潰壩發生時，在潰壩水體前端會產生一個 震波(compression wave)，並迅速形成一除直不連續面向前移動，同時 於潰壩水體後方產生一稀疏波往後方移動，如觀察特徵曲線圖則可發 現稀疏波會在特徵曲線圖上形成一扇形區域，如圖 3- 6 所示；模式模 χi+1 χi+1/2 χi λ2 i+1/2 hRHS h=0 λ1 i+1/2 χi+1 χi+1/2 χi λ2 i+1/2 hLHS h=0 λ1 i+1/2 h t a1 a2 b1 b2擬之數值解可分為兩組弱解，分為稀疏波解及震波解，兩解皆能滿足 方程式，這個不唯一性乃守恆律方程式的一個缺陷。因此本文參考 Einfeldt (1988)之熵修正概念，將(3-4)式之特徵速度修正如下： (3-34) 其中 、 為 Roe Speed，藉由此修正式對震波解附加條件， 以篩選出具物理意義之解。
圖 3- 6 一維潰壩案例示意圖(上：起始水深、中：特徵速度 曲線群、
下：特徵速度 曲線群；實線扇形為稀疏波範圍，虛線為震波位置。)
第四章 模式建置
由 1-3 節研究方法與流程，可知本文模式的架構主要可分為水理 與類神經網路兩大部份。本章除了定義水理與類神經網路的架構設定 外，並針對自訂之類神經網路，所參考的理論、輸入參數與學習案例 的選擇等，作更詳細的說明。4-1 理論依據
延續 2-5 節的介紹，可知在含高強度震波的不連續面處，捨去原 始平滑指示器中含二次微分項的部份，可增大含不連續面的計算元， 相對於鄰近計算元不平滑的程度，在線性組合數值通量時，能降低該 計算元所代表的權重值，使在震波傳遞的計算上提高其收斂性。以下 開始，定義原始 WENO 算則代號為 JSWENO (Jiang and Shu 1996)， 修正平滑指示器 WENO 算則為 ZSWENO (Zhang and Shu 2007)。因此利用 JSWENO 與 ZSWENO 算則分別對兩種含有不連續面解 的案例模擬並比較結果： (a) 水面分離案例：(仿 Bouchut 2004) 底床幾何形狀為： (4-1)
初始水流條件： (4-2) (4-3) 水流邊界條件為 上下游邊界 (4-4) 綜合以上條件，圖 4- 1 為初始條件設定圖。 Distance(m) 0 500 1000 1500 2000 2500 W ater Le vel( m) 0 2 4 6 8 10 12 14 Bed Level
Initial water surface
圖 4- 1 初始條件示意圖，虛線左側初始流量為-350cms，右側初 始流量為 350cms，水流各往左、右流動。 此案例為數模之假設案例，目的為測詴模式處理乾濕交界問題的 能力。圖 4- 2 為 JSWENO 算則採計格網數 400 格( 6.25 m)與 8000 格( 0.3125 m)，可蘭數採 0.75， m，於 t = 5.0 sec 之水位模擬結果，以下將 JSWENO 8000 格網點之模擬結果假定為解 析解以作為討論。由圖可看出 JSWENO 算則在水面分離的中心處，
會有些許的水面震盪發生，推測是由於在計算乾濕床交界近處，若其 最小水深 取較小值時，水深在稍大於 時，擁有高單寬流量 q， 一旦水深判斷小於 ，視為乾床其水深與單寬流量都判斷為零值， 而使得數值通量變化劇烈，造成通量傳遞計算上的難以收斂，需增加 格網點以降低彼此間通量的差異或加大 ，以減低水面震盪的高度 與振幅。圖 4- 3 為同等條件下，ZSWENO 採 400 網格數與假定解析 解的水位比較圖，可看出透過強化收斂性之修正平滑指示器，可以使 在水面分離中心處無水面震盪的產生，但也因收斂強化的結果而使水 位點位的連接上與 JSWENO 算則的模擬結果有些許的差異。 Distance(m) 0 500 1000 1500 2000 2500 W ater Le vel( m) 0 2 4 6 8 10 12 14 Bed Level JSWENO(400 grid) JSWENO(8000 grid) 圖 4- 2 JSWENO 模擬水位結果圖 ( m，t = 5.0 sec)
Distance(m) 0 500 1000 1500 2000 2500 W ater Le v el( m) 0 2 4 6 8 10 12 14 Bed Level ZSWENO(400 grid) JSWENO(8000 grid) 圖 4- 3 ZSWENO 模擬水位結果圖 ( m，t = 5.0 sec) (b) 穩態流流經障礙物之案例：(V´azquez-Cendon 1999) 採用網格數 200，模擬間距 0.125m，重力加速度為 9.81 m/s2 ，可蘭數採 0.75，摩擦損失忽略不計，總模擬時間 200 秒， 並假設此時到達穩態。 底床幾何形狀定義如下： ， ， (4-4) 水流起始條件為 起始水深(m)： 。 (4-5)
起始單寬流量(m2/s)： 。 (4-6) 水流邊界條件為 上游邊界 ；下游邊界 。 (4-7) 綜合以上條件，圖 4- 4 為初始條件設定圖。 Distance(m) 0 5 10 15 20 25 Value -0.1 0.0 0.1 0.2 0.3 0.4 0.5 Bed Level Water Level
Discharge per unit width
圖 4- 4 初始條件設定圖
表 4- 1 為 JSWENO 與 ZSWENO 算則結果對於解析解誤差表， 圖 4- 5 至圖 4- 8 分別為 JSWENO 與 ZSWENO 模擬水位與單寬流 量和解析解比較圖，其解析解可見 Goutal and Maurel (1997)一 書。
由以上的結果可看出，此案例 ZSWENO 相對於 JSWENO 算 則，雖然最大單寬流量誤差有一定程度的減小，但也因此降低了
不連續面附近的單寬流量守恆性，此外水深也在水躍前後附近點 的連接上，與解析解有較大的差異。 表 4- 1 JSWENO 與 ZSWENO 算則結果誤差表 JSWENO ZSWENO Eh 0.123 0.141 Eq 0.019 0.024 max_eh 0.016 0.021 max_eq 0.018 0.014 註： Eh：水深總絕對值誤差；Eq：單寬流量總絕對值誤差 max_eh：水深最大絕對值誤差；max_eq：單寬流量最大絕對值誤差 Distance(m) 0 5 10 15 20 25 Water Le v el( m) 0.0 0.1 0.2 0.3 0.4 Bed level Exact JSWENO 圖 4- 5 JSWENO 模擬水位與解析解之比較圖
Distance(m) 0 5 10 15 20 25 q(m 2 /s) 0.175 0.180 0.185 0.190 0.195 0.200 Exact JSWENO 圖 4- 6 JSWENO 模擬 q 與解析解之比較圖 Distance(m) 0 5 10 15 20 25 Water Le v el( m) 0.0 0.1 0.2 0.3 0.4 Bed Level Exact ZSWENO 圖 4- 7 ZSWENO 模擬水位與解析解之比較圖
Distance(m) 0 5 10 15 20 25 q(m 2 /s) 0.175 0.180 0.185 0.190 0.195 0.200 Exact ZSWENO 圖 4- 8 ZSWENO 模擬 q 與解析解之比較圖 經由以上兩個案例(a)與(b)之結果分析比較，推測 ZSWENO 算則 中在(a)這種較不受邊界條件影響的案例中表現較為良好，而對於(b) 案例似乎因為上游的流量條件給定，為保持不連續面單寬流量的守恆 性，而過度收斂導致不連續面前後的點位連接有較差的模擬結果。 因此吾人提出一個理論想法，對於原始平滑指示器中的二次微分 項 ，假設其對於平滑指示器的影響程度，與組合平滑指示器之五 點組成的幾何關係有其相關性，存在著一組 使在構築平滑指示器 時能在不連續面能保持其強化收斂的特性，又能盡量保持不連續面前 後點位數值解的精確度。以此概念將(2-40)式平滑指示器改寫為： (4-1)
，定義為各計算元 項的修正係數。 另假設 值與五點之幾何組合的相關性，能由一未知的函數式 來求解出，本文應用類神經網路類神經網路的黑盒分析概念，以(b) 案例為學習案例來建構類神經網路，以擬合此假設的函數式。
4-2 模式架構
本節就模式的架構，分其類神經網路與水理部份分別介紹。4-2-1 水理模式
在水理計算方面，綜合前面第四章的理論介紹，本文水理模式可 整理出以下結論： 統御方程式：一維淺水波方程式。 空間離散式：守恆形式差分式 WENO 算則結合有限體積法。 時間離散式：三階 TVD- Runge-Kutta 積分。4-2-2 類神經網路模式
以 4-1 節的理論為基礎，目標以類神經網路建構出一套組合平滑 指示器五點組成的幾何關係，與 值相關的函數式，以下分輸入與 輸出參數分別討論，網路概念圖如圖 4- 9 所示：(a) 輸入參數： 在計算數值通量為正時依(2-44)式，依序採用(4-2)式五點之值 作為輸入；數值通量為負時依(2-47)式，則依序採用(4-3)式五點 作為輸入： 、 、 、 、 (4-2) 、 、 、 、 (4-3) 由於需保證類神經網路的輸入值能位於一固定值域，故在進 入類神經網路前，定義 ，為輸入五點絕對 值之最大值，將各輸入值除以 正規化，以使輸入值值域 位於-1 ~ 1 之間，重新定義類神經網路輸入參數為 、 、 、 、 (4-4) 、 、 、 、 (4-5) 。當 時， 為 除以 0 的型態， 此為不定式而可做任意解解釋，但無論代入任何值皆不影響最 後平滑指示器皆等於 0，因此定義此情況下 皆為 1 值， 以防止程式在分母等於 0 無法運算的狀況。 另外由於空間離散式是採用守恆形式差分式 WENO 算則，由
(3-30)與(3-31)式可知無論在對流項或是源項的計算上，(4-4)與 (4-5)式之 。 (b) 輸出參數： 輸出值代表修正係數 ，故輸出層採計三個單元。由於 本文自訂之類神經網路，處理單元的轉換函數皆採用雙彎 曲函數，轉換後之輸出值 ，符合 的值域，故在輸 入值經類神經網路傳遞後，所得之三個輸出值即依序為修 正係數 、 、 ，不需再另外的處理。
類
神
經
網
路
輸入1 輸入2 輸入3 輸入4 輸入5 輸出1 輸出2 輸出3 前處理 圖 4- 9 網路概念圖4-3 學習案例
4-3-1 案例選擇
學習案例為前面所介紹之(b)案例。此案例為穿臨界流況，在穩態 時有水躍的情形發生，水面線由平滑的漸變段與不連續面所組成。圖 4- 10 至圖 4- 13 為穩態時(3-30)與(3-31)式之各點位 組合 的分佈值，上標 p 代表運算矩陣中第 p 個位置。由前文的介紹中， 每組組合內五點之間的分佈差異，是為其決定 WENO 算則權重的重要因子。由以上圖 4- 10 至圖 4- 13 可明顯看出震波在 通過障礙物時，容易產生較為劇烈的變化與震盪，而如何求取一套計 算 WENO 算則的權重法則，使數值模擬在震波的劇烈的變化下，能 收斂使其不震盪，並能保持精確性，此乃自訂之類神經網路所希望能 學習到的結果。學習過程中所用到之相關參數請見 4-3-2 節說明。 distance(m) 9.5 10.0 10.5 11.0 11.5 12.0 value -0.01 0.00 0.01 0.02 0.03 i-2 i-1 i i+1 i+2 圖 4- 10 穩態時 分佈圖distance(m) 9.5 10.0 10.5 11.0 11.5 12.0 12.5 v alu e -0.015 -0.010 -0.005 0.000 0.005 0.010 i-2 i-1 i i+1 i+2 圖 4- 11 穩態時 分佈圖 distance(m) 9.5 10.0 10.5 11.0 11.5 12.0 12.5 valu e -0.005 0.000 0.005 0.010 0.015 0.020 i-1 i i+1 i+2 i+3 圖 4- 12 穩態時 分佈圖 distance(m) 9.5 10.0 10.5 11.0 11.5 12.0 12.5 value -0.5 0.0 0.5 1.0 i-1 i i+1 i+2 i+3
4-3-2 參數設定
模式參數設定分為水理、類神經網路與基因演算法三部份說明。 (a) 水理模式： 在 3-3 節中探討了乾濕點交界的處理方法，而其中有定義一 最小水深容許值 ，本學習案例雖無處理乾床的問題，但在第 六章詴驗案例的結果，會對其做更深入的探討。此外定義平滑 指示器中一微小數值 。 (b) 類神經網路： 網路架構如圖 4- 14 所示，依葉怡成(2001)書中建議使用單層 隱藏層，處理單元數輸入層為 5，隱藏層為 4，輸出層為 3。處 理單元之轉換函數則全採取雙彎曲函數 ，最大循環數 = 1000 次。 輸入層 隱藏層 輸出層 Wij Wjk θk θj 圖 4- 14 類神經網路架構圖(c) 基因演算法： 一開始需給定初始基因組群，本案例使用總基因組數=15。將 類神經網路連結加權值 W 與門閥值 θ，以二進位法轉化為一基因 序列的組合，可參考附錄 B 遺傳演算法介紹之編碼部份，定義： 1. 連結加權值 W： 基因長度 L = 8、上限 與下限 。 2. 門閥值 θ： 基因長度 L = 8、上限 與下限 。 適合度函數主要反應學習案例中，模擬結果與解析解的誤差大 小。由於模擬結果有兩組參數為水深 h 與單寬流量 q，兩者之數 量級大致差十倍左右，在此採取三種不同的代表誤差，如表 4- 2 所示，來測詴水深 h 與單寬流量 q 在學習案例中，哪一個參數影 響為大。 定義 (4-6) ：第 組基因組代表誤差； ：第 組基因組水深總絕對值誤差； ：第 組基因組單寬流量總絕對值誤差； 為一常數。 表 4- 2 方案 A ~ C 選用常數與代表誤差表
Type A Type B Type C
常數 N 1 10 100
適合度函數定義為：
(4-7) 上式中， ：第 i 組基因組適合度值。當 ， 最大值 ；當 ， ，而由於對數曲線在極大負值 處，數值下降斜率已趨近於 0，為避免所有 皆過大，造成 皆過小造成適合度鑑別率不高的情況，故在此定義一適合度修正係數， ， 為一數學運算函數，為傳 回最接近參數 x 的整數值，藉此一修正係數，可以使多數 能 位於對數曲線變化斜率較大的值域，以提高優選的作用。 此外若由類神經網路所得的權重，無法使數值穩定傳遞，可能會 造成數值震盪與誤差的放大，因而導致程式意外中止，為避免此情形 發生，對於學習案例-穩態流流經障礙物問題，設下一中止條件： 當計算至水深 h <0 時，立即中止模式的水理計算，定義此基因組 所對應的適應度 ，使用另一組基因組重新開始運算。 後續在基因組交配、突變的流程中，其中所需參數設定如下： 1. 每一學習循環交配數 Ma_N = 5。 2. 交配機率 Ma_Rate = 0.9。 3. 突變機率 Mu_Rate = 0.01。
4. 採用菁英化，將此次學習循環適合度最佳值的基因組直接複製， 以待下一次循環使用。
第五章 模式驗證與結果分析
5-1 學習結果
經由表 4- 2 方案 A ~ C，類神經網路各學習測詴二次所得之最佳 學習結果誤差如下表 5- 1，水位與單寬流量與解析解比較圖則如圖 5- 1 至圖 5- 12 所示。 表 5- 1 類神經網路學習結果誤差表JSWENO Type A Type B Type C Test_01 Test_02 Test_01 Test_02 Test_01 Test_02 Er 0.117 0.117 0.287 0.287 1.682 1.676 Eh 0.123 0.096 0.096 0.101 0.103 0.151 0.151 Eq 0.019 0.021 0.021 0.019 0.018 0.015 0.015 max_eh 0.016 0.013 0.013 0.014 0.015 0.022 0.022 max_eq 0.018 0.019 0.019 0.018 0.018 0.012 0.012 由誤差表可看出，三種方案的 Er、Eh、Eq 的值，兩次學習測詴大 致都會收斂至各自一定的範圍內。從誤差與解析解的比較圖，可知 A 方案雖然在水深 h 有較好的模擬結果，但是單寬流量的守恆上，較原 本 JSWENO 算則有較差的結果；C 方案所得之模擬結果較類似於 ZSWENO 算則，雖然降低了在不連續面的最大單寬流量誤差，但也 造成了在不連續面附近的水位與單寬流量模擬有較差的結果；B 方案 所得之學習結果，在水位上無論是總絕對值誤差與單一最大誤差都有 其較佳的結果，而在單寬流量的模擬上，雖然最大誤差有些許增加， 但由模擬單寬流量圖可發現在不連續面上，與 JSWENO 算則相比有
較好的流量守恆性。故以 B 方案作其學習的最佳方案，並選用其中 代表誤差較小之第一次測詴，學習所得之類神經網路做為最終學習成 果，圖 5- 13 為其學習過程中代表誤差的收斂圖。定義其代號為 AWENO(ANN-WENO)算則，以下以此算則針對不同詴驗案例做模擬 ，並與 JSWENO 與 ZSWENO 算則之結果做分析與探討。 d is t a n c e ( m ) 0 5 1 0 1 5 2 0 2 5 W a te r L e v e l 0 .0 0 .1 0 .2 0 .3 0 .4 b e d le v e l E x a c t T y p e - A ( 1 ) 圖 5- 1 方案 A (1)模擬水位與解析解之比較圖 d is t a n c e ( m ) 0 5 1 0 1 5 2 0 2 5 q (m 2/s ) 0 .1 7 5 0 .1 8 0 0 .1 8 5 0 .1 9 0 0 .1 9 5 0 .2 0 0 E x a c t T y p e - A ( 1 ) 圖 5- 2 方案 A (1) 模擬 q 與解析解之比較圖
d is t a n c e ( m ) 0 5 1 0 1 5 2 0 2 5 W a te r L e v e l 0 .0 0 .1 0 .2 0 .3 0 .4 b e d le v e l E x a c t T y p e - A ( 2 ) 圖 5- 3 方案 A (2)模擬水位與解析解之比較圖 d is t a n c e ( m ) 0 5 1 0 1 5 2 0 2 5 q (m 2/s ) 0 .1 7 5 0 .1 8 0 0 .1 8 5 0 .1 9 0 0 .1 9 5 0 .2 0 0 E x a c t T y p e - A ( 2 ) 圖 5- 4 方案 A (2) 模擬 q 與解析解之比較圖 d is t a n c e ( m ) 0 5 1 0 1 5 2 0 2 5 W a te r L e v e l 0 .0 0 .1 0 .2 0 .3 0 .4 b e d le v e l E x a c t T y p e - B ( 1 ) 圖 5- 5 方案 B (1)模擬水位與解析解之比較圖
d is t a n c e ( m ) 0 5 1 0 1 5 2 0 2 5 q (m 2/s ) 0 .1 7 5 0 .1 8 0 0 .1 8 5 0 .1 9 0 0 .1 9 5 0 .2 0 0 E x a c t T y p e - B ( 1 ) 圖 5- 6 方案 B (1) 模擬 q 與解析解之比較圖 d is t a n c e ( m ) 0 5 1 0 1 5 2 0 2 5 W a te r L e v e l 0 .0 0 .1 0 .2 0 .3 0 .4 b e d le v e l E x a c t T y p e - B ( 2 ) 圖 5- 7 方案 B (2)模擬水位與解析解之比較圖 d is t a n c e ( m ) 0 5 1 0 1 5 2 0 2 5 q (m 2/s ) 0 .1 7 5 0 .1 8 0 0 .1 8 5 0 .1 9 0 0 .1 9 5 0 .2 0 0 E x a c t T y p e - B ( 2 ) 圖 5- 8 方案 B (2) 模擬 q 與解析解之比較圖
