行政院國家科學委員會專題研究計畫 期中進度報告
子計畫:參數即時校正河川洪水演算模式之研發與應用
(2/3)
計畫類別： 整合型計畫 計畫編號： NSC93-2625-Z-002-003- 執行期間： 93 年 08 月 01 日至 94 年 07 月 31 日 執行單位： 國立臺灣大學生物環境系統工程學系暨研究所 計畫主持人： 許銘熙 計畫參與人員： 吳建均, 傅金城, 林洙宏 報告類型： 精簡報告 處理方式： 本計畫可公開查詢中 華 民 國 94 年 5 月 30 日
行政院國家科學委員會專題研究計畫期中進度報告
參數即時校正河川洪水演算模式之研發與應用(2/3)
Development and Application of Flood Routing Model with Real-Time
Parameter Correction (2/3)
計畫編號：NSC 93 - 2625 - Z - 002 - 003
執行期限：93 年 8 月 1 日至 94 年 7 月 31 日
計 畫 主 持 人 ：許 銘 熙 Ming-Hsi Hsu
研 究 助 理：吳 建 均 Jian-Cyun Wu
研 究 助 理：傅 金 城 Jin-Cheng Fu
研 究 助 理：林 洙 宏 Shu-Horng Lin
執 行 單 位：國立台灣大學生物環境系統工程學系
一、中文摘要 本研究計畫屬於「即時洪水預報模式 之研發與應用」整合型計畫之子計畫三─ 「參數即時校正河川洪水演算模式之研發 與應用」，全程計畫為期三年，目前為第 二年。 蘭陽溪所處之蘭陽平原因地勢低窪， 夏秋之際常因梅雨季或颱風所帶來的豪雨 造成淹水災害，因此，若能提供沿岸的洪 水位剖面線及特定地點的洪水位預報資 訊，作為淹水預警之參考，將可有效減低 洪水所帶來的災害損失。故本計畫之目的 即為針對蘭陽溪的特殊地理環境及水文特 性，建立合適的參數即時校正河川洪水演 算模式。 第一年度在現地踏勘及資料蒐集後， 根據動力波理論及數值方法建立河川洪水 演算模式及潮位預報模式。 第二年度接續第一年度的研究成果， 進一步的藉由觀測水位的獲得，建立參數 即時校正河川洪水演算模式，並與各子計 畫進行模式演算成果之串聯銜接，執行整 體洪水預報模式測試工作。 關鍵詞：洪水預報、參數校正。 AbstractThis project is the 3rd-subproject of the integrated project “Development of a
Real-Time Forecasting Model and Its Application”.
The Lan-Yang Creek basin is located in the low-lying Lan-Yang Plain, where the inundation disasters occur frequently in summers and falls because of the torrential rain. If the real-time forecast information, including the water stage at the significant locations and longitudinal profiles along the river is available before flooding, the damages would be effectively mitigated. The purpose of this study is to develop a flood routing model with real-time parameter correction for the Lan-Yang Creek.
In the first year, a numerical flood routing model and tidal stage forecasting technique have been developed for the Lan-Yang Creek.
According to the first year of development plans, the project establish the flood forecasting with real-time parameter correction taking the observed river stages in the second year. Besides, we cooperate with the other subprojects to examine and compare the real-time forecasting model.
Keywords: Flood forecasting, Parameter
蘭陽溪為台灣東北部宜蘭縣境內主要 河川，發源於南湖大山北麓(標高 3,536 公 尺)，本流向東北流至破布烏始入平地，至 蘭陽橋附近形成幹流，再延東北流至東港 注入太平洋，流域面積約 978 平方公里， 主幹流長度約 73 公里。流域下游為蘭陽平 原，地勢平坦，土壤肥沃，向以農業為經 濟基礎，近年來因受工商業日漸發展之影 響，使得週邊交通網路逐漸發達，勢必大 幅提昇宜蘭地區之繁榮。 但由於蘭陽溪流域年平均雨量高達 2,656 公釐，並因流域土質屬排水性較差之 砂土性質，導致洪水宣洩不良而氾濫成 災，尤其是每年夏、秋的颱風期間，受到 暴雨來襲而造成災害。近年來相關單位雖 以工程治理作為防洪手段，但受到工程措 施仍有其設計風險存在，並無法完全消除 水患，故應輔以非工程之防洪措施，如洪 水預警及洪災保險措施之規劃等，以彌補 工程方法之不足。 本子計畫針對蘭陽溪河系地文與水文 特性，建立河川洪水預報模式。主要目的 係承接上游逕流模式(子計畫一)之預報流 量，進行下游河系之洪水預報模擬，以提 供未來 1~6 小時之洪水資訊。本子計畫以 一維緩變量流動力波模式作為河川洪水預 報模式之主體。並在納入研究區域內各水 文控制點之即時觀測水位後，立即修正河 川洪水模式之輸入資料，並期望經由即時 校正回饋演算後所得之預報水位能更接近 觀測水位，以即時提供蘭陽溪河系沿岸的 洪水位剖面線及若干特定地點的洪水位資 訊，作為蘭陽溪、羅東溪及宜蘭河等兩岸 洪水平原及區域淹水預警之參考，並支援 決策系統，以達到流域防洪減災之目的。 研究過程中，本計畫需銜接子計畫一 的逕流量作為模式上游邊界，並將子計畫 二之研究成果及觀測水位作為模式內邊界 條件予以調整流況，最終研究成果則配合 子計畫四將其演算成果予以展示，作為預 報研究區域內重要地點的洪水位及河川沿 岸的水面剖線之參考，並加強各子計畫間 資料之整合，以達成總計畫所規劃之研究 目標。 三、文獻回顧 有關參數即時校正河川洪水演算模式 的相關研究，本計畫分成為河川變量流模 式之發展與應用、模式流況修正及模式參 數校正之研究等三方面，茲分別說明如下。 3.1 河川變量流模式發展及應用之研究 河川中水流之流動係隨時間及位置而 改變，其流況通常可依一維連續及運動方 程式所組成之緩變量流方程式來描述，但 由於此變量流方程式相當複雜，一般情況 下無法得到解析解，必須利用數值方法求 解以應用至實際狀況。近年來隨著計算機 能力的增強及相關理論的進展，河川洪水 演算模式已成為洪水預警有力的工具之 一。此處將簡介國內、外學者所建立的河 川洪水演算模式。 在國外方面，Mahmood 與 Yevievich (1975)和 Abbott(1976)分別以較完整的理 論分析渠道變量流及其水理。Englund 和 Fredsoe(1976) 發 展 可 作 為 模 擬 一 維 之 河 川、渠道、灌溉系統、降雨逕流及流移-延散現象的一維河川及渠道之模式系統 (MIKE 11)。Chang(1976)發展適用於沖積 河川及間歇性河川的一維變量流動床模式 (FLUVIAL-12)。Yen(1985)以迪聖凡納方程 式(de St. Venant equation)為主體，分別將 動力波、準變量波、擴散波及運動波之差 分方程式作一概述，並且探討其適用性。 Schaffranek(1987a,b)採用四點權重、隱性 之有限差分近似法發展分支型式的網路動 力流動模式(branch network dynamic flow model, BRANCH)。Einhelling 等(1990)發 展可用於辮狀沖積河槽水面及河床演變之 長 期 模 擬 的 一 維 變 量 水 流 及 沉 滓 模 式
(CHARIMA) 。 Parker(1991) 多 次 修 改 FLUVIAL-12 模式，使其功能更為完善。 Wiley(1993)發展具有資料校正、可預測流 量及流速功能的網路型態渠道之河川流量 演算模式(DAFLOW)。Nguyen 和 Kawano (1995)以雙掃法模擬求解無迴水影響渠流 的洪水演算。Bajracharya 和 Barry (1997) 以截尾誤差分析(truncation error analysis) 求解運動波並發展為洪水演算模式。 國內方面，許時雄(1969)建立了淡水 河下游接近河口感潮段之演算模式。周乃 昉(1979)以變量流理論建立非線性四點完 全隱式法來模擬淡水河洪水時的水理。顏 清連(1979~1983)使用非線性四點完全隱 式差分法配合牛頓-拉福生法求解完整的 迪聖凡納方程式發展變量流模式，並針對 台灣河川之地文與水文特性，以理論分析 及數值試驗等方式探討河川坡度、自然河 川斷面不規則變化、曼寧糙度、匯流洪峰 比及上、下游邊界條件對數值模式的穩定 性影響。在此之後(1983~1987)並以該變量 流模式從事淡水河洪水模擬及預報為主的 一系列研究。蔡長泰(1977~1991)以線性完 全隱式差分法配合雙掃法求解完整的迪聖 凡納方程式發展洪流演算模式，該模式首 先以朴子溪做為模式的驗證與探討，之後 更將該模式應用於濁水溪及曾文溪洪水預 報 。 賴 經 都 (1995) 利 用 多 方 式 特 徵 法 (MMOC)以耦合求解方式求解單一河道變 量流水理的問題。李鴻源等(1996)發展出 一個適用於台灣獨特蜿蜒河川的河川網路 之擬似二維輸砂模式(NETSTARS)，能夠 更適切且合理地模擬河川中主支流交匯狀 況下之複雜河系的運動行為。許銘熙(1984, 1996~1998)延續顏清連所建立之模式，並 在考量地形變化過劇及提昇演算精度的情 況下，增加了運動波的演算及即時水位校 正的功能，該模式目前已納入淡水河洪水 預警系統。 3.2 模式流況修正之研究 為使計算水位能在質量及動量守恆下 趨近於觀測水位，則必須要廣泛的蒐集相 關於如何將觀測資料輸入數值模式的研 究，以下則概略的介紹國內、外學者在此 方面的研究成果。 Kuchment 和 Nazarov(1999)以對流擴 散方程式建立河川洪水預報模式，當模式 獲得觀測流量時，立即以線性內插的方式 配合邊界條件修改整個河道的流況，以作 為 演 算 下 一 時 刻 河 川 流 況 的 初 始 值 。 Camacho 和 Lee(1999)採用多線性離散數 值方法作為河道演算的基礎，並配合河川 觀 測 流 量 的 獲 得 而 發 展 為 預 報 模 式 。 Kothyari 和 Singh (1999)以時雨量和日流 量配合統計的多變數迴歸法建立了預測日 流量的關係式，並提出納入即時觀測資料 對於預報的準確度有莫大的影響。Hsu 等 (2000)先以一維緩變量流理論建立了淡水 河洪水預報模式，然後對於每小時回傳的 觀測水位除了取代數值模式的計算值外， 並將此觀測水位的條件取代演算過程中的 動量方程式，雖然此方式可以在保持質量 守恆的情況下使計算水位等於觀測水位， 但將會造成水面剖線的不連續。Carmeron 等(2002)配合觀測值的輸入以類神經網路 建立預報模式，並在經過驗證後發現對於 較小的洪水事件，觀測值可以使預報流量 有較為合理的預測結果。Hsu 等(2300)以 觀測水位配合最小平方法即時校正變量流 模式流況，且其應用於短延時(1~2 小時) 預報具有良好的效果。 3.3 模式參數校正之研究 系統參數對於數值模式的演算結果具 有重大的影響力，因此大多數的數值模式 在應用前都必需對數值模式的參數進行檢 定，但受限於該過程甚為耗力費時且可能 受到主觀因素的影響，故如何作有效率的
參數檢定以符合洪水預報的時效性，則必 須廣泛的探討國內外學者所使用的參數檢 定方式及相關的即時修正理念，以下則概 略的介紹在此方面的研究成果。 Chow (1973)藉由實驗與現地調查資 料等建立了曼寧 n 值選用的標準，並表列 供研究設計參考。Chiu 和 Emmanuel(1978) 以 卡 門 濾 波 配 合 觀 測 水 深 推 估 曼 寧 n 值。Wormleaton 和 Karmegam (1983)以地 文因子(渠底寬、邊坡坡度)及水力因子(水 力半逕與通水面積)，在目標函數為流量及 水深誤差最小的條件下以最佳化的最小平 方法，求解最適當之曼寧 n 值以縮減計算 水 位 及 流 量 與 觀 測 值 間 的 誤 差 。 Jarrett (1984)在蒐集急流河川的現地資料後，以 多變數迴歸法建立摩擦坡度、水力半徑和 曼寧 n 值的關係式，故只需量測河川摩擦 坡度和水力半徑即可預測曼寧 n 值。Fread (1992) 藉由試誤檢定方法配合洪水事件 的觀測值，推算出最適合該河川的曼寧 n 值。Alain (1992) 使用曼寧公式加上疊代 法來估算以迪聖凡納方程式所建立的洪水 演算模式中的曼寧 n 值，並在考慮到洪水 平原的影響下應用至實際河川。Crissman 和 Chiu(1993)為探討 Niagara River 在冬天 時冰對於河川阻力的影響，以迪聖凡納方 程式配合卡門濾波建立河床糙度曼寧 n 值可隨時間變化的預測模式。Sepaskhah 和 Bondar(2002)依據灌溉、入流量和植生 覆蓋率而推求出兩個可以估算曼寧 n 值 的方程式。Wong 和 Zhou(2003)先建立了 水深及曼寧 n 值之間的關係，並在執行運 動波模式的演算過程中，藉由水深的變化 而即時修正曼寧 n 值。Liu 等(2003)在執 行水理演算的過程中，藉由河川水深及植 生作物的浸沒高度而計算出河道的曼寧 n 值以即時反應河道阻力。 綜合整理以上所述前人之研究，可以 獲得如下之結論： 1.為適應蘭陽溪流域的特殊地理環境及配 合數值模式原始程式的取得，故本計畫 以顏與許(1996-1998)經由多方面測試所 發展的變量流模式為基礎，進一步發展 參數即時校正河川洪水演算模式。 2.在觀測資料輸入數值模式的研究方面， 本計畫在考量水面剖線連續性、質量守 恆及時效性的狀況下，採用最小平方法 配合觀測水位進行預報初始值的設定。 3.在參數校正方面，本計畫先使用試誤法 進行模式參數的檢定，作為洪水預報期 間，參數校正的初始值。然後，再使用 最佳化理論建立觀測水深和河段曼寧係 數間的關係，作為模式的參數校正方式。 四、模式理論 本計畫於數值模式方面主要係針對顏 與許(1996-1998)等人所建立完成的河川洪 水預報模式為基礎，再發展成參數即時校 正河川洪水演算模式。並藉由即時觀測水 深，以修正預報初始值及系統參數-曼寧 n 值，而獲得合理的即時流況及精確的洪水 預報資訊。有關參數即時校正河川洪水演 算模式之控制方程式、邊界條件、數值方 法、修正預報初始值及參數校正之方式， 茲說明如后。 4.1 河川洪水演算模式 本 子 計 畫 根 據 迪 聖 凡 納 氏(de Saint Venant)所導出之一維緩變量流之動力波 理論來描述河川中水流動態。本子計畫為 考慮沿岸殘流域流量及感潮河段之影響， 將研究區域範圍之上游邊界延伸至蘭陽溪 之家源橋、羅東溪之鼻頭橋與宜蘭河之員 山橋。然而上游河道延伸至新上游邊界 後，進行動力波數值演算過程會發生數值 不穩定現象，經探究其原因為河道坡度過 於陡峭所致，在蘭陽溪上游葫蘆堵大橋至 家源橋之河段 (斷面編號 27~74)在 27,878
公 尺 長 度 內 高 程 由 10.78 公 尺 陡 升 至 367.21 公尺；另外在羅東溪上游河段，斷 面編號 009~014(鼻頭橋)之河段在 3,348 公 尺長度內高程由 9.87 公尺上升至 51.605 公尺。本子計畫為使河川演算模式能達到 數值穩定的範圍內，於防汛期間穩定提供 正常預報功能，故於蘭陽溪河段家源橋至 牛鬥橋；羅東溪河段鼻頭橋至 009 斷面等 兩處河段，以較不受地形影響之運動波模 式演算，再銜接動力波模式，完成蘭陽溪 河系之洪水演算。 4.1.1 動力波模式 1.理論背景 將壓力、重力、摩擦力、局部加速度 及對流加速度等項完整考慮之水流運動方 程式，則稱為動力波。其水流連續及運動 方程式依次為 0 2 1+ = − + q q x Q t A ∂ ∂ ∂ ∂ (1) 0 2 1 1 2 = + − − − − + A Q q V q S x Y S gA A Q x t Q f o ∂ ∂ ∂ ∂ ∂ ∂ (2) 其中，A 為通水斷面積；Q 為河川流量；q1 為單位河川長度之側入流量；q2為單位河 川長度之溢出流量；t 為時間座標；x 為沿 水流方向之距離座標；g 為重力加速度；S0 為河道底床坡度；Y 為河道斷面之水深； f S 為摩擦坡度；V1為河道側流平均流速在 主流方向之分量。在數值演算的過程中， 可以藉由曼寧公式(Manning’s formula)將 摩擦坡度表示為 2 3 / 4 2 A R Q Q n Sf = (3) 式中，R 為水力半徑；n 為曼寧係數，代 表河川之糙度值。曼寧 n值通常經由颱洪 事件的檢定與驗證後決定，在一般數值模 式演算的過程中通常假設曼寧n值不會隨 著時間而變化。 (1)與(2)式分別表示水流之連續及運 動方程式，若忽略河川斷面之底床沖淤， 則通水面積(A)與水深(Y)間具有一定函數 關係，故(1)及(2)式包括水深(Y)及流量(Q) 兩組未知的流況應變數。然而(1)與(2)式均 為雙曲線型偏微分方程式，在天然河川中 其流況應變數(Q、Y)通常無法直接求得解 析解，而必須利用數值方法進行求解。 2.數值方法 本計畫利用非線性完全隱式法(fully implicit scheme)求解各時段在河川各斷面 的水深及流量。若以 f 表示欲求解之流況 應變數(包括 Q 或 A)，依四點完全隱式差 分法寫成的差分式可被表示為(Amein 和 Fang, 1970)
(
)
(
)
(
) (
)
(
t)
j t j t j t j t j t j t j t j f f f f t t f f f x x f f f f − + − ∆ = ∂ ∂ − ∆ = ∂ ∂ + = + + + + + + + + + + 1 1 1 1 1 1 1 1 1 1 2 1 1 2 1 (4) 式中，t、j分別代表時間與空間座標指數； j j x x x= − ∆ +1 為河川相鄰斷面之距離；∆t為 演算時距。 將(4)式代入(1)式及(2)式中，即可將任 意的兩相鄰斷面間(xj,xj+1)之連續方程式 及運動方程式，表示為前進時間段的流量 及水深之函數關係，分別如Cj及Mj函數 所示(
)
(
, , , , , , ,)
0 0 , , , , , , , 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 = = + + + + + + + + + + + + + + + + t j t j t j t j t j t j t j t j j t j t j t j t j t j t j t j t j j Y Q Y Q Y Q Y Q M Y Q Y Q Y Q Y Q C (5) 式中， j=1,2,3,…., J為河道斷面指標，共 計有 J 個斷面；Cj及Mj分別為藉由第 j 個和第(j+1)個斷面所組成的連續方程式及 運動方程式；t 為現在時刻；而 (t+1) 則 為下一時刻((5)式之詳細型式詳見參考文 獻(許銘熙, 1984))。 每兩相鄰斷面可寫出一組如(5)式之 聯立方程組，若河道有 J個斷面，則欲求 解之流況應變數共有 2J 個，經由(5)式可 以得到 2(J-1)組之差分方程式，再配合河 道上、下游端之邊界條件，可建立成 2J 個聯立方程式，便可推求出 2J 個流況應 變數。將(5)式配合上、下游邊界條件，依河道斷面次序排列展開可以得到
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
,)
0 0 , , , , , , , 0 , , , , , , , 0 , , , , , , , 0 , , , , , , , 0 , , , , , , , 0 , , , , , , , 0 , 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 1 1 1 1 1 2 1 2 1 1 1 2 2 1 1 1 1 1 2 1 2 1 1 1 1 1 0 = = = = = = = = + + − − + − + − + + − − − + − + − + + − + + + + + + + + + + + + + + + + + + + + + + + + + + t J t J J t J t J t J t J t J t J t J t J J t J t J t J t J t J t J t J t J J t j t j t j t j t j t j t j t j j t j t j t j t j t j t j t j t j j t t t t t t t t t t t t t t t t t t Y Q M Y Q Y Q Y Q Y Q M Y Q Y Q Y Q Y Q C Y Q Y Q Y Q Y Q M Y Q Y Q Y Q Y Q C Y Q Y Q Y Q Y Q M Y Q Y Q Y Q Y Q C Y Q C Μ Μ Μ Μ (6) (6)式中共有 2J 個聯立方程式，可解 得 2J 個流況應變數。但因為(6)式係為非 線性聯立方程式，無法直接求解，故須利 用疊代法求解，本計畫採用牛頓-拉福生法 (Newton-Raphson method)逐次疊代修正以 逼近真解，即(6)式之解可以表示成( ) ( ) ( )
( ) ( ) ( )
m j m t j m t j m j m t j m t j Y Y Y Q Q Q δ + = δ + = + + + + + + 1 1 1 1 1 1 (7) 式中，m = 0,1,2,3,…為疊代修正之次數； j=1,2,3,…., J 為斷面指標，共計有 J 個斷 面； t+1 j Q 和Yjt+1分別為(t+1)時刻的未知流 況應變數；δQj及δYj表示每次疊代修正 量。若假設任意一組之疊代修正起始值 0 1 ] [Qij+ 及 0 1 ] [Yji+ ，將其代入(6)式中，可以 發現左邊之函數值並不為零，而得到一組 殘餘值rj，此結果即表示假設之疊代修正 起 始 值 1 0 ] [Qij+ 及 1 0 ] [Yji+ 並非為(6)式之真 解，而需再重新的進行疊代修正求解。當 m = 0 時, 則採用 t 時刻的已知流況變數 作為起始值，即( )
t j t j Q Q + = 1 1 和( )
t j t j Y Y + = 1 1 。 根據牛頓-拉福生法，每次的疊代修正 量δQj及δYj可由下式的推求獲得 − − − − − − − − = ⋅ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ − − − − − − − − − − − − − − − − + + + + J J J j j J J J j j J J J J J J J J J J J J J J J J J J J J j j j j j j j j j j j j j j j j r r r r r r r r Y Q Y Y Q Q Y Q Y M Q M Y M Q M Y M Q M Y C Q C Y C Q C Y M Q M Y M Q M Y C Q C Y C Q C Y M Q M Y M Q M Y C Q C Y C Q C Y C Q C 2 1 2 2 2 2 1 2 3 2 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 2 1 1 1 1 1 2 1 2 1 1 1 1 1 1 0 1 0 Μ Μ Μ Μ Ο Ο Ο Ο Ο Ο δ δ δ δ δ δ δ δ (8) 或將(8)式表示成 i j j i b a, ⋅δ = (9) 式中，i = 1,2,3,…,J為方程式指標，方程 式總數共 2J， j = 1,2,…,J為流況應變數指 標，流況應變數總數共 2J，ai,j 為係數方 矩陣(square matrix)，δj為待求解之疊代修 正量的行矩陣(column matrix)，b
i為殘餘 值的行矩陣。 (9)式係為線性方程式組，可以使用樞 紐列高斯消去法(Gauss elimination method with pivot row)求解河道各斷面之流況疊 代修正量(δQj及δYj)，以增加計算之速 度。由(8)式計算所得各斷面之流況疊代修 正量(δQj及δYj)代入(7)式後，可得到疊代 修 正 後 之 各 斷 面 流 況 應 變 數 ( t+1 j Q 及 1 + t j Y )，此時檢驗求得之相對疊代誤差量 1 + t j j Q Q δ 和 t+1 j j Y Y δ ，若小於一容許值 (本子計畫採用 0.01)，則此時殘餘值(rj) 將接近於零，代表各斷面之流況應變數 (Qtj+1及 1 + t j Y )已接近(6)式之真解。否則必需 將各斷面之流況應變數( t+1 j Q 及 t+1 j Y )重新 代入(8)式中，進行下一次的疊代修正，直 到滿足所需精度為止。依此步驟，可逐步 解得每一時刻河道中各斷面之流量與水 深。 若應用至河系時，則假設在匯流處的 水位相等及流量進出平衡的情況下，即能 建立匯流點下游河道的上游邊界條件，故 由此觀念，可使計算工作從上游河道通過 匯流點而演算至下游河道，建立一無參數 校正之河川洪水演算模式[許, 1984]。4.1.2 運動波模式 1.理論背景 水流運動方程式，若忽略局部加速 度、位移加速度及壓力梯度之影響時，即 僅考慮重力與摩擦力平衡時，其與連續方 程式組合之運動現象稱之為運動波，即： 1 λ q t A x Q+ = ∂ ∂ ∂ ∂ (10) S0 =Sf (11) 式中，A 為河川通水斷面積；Q 為河川流 量；
q
λ1為單位河川長度之側入流量；t 為 時間座標；X 為沿水流方向之距離座標； So為河床坡度；Sf為摩擦坡度。 式中 Sf若以曼寧公式表示，則(11)式 可以改寫為： Q nAR SO =1 2 3 1 2 (12) 或A=αQβ (13) 其中， , 0.6 3 2 = = β α o S nP 式中 P 為濕周，n 為曼寧糙度。將(13)式 代入(10)式中可得： ∂ ∂ αβ β ∂∂ Q x Q Q t q + * −1 = 1 λ (14) 此即為運動波模式所用以計算之方程 式。若α 、β已知，則(14)式中僅有一未 知數Q，可利用數值方法求解。 2.數值方法 (1)線性法： 線性法所採用之差分式為： x Q Q x Q ij i j ∆ − = ++11 +1 ∂ ∂ (15) ∂ ∂ Q t Q Q t j i j i = ++1− + 1 1 ∆ (16) 式中，i,j 分別代表空間及時間之座標指 數；∆x為河道相鄰斷面間之距離；∆t為 演算時距。為求得線性解，再將Q及q表 示成： Q Qj Q i j i ≈ +1+ + 1 2 (17) q q j q i j i λ1 λ λ 1 1 1 1 1 2 ≈ ++ + + (18) 將式(15)、(16)、(17)、(18)代入式(14) 中，可得線性運動波之差分式： Q Q X Q Q Q Q t q q j i j i j i j i j i j i j i j i + + + + + − + + + + + + − + + − = + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 ∆ αβ ∆ β λ λ (19) 由式(19)，可解得未知項 1 1 + + i j Q ： Q Q Q X Q Q Q t q q t x Q Q j i j i j i j i j i j i j i j i j i j i ++ ++ + + + + − + + + + + − = − + + + + + + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 ∆ ∆ ∆ ∆ αβ αβ β β λ λ (20) (2)非線性法： 將式(4-15)表成有限差分之型式如下： Q Q X A A t q q j i j i j i j i j i j i + + + + + + + + + − + − = + 1 1 1 1 1 1 1 1 1 1 1 2 ∆ ∆ λ λ (21) 又由式(13)式可得：( )
Ai Q j j i ++1 = ++ 1 1 1 α β (22)( )
Aij Q j i +1=α +1 β (23) 將式(22)、(23)代入式(24)中，可得非 線性運動波之差分式：( )
( )
∆ ∆ ∆ ∆ ∆ t XQ Q t XQ Q t q q j i j i j i j i j i j i + + + + + + + + + + = + + − 1 1 1 1 1 1 1 1 1 1 1 2 α β α β λ λ (24) 在式(24)中，等式右邊均為已知項， 若令其為C~，即得：( )
~ C t X Q Q t q q j i j i j i j i = + + − + + + + + ∆ ∆ ∆ 1 1 1 1 1 1 1 2 α β λ λ (25) 使 用 疊 代 法 計 算 時 ， 其 殘 數 誤 差( )
1 1 + + i j Q f (residual error)可表為：f Q t X Q Q C j i j i j i ( ++1)= ++ + ( ++ ) − ~ 1 1 1 1 1 ∆ ∆ α β (26) 應用牛頓疊代法：
( ) ( )
( )
( )
1 1 1 1 1 1 1 1 1 + + + + + + + + +=
−
′
i j i j k i j k i jQ
f
Q
f
Q
Q
(27)( )
( )
f Q t X Q j i j i ' ++11 = ∆ ++11 −1 ∆ αβ β (28) 式中，k代表疊代次數。當( )
1 1 1 + + + k i j Q f 小於允 許誤差時，可求得未知項 11 + + i j Q 。在應用非 線性法求解時，未知項( 11 + + i j Q )初值的假設對 疊代的收斂與否影響很大，故本模式在 ( 1 1 + + i j Q )的初值假設時先採用線性法所得之 結果，再經由非線性法修正，以求得最近 似的數值解。 4.2 參數校正方式 參數校正是期望經由河川洪水演算模 式中的系統參數之校正，以縮小預報水深 與觀測水深之間的水深差。在河川洪水預 報模式中與洪水位有直接關係的系統參數 為曼寧n值。故本子計畫期望藉由調整各 河段曼寧n值，能夠予以降低預報水深與 觀測水深間的水深差平方和。因此在水深 (Y)與曼寧係數為單值關係的假設下，目標 函數可定義為(
)
∑
= Ω Ω = − K k k k Y Y f 1 2 2 1 ) ( min (29) 式 中 ， f (Ω) 為 目 標 函 數 而 變 數{
n1,n2,n3,..., nK}
∈ Ω ；K為流域內可提 供觀測水深之水文站總數；k 為水文站指 標；Yk為計算水深；Yk則為觀測水深。 因(29)式為一非線性方程式，故本計 畫 利 用 最 佳 化 理 論 之 Gauss-Newton method (茨 木 俊 秀 等, 1997; Jorge 和 Stephen, 1999)。將校正前、後的曼寧係數 表示為[ ]
[ ]
[
] [
m]
T m m T m m m [ ] [ ] [ ] [ ] 1 1= Ω − Φ ⋅ Φ ⋅ Φ ⋅ Ε Ω + − (30) (30)式中 m K K m Y Y Y Y Y Y Y Y − − − − = Ε Μ 3 3 2 2 1 1 ] [ (31) m K K K K K K K K j i m n Y n Y n Y n Y n Y n Y n Y n Y n Y n Y n Y n Y n Y n Y n Y n Y n Y ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ = ∂ ∂ = Φ Λ Μ Ο Μ Μ Μ Λ Λ Λ 3 2 1 3 3 3 2 3 1 3 2 3 2 2 2 1 2 1 3 1 2 1 1 1 , ] [ (32) 式中，[Φ]m矩陣中之 i 和 j 分別表示水 文站和校正曼寧係數河段的編號。並定義 水文站編號 k=1,2,3,….,K，是從上游往下 游依序排列。 故(13)式之敏感係數，經由配合水文 站編號和可校正曼寧係數河段的編號，可 分成為三類 1.當i = j時 表示水文站與配合曼寧係數校正的河 段在同一河段位置，則本子計畫在假設演 算河道為寬淺斷面且為定量流的情況下， 則計算流量、計算水面寬度及摩擦坡度之 變化可被忽略，且水力深度相近於計算水 深。可將曼寧公式改寫為 5 / 3 2 / 1 = f BS Qn Y (33) 經由(33)式對曼寧係數偏微後，可以 將敏感係數表示成 j i j i f j i n Y BS Qn n n Y 5 3 , 5 / 3 2 / 1 , = ∂ ∂ = ∂ ∂ (34) 2.當i < j時 表示水文站位於配合曼寧係數校正的 河段上游。本計畫考量校正下游河段的曼 寧係數時所形成的水深變化，將會以迴水 的方式影響上游水文站的水深。故可將水文站受到迴水影響後的水深變化表示成
∏
− = ⋅ = ∂ ∂ 1 , 5 3 j i j j j i D n Y n Y λ λ (35) 式中， 1 2 1 1 − + + − = =∏
i i i j j i D D D D D Λ λ λ , Dλ，詳細推 導參見附錄一。 3.當i > j時 表示水文站在配合曼寧係數修正的河 段下游。在亞臨界流及相同邊界情況下， 上游河段的曼寧 n 值變化對於下游水文 站的水深影響甚微，則該敏感係數可被設 定為 0 , = ∂ ∂ j i n Y (36) 本計畫整理(34)、(35)及(36)式後，可 將敏感係數彙整為 = ∂ ∂∏
− = 0 3 5 3 5 1 , j i j j j i j i D n Y n Y n Y λ λ for i = j for i < j for i > j (37) 將(37)式之敏感係數代入(32)式，即可 推求出其敏感係數矩陣。再藉由(30)式， 即可疊代修正得到另一組更合適之河段曼 寧n值，再代入數值模式演算，而得到修 正後的計算水深。因此，藉由 t 時刻的觀 測值同時的修正預報初始值及曼寧係數， 而建立參數即時校正河川演算模式，其演 算流程如圖1所示。 4.3 邊界條件 河川演算模式之邊界條件的型式約可 分為三種；水位歷線、流量歷線和率定曲 線，大多數在從事河川變量流數值演算 時，上游端之邊界條件通常採用流量歷線 或水位歷線，下游之邊界條件則可以採用 流量歷線、水位歷線、或流量與水位關係 之率定曲線。 本計畫在河道上游端的邊界條件採用 子計畫一的預報流量，下游端的邊界條件 則使用本計畫所研發的潮位預報模式。 本計畫所研發的潮位預報模式，主要 以調和分析所得之天文潮位為基礎，在納 入颱洪期間氣象潮之影響，進行河口潮位 預報。茲列式如下： Η( )t
=η
P ( )t
+η
f ( )t
(38) 式中，Η(t)為河口預報潮位；ηp(t)為調和 分析之天文潮位；ηf(t)為氣象潮。在本計 畫中，假設氣象潮在預報未來的時刻內均 為已知定值ηf(t)，該值可由預報當時河口 之觀測水位與調和分析天文潮位之差值求 得。分述如下： 1.天文潮 利用調和分析分離出不同之週期分潮 對潮差之影響，並採用其中影響潮差較重 要之九個分潮作為調和分析計算之基準， 此九個分潮分別為M2、S2、N2、K1、Sa、 O1、K2、P1及M4潮，其方程式為： ηp η k ωk θk k m t a t ( )= + cos( − ) =∑
1 (39) 式中，ηp(t)為綜合潮之水位；η為觀測平 均水位；ωk為分潮之角速度；a 為調和k 分析所得之分潮振幅；θk為調和分析所得 之相位角；m為所使用之分潮數目(m=9)。 本計畫為配合模式之順利進行，蒐集民國 92年1~6月份之河口逐時潮位資料用以進 行天文潮之調和分析。 2.氣象潮 在颱洪期間若利用調和分析進行潮位 預報時，須納入氣象潮之影響，但氣象潮 之變因頗多，若完全納入考慮於實行上有 諸多困難。故目前之作法為將預報當時調 和分析天文潮位與觀測潮位之差值視為氣 象潮之影響，並在每次預報過程中將氣象 潮之影響視為此定值，如(38)式所述。 五、資料之整理與輸入 5.1 地文資料 1.蘭陽溪河系各河道斷面資料 在河川洪水演算過程中，河道斷面的形狀影響模擬結果甚大，為求符合實際現 況並提高模擬精度，本計畫採用水利署第 一河川局及宜蘭縣政府於民國 74、85 和 92年所測的蘭陽溪河系最新斷面資料，包 括蘭陽溪、羅東溪與宜蘭河，合計共 146 個斷面，其河道斷面縱剖面圖如圖2~圖4 所示。 2.堤防、抽水站及水工構造物資料 洪水模式演算過程中，需要輸入各斷 面的堤防高程資料，本計畫乃依據經濟部 水利署經由現地量測所得資料輸入數值模 式內，並假設在數值模式的演算過程中無 潰堤之情況發生，其堤防高程如圖2~圖4 所示。而流域內護岸及橋樑等水工構造物 則詳列如表1所示。至於抽水站之詳細資 料則尚待蒐集。 3.河床糙度 在進行河川洪水演算時，必須先將各 斷面之曼寧粗糙係數輸入模式中以進行數 值模擬。本模式目前所使用之曼寧粗糙係 數係依據經濟部水利署『蘭陽溪治理規劃 報告』報告中之分析資料，其中在下游河 道斷面之曼寧係數約為0.030~0.035、中游 者約為 0.035~0.038、上游者則為 0.040。 此外，本計畫在有橋樑的斷面處將曼寧粗 糙係數酌予提高，以反應水流受橋墩阻礙 之影響，各河道曼寧粗糙係數值如表2所 示，往後將經研究後再依據實測資料予以 檢定。 5.2 水文資料 1.邊界條件 本模式所需之邊界條件包括各河段上 游之入流歷線、各殘流域側入流及抽水站 等邊界條件。其中上游之邊界條件需要蘭 陽溪之家源橋、羅東溪之鼻頭橋與宜蘭河 之員山橋等3個水文站之流量歷線，故本 模式必須與子計畫一之降雨逕流模式銜 接。各殘流域之降雨逕流亦依子計畫一之 殘流域劃分，將其模擬結果與蘭陽溪河系 下游沿岸抽水站之最大抽水量，以側入流 方式納入本模式演算。在第一年度之模擬 初期考慮如表3所示之2、10和50年重現 期之逕流量；而側入流邊界條件則暫不予 以考量。本年度在各殘流域之降雨逕流方 面，則配合子計畫一之殘流域劃分，將其 模擬結果，以側入流方式納入本模式演 算。目前在上游邊界條件與殘流域側入流 之處理，已與子計畫一協調整合並完成模 式之銜接與初步測試。 2.下游河口潮位歷線 本子計畫已利用調和分析法發展蘭陽 溪河口潮位預估模式，未來將再配合中央 氣象局的實測資料以校正下游河口的邊界 條件。目前暫依據民國 93 年 1~6 月蘭興 橋分析結果的平均潮位0.24公尺，作為初 期數值模擬之下游邊界。 此外，本子計畫於第二年度研究期間 經由水利署第一河川局的協助，已收集到 興蘭橋潮位站 1985~1998 年紙本記錄及 1999~2004電子檔潮位資料。 六、初步研究結果 6.1 重現期洪水 本子計畫所建構之河川洪水演算模式 的模擬上游邊界包括：蘭陽溪上游之家源 橋、羅東溪上游之鼻頭橋與宜蘭河上游之 員山橋，下游邊界為河口。至於作為參數 校正的水文站則分別為蘭陽溪的牛鬥橋及 蘭陽大橋、宜蘭河的西門橋。其相對位置 示意圖見圖5所示。 經採用2、10和50年重現期之逕流量 及河口平均潮位0.24公尺等邊界條件，其 初步模擬於蘭陽溪、宜蘭河及羅東溪之水 位剖面線繪如圖6~圖8所示。因第一階段 研究主要在於模式架構之建立及程式之測 試，經由2、10和50年重現期洪水之測試 後，證實以運動波銜接動力波之河川演算 模式，可應用於蘭陽河系。
6.2 河口潮位 本子計畫於第一年度採用92年01月 01日至92 年06月30 日之興蘭橋潮水位 站潮位資料進行調和分析，所得各分潮之 參數如表 4。經由此法所完成之潮位預報 模式，進行模擬預報 92 年 07 月 01 日至 92年09月20日之河口潮位，潮位預報之 結果繪如圖 9~圖 11 所示，圖中顯示在 7 月份及8月份的潮位預報值與觀測潮位甚 為近似，顯示該潮位預報模式對於潮位變 化之趨勢確可掌握。但在9月份時因為儀 器故障導致實測資料部份缺失，因而降低 了潮位預測精度。 本子計畫於第二年度持續蒐集了 93 年01月01日至93年06月30日之興蘭橋 潮水位站潮位資料進行調和分析，所得各 分潮之參數如表5。並進行模擬預報93年 07 月 01 日至 93 年 09 月 20 日之河口潮 位，潮位預報之結果繪如圖 12~圖 14 所 示，圖中顯示潮位預報值與觀測潮位之趨 勢甚為一致。經由上述兩組測試，顯示該 潮位預報模式對於河口潮位變化之趨勢確 可充分掌握，若於汛期再配合氣象潮位之 影響，整體而言該模式確可提供符合現況 的潮位預報值。 6.3 歷史颱洪事件測試 本子計畫於第二年度，蒐集蘭陽溪歷 年颱洪期間的水位資料，茲彙整如表6所 示。由表6中顯示安珀颱洪事件的資料較 為齊全，故選用安珀颱洪事件進行測試。 本子計畫採用子計畫一(逕流模式)提供之 安珀颱洪(86 年)之上游邊界流量歷線及下 游河口潮位實測資料，進行河川演算模式 之模擬測試。 安珀(AMBER)颱風於民國 86 年 8 月 生成於呂宋島東方海面，於 8月 29日 03 時 50 分左右在花蓮秀姑巒溪北方附近登 陸後，強度減弱，繼續以西北方向前進， 於8月29日20時左右由馬祖附近進入大 陸，其路徑如圖15所示。 本子計畫選取該場颱風洪峰期間(86 年8月29日6時至8月29日24時)，進 行河川演算模式模擬測試，於牛鬥橋及蘭 陽橋之模擬水位結果結果如圖 16 及圖 17 所示，圖中顯示模擬結果均有高估之現 象，其中，牛鬥橋模擬水位與實測水位之 均方差(RMSE)為 2.65 公尺；蘭陽橋模擬 水位與實測水位之均方差(RMSE)為 1.97 公尺。 6.4 歷史颱洪事件整合預報測試 為求進一步測試整合模式，本階段主 要工作即針對安珀颱洪事件，就降雨預 報、地表逕流預報及河川洪水位預報進行 模式串接，在模式銜接測試前，並就各模 式之資料傳輸格式與其他子計畫進行檢討 並彙整，最後由總計畫訂定模式間傳輸檔 案之格式，以利模式順利銜接。 本階段中以降雨及水位資料較為完備 的安珀颱洪(86年8 月)進行各個模式串接 整合之預報測試。安珀颱洪之整合模擬預 報於牛鬥橋及蘭陽橋之預報水位結果如圖 18 及圖 19 所示。圖中顯示預報結果均有 高估之現象，究其原因包括：(1)河川水位 預報模式為整合模式之最下游端，因承接 上游端模式(降雨及逕流預報模式)之預報 誤差，造成誤差放大；(2)受限於歷史颱洪 事件實測資料短缺，河川水位預報模式之 參數檢定及驗證工作尚不完整，未來將持 續進行以求改善；(3)河川斷面資料尚待進 一步更新(目前蘭陽溪斷面資料為 85 年所 測，羅東溪斷面資料為74年之設計規劃斷 面，宜蘭河斷面資料最新，為92年所測)， 因河川斷面資料資料部份過於老舊恐不合 時宜且年份不一，故造成誤差；(4)河川水 位預報模式未來將即時觀測水位納入進行 回饋演算後，亦可提高預報精度。 七、結語與展望 7.1 結語 1.本計畫第一年度建立了蘭陽溪河川洪水
演算模式，並以2、10及50年重現期流 量進行模式的演算功能檢驗，經由模擬 過程顯示，以動力波理論為基礎所建立 的河川洪水演算模式確實能反映出蘭陽 溪河系於洪水來臨時之水位變化。 2.本計畫除了建立蘭陽溪河川洪水演算模 式之外，亦完成了修正預報初始值及參 數校正的理論推導，作為河川洪水演算 模式修改的依據。 3.第一年度因受限興蘭大橋於92年9月份 潮位站儀器故障及更換，使得 9 份月的 潮位預報結果較為不理想。故本子計畫 於第二年度蒐集93年1月至6月的潮位 觀測值，重新進行參數的檢定以符合流 域的現況。 4.為能使模擬結果能更符合現況，蘭陽溪 河系部份殘流域之側入流歷線，已由子 計畫一提供資料並進行整合性的模擬。 5.由安珀颱洪之模擬測試及預報測試結 果，可知各模式間在訂定傳輸檔案之統 一格式後，確可順利銜接並模擬及預報。 6.為能精確預測蘭陽溪河系的洪水流況， 應再蒐集過去蘭陽溪河系相關研究報 告，並將其結果進行分析與比較，以提 高河川洪水演算模式之完整性與精確 度，而達到防洪預報以減少洪災之功能。 7.2 展望 1.為能精確地作好洪水預報，應進一步以 河川洪水演算模式模擬其他數場颱洪事 件，加強各項參數之檢定與驗證。 2.為能提供更精確之洪水預報，河川洪水 演算模式內容尚需進行即時水位回饋演 算，並配合修正預報初始值及參數校正 等理論作進一步之修改，以符合洪水預 報之要求，亦為本計畫未來之主要研發 探討的工作。 3.為提昇潮位預報的精度，需再蒐集94年 1 月至 6 月的潮位資料進行調和參數的 分析。 4.由歷史颱洪模擬及預報結果顯示均有高 估之現象，其誤差原因如文中所述，未 來將持續蒐集相關資料，以提高預報精 度。 5.未來將加強與各子計畫間之資料整合， 如銜接子計畫一的逕流量作為模式上游 邊界，將子計畫二之研究成果予以調整 模式流況，配合子計畫四將其演算成果 予以展示，作為預報研究區域內重要地 點的洪水位及河川沿岸的水面剖線之參 考。 八、參考文獻 李鴻源、楊錦釧、葉克家、楊志達、謝慧 民，1996，瓣狀河系沖淤模式之發展 (四)-NETSTARS 模式使用者手冊，臺 大土木研究所研究報告第8501號。 茨木俊秀，福島雅夫，1997，最優化方法， 世界圖書出版公司。 許銘熙，1984，河川變量流非線性隱式法 模式之穩定性，臺大土木工程研究所 博士論文。 蔡長泰，1977，隱式法模擬變量流之應用 與比較，中華民國第一屆力學會議論 文集，285-307。 蔡長泰、歐善惠、顏沛華、宋長虹、陳聰 智，1984，濁水溪洪水預報系統之研 究(三)，行政院國科會防災科技研究報 告74-38。 蔡長泰、郭惠貞、蔡柏棋、蔡智恆，1991， 曾文溪洪水預報模式之研究(二)，行政 院國科會防災科技研究報告81-61號。 賴經都、王燦文、賴坤勇，1995，以多方 式特徵法建立多成分河流模式，行政 院 國 科 會 專 題 研 究 報 告 ，NSC 83-0410-E002-038。 許銘熙，1984，河川變量流非線性隱式法 模式之穩定性，臺大土木工程研究所 博士論文。 顏清連、王如意、許銘熙、楊德良、李天 浩，1996-1998，淡水河整體洪水預報 系統模式之研發第一至六階段報告， 台灣省水利處委託研究計畫，國立台 灣大學水工試驗所。
台灣省水利局，民國71年，宜蘭河治理計 畫，第一工程處。 台灣省水利局，民國83年，蘭陽溪治理規 劃報告，第一工程處。 經濟部水利處，民國89年，蘭陽溪治理規 劃報告，第一河川局。 Abbott, M. B., 1976, Computational Hydraulics, Pitman Publishing Ltd. London.
Ahmed, E. I., and Kenneth, V. H., 1979, Experimental investigation of flow over side weirs, Journal of Hydraulics Division, ASCE, 102, 1255-1367.
Amein, M., Fang, C. S., 1970. Implicit flood routing in natural channel, Journal of Hydraulics Division, ASCE, 96, 2481-2500.
Bajracharya, K., and Barry, D. A., 1997. Accuracy criteria for linear diffusion wave flood routing, Journal of Hydrology, 195, 200-217.
Camacho, L. A., Lees, M. J., 1999, Multilinear discrete lag-cascade model for channel routing, Journal of Hydrology. 226, 30-47.
Carmeron, D., Kneale, P., See, L., 2002. An evaluation of a traditional and neural net modeling approach to flood forecasting for an upland catchment, Hydrological Processes, 16, 1033-1046.
Chiu, C. L., Emmanuel, O. I., 1978. Kalman filtering in open channel flow estimation, Journal of the Hydraulics Division, 104(8), 1137-1151.
Chow, V. T., Maidment, D. R., Mays, L. W., 1988. Applied Hydrology, McGraw-Hill Book Company.
Chow, V. T., 1973. Open-Channel Hydraulics, McGraw-Hill Book Company.
Crissman, R. D., and Chiu, C. L., 1993. Uncertainties in flow modeling and forecasting for Niagara River, Journal of Hydraulic Engineering, ASCE, 119(11), 1231-1250.
Einhelling, R., Holly, F. M., Hsu, S. H., Schwarz, P., Schaefer J., and Yang, J. C., 1990. Charima-numerical simulation of unsteady water and sediment movement in multiply connected network of
mobile-bed channels. IHR Report No. 343, Iowa institute of Hydraulic Research, The University of Iowa, Iowa City, Iowa 52242 USA.
Englund, F., and Fredsoe, J., 1976. A sediment transport model for straight alluvial channels. Nordic Hydrology, 7, 296-306.
Fread, D. L., 1992. Flood routing models and Manning n, Channel Flow
Resistance：Centennial of Manning's
Formula, 421-435.
Hsu, M. H., Fu, J. C., and Liu, W. C. 2003. Flood routing with real-time stage correction method for flash flood forecasting in the Tanshui River, Taiwan. Journal of Hydrology, 283, 267-280. Jarrett, R. D., 1984. Hydraulics of
high-gradient streams, Journal of Hydraulic Engineering, 110(11), 1519-1539.
Jorge, N., and Stephen, J. W., 1999. Numerical optimization, McGraw-Hill Book Company.
Kothyari, U. C., and Singh, V. P., 1999. A multiple-input-output model for flow forecasting, Journal of Hydrology, 220, 12-26.
Kuchment, L. S., Nazarov, N. A., 1999. A model of transient flow and a procedure of short-term forecasting of discharges for the Lower Vyatka river basin, Water Resources, 26(4), 361-367.
Liu, W. C., Hsu, M. H., and Wang, C. F., 2003. Modeling of flow resistance in mangrove swamp at mouth of tidal Keelung river, Taiwan, Journal of Waterway, Port, Coastal and Ocean Engineering, ASCE, 129(2), 86-92.
Nguyen, Q. K., and Kawano, H., 1995. Simultaneous solution for flood routing in channel networks, Journal of Hydraulic Engineering, ASCE, 121(10), 744-749.
Mahmood, K., and Yevievich, V., 1975. Unsteady flow in open channels, Vol. I & II , Water Resources Publications, Inc., Fort Collins, Colorada.
Parker, G., 1991. Selective sorting and abrasion of river gravel, Journal of Hydraulic Engineering, ASCE, 117(2),
131-149.
Schaffranek, R. W., 1987a. Flow model for open-channel reach or network, U.S. Geological Survey Professional Paper 1384, 12.
Schaffranek, R. W., 1987b. A flow simulation model of the tidal Potomac River, U.S. Geological Survey Water-Supply Paper 2234-D, 41.
Sepaskhah, A. R., and Bondar, H., 2002. Estimation of Manning roughness coefficient for bare and vegetated furrow,
Wiley, J. B., 1993. Simulation flow and solute transport and mitigation of a hypothetical soluble-contaminant spill for the new river in the new river gorge national river, U.S. Geological Survey Water-Resources Investigations Report 93-4105, 39.
Wong, T. S. W., and Zhou, M. C., 2003. Kinematic wave parameters and time of
travel in circular channel revisited, Advances in Water Resources, 26, 417-425.
Wormleaton, P. R., and Karmegam, M., 1983. Parameter optimization in flood routing, Journal of Hydraulic Engineering, 110(12), 1799-1814.
Yen, C. L., 1985. An overview on flood routing models, Proceedings of the ROC-Japan Joint Seminar on Multiple Hazards Mitigation National Taiwan University, Taipei, Taiwan, 559-582. Yen, C. L., Hsu, M. H., 1982. A numerical
model for unsteady flow in river system, Proceedings of the thirds congress of the Asian and Pacific regional division of the International Association for Hydraulic Research, Bandung, Indonesia, 91-101.
表1 蘭陽溪河系沿岸之堤防、護岸及橋樑 左岸 右岸 河川 堤防名稱 護岸名稱 堤防名稱 護岸名稱 橋樑名稱 狀圍 大洲 中溪洲 家源橋 員山 中溪洲 紅柴林 牛鬥橋 再連 三星二號 蘭陽大橋 常嶺 三星舊堤 興蘭大橋 崙埤 三星一號 廣興橋 松羅 破布烏 清洲橋 牛鬥 東疊 碼崙 2 號 碼崙 碼崙 1 號 蘭 陽 溪 五結 柯子林 寒溪左岸護岸 廣興 阿里史 鼻子頭 小南澳 小埤 淋漓 寒溪一號 羅 東 溪 寒溪 後防潮堤 新生護岸 狀圍防潮堤 興蘭大橋 公館 狀圍 狀圍大橋 梅州 金六結 中央大橋 永和 七張橋 慈安橋 鐵路橋 宜蘭橋 中山橋 宜 蘭 河 員山大橋 表2 蘭陽溪河系各斷面曼寧系數值 斷面之曼寧 n 值採用表 蘭陽溪 羅東溪 代表斷面 採用值 代表斷面 採用值 1~12 0.030 1~3 0.030 13~15 0.030 4~5 0.030 16~18 0.033 6~7 0.030 19 0.033 8~9 0.031 20~21 0.035 10~11 0.031 22~25 0.035 12~14 0.031 26~27 0.035 15 0.031 28~29 0.035 16 0.032 30~32 0.035 17~18 0.032 33~34 0.035 19~24 0.033 35~36 0.035 37~40 0.035 41~51 0.035 52~57 0.035 58~59 0.038 60~61 0.038 62~63 0.038 64~65 0.038 66~67 0.038 68~69 0.038 70~73 0.038 74~75 0.040 76~77 0.040
表3 蘭陽溪河系各重現期洪水量(單位：cms) 重現期距 2 5 10 20 25 50 100 200 蘭陽溪部分 蘭陽大橋 3100 4700 5700 6600 6900 7700 8500 9300 羅東溪合流前 2600 4000 4800 5600 5800 6500 7200 7800 清水溪合流前 2100 3400 4200 4900 5100 5800 6400 7000 羅東溪流域(出口) 700 1000 1200 1300 1400 1500 1600 1800 宜蘭河部分 大礁溪匯流點 330 450 530 640 --- 800 --- --- 中山橋 660 930 1100 1310 --- 1610 --- --- 河口 660 930 1100 1310 --- 1610 --- --- 表4 蘭陽溪河口興蘭橋水位九個分潮之調和分析參數值 分潮編號 分潮名稱 分潮週期(時) 分潮振幅(公分) 相位角(度) 1 M2 12.42 23.064 136.194 2 S2 12.00 9.532 -125.633 3 N2 12.66 3.824 -176.20 4 K1 23.93 14.618 -128.835 5 Sa 8765.28 9.599 -155.728 6 O1 25.82 12.819 149.937 7 K2 11.97 5.883 35.335 8 P1 24.07 4.138 -141.092 9 M4 6.21 1.562 -112.866 Mean(平均水位高程) ：23.505 公分 RMS(殘餘水位差) ：14.905 公分 ER(潮波能量比) ：72.32% 註：92年1~6月資料分析結果 表5 蘭陽溪河口興蘭橋水位九個分潮之調和分析參數值 分潮編號 分潮名稱 分潮週期(時) 分潮振幅(公分) 相位角(度) 1 M2 12.42 26.471 40.25 2 S2 12.00 11.691 -131.03 3 N2 12.66 5.923 171.75 4 K1 23.93 17.528 -129.93 5 Sa 8765.28 17.763 -137.06 6 O1 25.82 13.683 50.85 7 K2 11.97 1.654 38.83 8 P1 24.07 4.4314 37.99 9 M4 6.21 4.1665 -124.30 Mean(平均水位高程) ：27.803 RMS(殘餘水位差) ：10.439 ER(潮波能量比) ：87.89% 註：93年1~6月資料分析結果
表6 蘭陽溪流域歷史颱洪水位資料記錄調查統計 水位站 水為站名稱 颱風名稱 家源橋 牛鬥橋 蘭陽大橋 西門橋 興蘭大橋 ELLIE,愛麗 --- 3 3 --- --- NAT,耐特 --- 3 3 3 --- NAT,耐特 --- --- --- --- --- RUTH,露絲 --- --- 3 3 --- POLLY,寶莉 --- --- --- --- --- TED,泰德 --- --- --- --- --- TIM,提姆 --- --- 3 3 --- DOUG,道格 --- --- 3 3 --- GLADYS,葛拉絲 --- --- 3 3 --- SETH,席斯 --- --- 3 3 --- HERB,賀伯 3 3 3 3 3 GLORIA,葛樂禮 --- 3 3 --- 3 WINNIE,溫妮 --- --- 3 3 3 AMBER,安珀 3 3 3 3 3 YANNI,楊妮 --- --- 3 3 3 註：1.---表無水位資料。2.3有水位資料。 以 Newton-Raphson 疊代法 求解各斷面之流況 配合邊界預報值 計算河川未來流況 否 否 是 是 是否接收到 即時觀測水位 結束 演算開始 輸入資料 1.即時邊界觀測值 2.前一時刻之流況 計算 水深差平方和 是否小於容許值 重新計算 各斷面之即時水位與流量 否 否 是 是 預報初始值設定 使用最小平方法 即時參數修正 使用最佳化理論 水深差平方和 是否符合要求精度? 否 否 是 是 圖1 參數即時校正河川洪水演算流程圖 -50 0 50 100 150 200 250 300 350 400 450 0.0 5000.0 10000.0 15000.0 20000.0 25000.0 30000.0 35000.0 40000.0 45000.0 距離(m) 高程( m ) 底床高程 左岸堤防高 右岸堤防高 -5.00 0.00 5.00 10.00 15.00 -200 0 200 400 600 8001000 1200 1400 河口 206.00 208.00 210.00 212.00 214.00 216.00 218.00 220.00 222.00 224.00 0 50 100 150 200 250 300 牛鬥橋 家源橋 蘭陽橋 堤防高 底床高程 堤防高 底床高程 圖2 蘭陽溪縱斷面剖面圖 -5 0 5 10 15 20 25 0 1000 2000 3000 4000 5000 6000 7000 距離(m) 高程( m ) 左岸堤高 右岸堤高 底床高程 0.00 2.00 4.00 6.00 8.00 10.00 12.00 0.00 50.00 100.00 150.00 200.00250.00 300.00 西門橋 員山橋 堤防高 底床高程 堤防高 底床高程 圖3 宜蘭河縱斷面剖面圖 0 10 20 30 40 50 60 70 0 1000 2000 3000 4000 5000 6000 7000 8000 距離(m) 高程 (m ) 左岸堤高 右岸堤高 底床高程 51.00 52.00 53.00 54.00 55.00 56.00 57.00 58.00 59.00 0.00 50.00 100.00 150.00 200.00 250.00 300.00 350.00 鼻頭橋 堤防高 底床高程 0 10 20 30 40 50 60 70 0 1000 2000 3000 4000 5000 6000 7000 8000 距離(m) 高程 (m ) 左岸堤高 右岸堤高 底床高程 51.00 52.00 53.00 54.00 55.00 56.00 57.00 58.00 59.00 0.00 50.00 100.00 150.00 200.00 250.00 300.00 350.00 鼻頭橋 0 10 20 30 40 50 60 70 0 1000 2000 3000 4000 5000 6000 7000 8000 距離(m) 高程 (m ) 左岸堤高 右岸堤高 底床高程 51.00 52.00 53.00 54.00 55.00 56.00 57.00 58.00 59.00 0.00 50.00 100.00 150.00 200.00 250.00 300.00 350.00 鼻頭橋 堤防高 底床高程 堤防高 底床高程 圖4 羅東溪縱斷面剖面圖
上游邊界 回饋站 下游邊界 家源橋 興蘭橋 蘭陽溪 宜蘭河 牛鬥橋 蘭陽橋 鼻頭橋 員山橋 中山橋 羅東溪 上游邊界 回饋站 下游邊界 家源橋 興蘭橋 蘭陽溪 宜蘭河 牛鬥橋 蘭陽橋 鼻頭橋 員山橋 中山橋 羅東溪 圖5 蘭陽溪河川洪水演算模式示意圖 -50 0 50 100 150 200 250 0 5000 10000 15000 20000 25000 30000 35000 離河口距離(m) 高程 (m) 2年 10年 50年 堤防高 圖6 蘭陽溪各重現期洪水位剖面圖 0 5 10 15 20 25 0 5000 10000 15000 20000 離河口距離(m) 高程 (m) 2年 10年 50年 堤防高 圖7 宜蘭河各重現期洪水位剖面圖 0 5 10 15 20 25 0 500 1000 1500 2000 2500 3000 3500 離河口距離(m) 高程 (m ) 2年 10年 50年 堤防高 圖8 羅東溪各重現期洪水位剖面圖 -60 -40 -20 0 20 40 60 80 100 120 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 日 潮位 (c m ) 預測 實測 圖9 92年7月潮位歷線比較圖 -40 -20 0 20 40 60 80 100 120 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 日 潮位 (c m ) 預測 實測 圖10 92年8月潮位歷線比較圖 -40 -20 0 20 40 60 80 100 120 140 160 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 日 潮位( cm ) 預測 實測 圖11 92年9月潮位歷線比較圖 七月潮位 -60 -40 -20 0 20 40 60 80 100 120 140 日 潮位 (c m ) 預測1 實測 圖12 93年7月潮位歷線比較圖
八月潮位 -100 -50 0 50 100 150 200 日 潮位( cm ) 預測1 實測 圖13 93年8月潮位歷線比較圖 九月潮位 -60 -40 -20 0 20 40 60 80 100 120 140 日 潮位 (c m ) 預測1 實測 圖14 93年9月潮位歷線比較圖 圖15 安珀颱風路徑圖 牛鬥橋之模擬結果 207.00 208.00 209.00 210.00 211.00 212.00 213.00 214.00 215.00 216.00 217.00 0 2 4 6 8 10 12 14 16 18 20 時間（小時） 水 位︵ 公 尺︶ 觀測值 模擬值 圖16 安珀颱風牛鬥橋模擬水位 蘭陽橋之模擬結果 0.00 1.00 2.00 3.00 4.00 5.00 6.00 7.00 8.00 0 2 4 6 8 10 12 14 16 18 20 時間(小時) 水 位 ( 公 尺 ) 觀測水位 模擬水位 圖17 安珀颱風蘭陽橋模擬水位 205.00 207.00 209.00 211.00 213.00 215.00 217.00 06: 00 07: 00 08: 00 09: 00 10: 00 11: 00 12: 00 13: 00 14: 00 15: 00 16: 00 17: 00 18: 00 19: 00 20: 00 21: 00 22: 00 23: 00 00: 00 8月29日 8月30日 水 位 ( 公 尺 ) 實測水位 預報 1 預報 2 預報 3 預報 4 預報 5 預報 6 預報 7 預報 8 預報 9 預報 10 預報 11 預報 12 預報13 預報14 圖18 安珀颱風牛鬥橋預報水位 -0.50 0.50 1.50 2.50 3.50 4.50 5.50 6.50 7.50 8.50 06 :00 07 :00 08 :00 09 :00 10 :00 11 :00 12 :00 13 :00 14 :00 15 :00 16 :00 17 :00 18 :00 19 :00 20 :00 21 :00 22 :00 23 :00 00 :00 8月29日 8月30日 水 位 ( 公 尺 ) 實測水位 預報 1 預報 2 預報 3 預報 4 預報 5 預報 6 預報 7 預報 8 預報 9 預報 10 預報 11 預報 12 預報13 預報14 圖19 安珀颱風蘭陽橋預報水位 附錄一 敏感係數 假設河道中無側入流及溢流，則可將 水面剖線表示成 2 1 F S S dx dY o f − − = (A1-1) 式中，F為福祿數。 藉由兩個相鄰水位站i 和(i+1)，則可 以將(A1-1)表示為 i i i Y Y = +1−ψ (A1-2)
(
)
(
)
∆ ⋅ + − + − = + + i i i i f i f o i X F F S S S 2 1 2 1 2 1 1 2 1 ψ (A1-3) 式中，ψi 為兩水文站的水深差微分式。當 j > i時，經由連鎖律可將(A1-2) 對 曼寧係數的偏微分寫成 ( ) i i j i i i j i i i j i i i j i j i i i j j i Y n Y Y n Y Y n Y Y n Y n Y Y n n Y ∂ ∂ ∂ ∂ − ∂ ∂ − ∂ ∂ = ∂ ∂ ∂ ∂ − ∂ ∂ ∂ ∂ − ∂ ∂ = − ∂ ∂ = ∂ ∂ + + + + + + ψ ψ ψ ψ ψ 1 1 1 1 1 1 1 (A1-4) (A1-4)式中的微分項 1 + ∂ ∂ i i Y ψ 和 i i Y ∂ ∂ψ 可 被改寫為 i j i j i D n Y n Y ∂ ∂ = ∂ ∂ +1 (A1-5) 而Di為 ∂ ∂ + ∂ ∂ − = + i i i i i Y Y D ψ ψ 1 1 1 (A1-6) 同 樣 的 ， 可 以 將 j i n Y ∂ ∂ +1 表 示 成 具 有 j i n Y ∂ ∂ +2 的形式 i i j i i j i D D n Y D n Y ∂ ∂ = ∂ ∂ + + + 1 2 1 (A1-7) 重覆上述步驟，則 i i j j j j j i D D D D n Y n Y 1 2 1 − + − ∂ ∂ = ∂ ∂ Λ (A1-8) 將(18)代入(A1-8)則
