行政院國家科學委員會專題研究計畫 成果報告
山崩及土石流災害之探討--子計畫:以碎形維度理論結合數
位影像分析應用於土石流潛勢之判定(III)
研究成果報告(完整版)
計 畫 類 別 : 整合型 計 畫 編 號 : NSC 95-2625-Z-002-012- 執 行 期 間 : 95 年 08 月 01 日至 96 年 07 月 31 日 執 行 單 位 : 國立臺灣大學土木工程學系暨研究所 計 畫 主 持 人 : 陳榮河 計畫參與人員: 博士班研究生-兼任助理:黃永鈴、黃永鳳 碩士班研究生-兼任助理:紀泰安、許珮甄 處 理 方 式 : 本計畫可公開查詢 中 華 民 國 96 年 10 月 31 日I
目錄及頁次
目錄及頁次...I 表目錄...IV 圖目錄...V 第一章 緒論...1 1.1 研究動機...1 1.2 研究目的...2 第二章 文獻回顧...5 2.1 土石流發生機制...5 2.2 孔隙水壓引發型土石流...7 2.3 降雨引發型土石流...9 2.4 水流入滲之效應...10 第三章 碎形維度理論與相關推導...13 3.1 碎形理論...13 3.2 方格覆蓋維度...14 3.3 現地與篩分析試驗之碎形維度差異...19 3.4 土壤保水曲線碎形模式...19 3.5 案例驗證...23 3.6 滲透曲線碎形模式...29 3.6.1 理論推導...29 3.6.2 滲透曲線碎形模式探討...32 第四章 影像於篩分析之應用...34II 4.1 二維影像分析之研究...34 4.2 數位影像理論...36 4.2.1 影像分割...36 4.2.2 分水嶺分割法...38 4.2.3 門檻值法...43 4.3 二維影像分析模型試驗...48 第五章 現地採樣及粒徑影像分析...55 5.1 研究區域地理位置...55 5.2 研究區域地形概況...57 5.3 研究區域地質概況...62 5.4 粒徑影像分析...66 第六章 單向度水分傳輸實驗...69 6.1 實驗設備...70 6.2 雨量校正試驗...76 6.3 實驗方法與步驟...80 6.4 實驗結果與討論...83 6.4.1 No.1實驗...83 6.4.2 No.2實驗...86 6.4.3 No.3實驗...88 6.4.4 No.4實驗...90 6.4.5 No.5實驗...92 6.4.6 綜合討論...95 6.5孔隙水壓傳遞機制模式之應用...96 第七章 具體完成項目...99
III 7.1 現地採樣...99 7.2 完整粒徑分佈曲線繪製...101 7.3 物理模型試驗進行...103 7.4 土石流發生孔隙水壓傳遞機制模式建立...104 7.5 成果報告彙整...105 參考文獻...107
IV
表目錄
表 3.1 估算保水曲線之基本參數...29 表 4.1 各種體積估算結果比較...53 表 5.1 台灣北部漸新世地層分層與命名表...68 表 6.1 水頭高與水量之關係...78 表 6.2 供給水量與模擬降雨強度之關係...79 表 6.3 本實驗之材料性質...80V
圖目錄
圖1.1 研究流程圖...4 圖2.1 天然土石壩潰決的三種型態...7 圖2.2 雨量記錄及在nest3 處所量測孔隙水壓紀錄...9 圖2.3 土層內部水分累積示意圖...12 圖3.1 在不同網格尺規下同一顆粒之覆蓋方格數示意圖...15 圖3.2 方格覆蓋法邏輯與篩分析過程之比對...16 圖3.3 篩分析試驗各號篩停篩結果示意圖...16 圖3.4 方格覆蓋法計算過程示意圖...18 圖3.5 碎形毛細管模式示意圖...26 圖3.6 雙對數關係表示粒徑分佈曲線示意圖...27 圖3.7 篩分析曲線雙對數圖...27 圖3.8 改良篩分析曲線對數圖...28 圖3.9 保水曲線預測值與現場實驗值比較圖...28 圖3.10 不同顆粒分布下相關水力滲透係數與水力勢能之關係圖...33 圖3.11 不同顆粒分布相關水力滲透係數與正規化體積含水率關係圖..33 圖4.1 海砂粒徑分佈之累積頻率曲線...35 圖4.2 義大利Ravaneto 區域影像拍攝位置之斷面全景...35 圖4.3 影像分析照片與再製圖...36 圖4.4 一灰階影像計算形態學梯度影像過程...39 圖4.5 形態學梯度運算處理過程示意圖...40 圖4.6 流域、極小區域、分水嶺線概念圖...41 圖4.7 流域築堤概念圖...42VI 圖4.8 門檻值法之灰階直方圖...45 圖4.9 彩色影像轉灰階影像...49 圖4.10 直接進行二值化運算...49 圖4.11 分水嶺演算法...50 圖4.12 分水嶺演算後經二值化處理...50 圖4.13 標註個別顆粒物體...51 圖4.14 分析個別顆粒物體...51 圖4.15 體積幾何示意圖...52 圖5.1 北橫公路復興至巴陵段研究區地理位置圖...56 圖5.2 北橫公路復興至巴陵段研究區域地形圖...58 圖5.3 研究區域一之土石流上游全貌...59 圖5.4 研究區域一之二處主要源頭之近照...60 圖5.5 研究區域一採樣區域照片...60 圖5.6 研究區域二採樣區域照片...61 圖5.7 研究區域三採樣區域照片...61 圖5.8 台灣地質分區圖...62 圖5.9 雪山山脈帶與西部麓山帶和脊樑山脈帶分界圖...63 圖5.10 北橫公路復興至巴陵段研究區域地質圖...65 圖5.11 土石流材料...66 圖5.12 影像分析二值化之結果...67 圖5.13 IPTK分析所得之顆粒編號之標註...67 圖5.14 影像分析與篩分析之比較...68 圖6.1 模型試驗配置示意圖...69 圖6.2 實驗設置全景...70
VII 圖6.3 上部供水設施近照...73 圖6.4 微型套管套裝不銹鋼針頭...73 圖6.5 不銹鋼針頭形成之均勻降雨...74 圖6.6 透過壓克力裝置觀察模擬降雨情況...74 圖6.7 下部供水設施近照...75 圖6.8 水量校正實驗...77 圖6.9 水量與水頭高...78 圖6.10 入滲率與水頭高示意圖...79 圖6.11 本研究所使用之材料粒徑分布曲線...81 圖6.12 本研究所使用材料之碎形維度關係圖...81 圖6.13 No.1實驗之基質吸力與深度之關係圖...85 圖6.14 No.1實驗之體積含水率與深度之關係圖...85 圖6.15 No.2實驗之基質吸力與深度之關係圖...87 圖6.16 No.2實驗之體積含水率與深度之關係圖...87 圖6.17 No.3實驗之基質吸力與深度之關係圖...89 圖6.18 No.3實驗之體積含水率與深度之關係圖...89 圖6.19 No.4實驗之基質吸力與深度之關係圖...91 圖6.20 No.4實驗之體積含水率與深度之關係圖...92 圖6.21 No.5實驗之基質吸力與深度之關係圖...94 圖6.22 No.5實驗之體積含水率與深度之關係圖...94 圖7.1 北橫公路復興至巴陵段沿途採樣地點...100 圖7.2 影像分析結果示意圖...101 圖7.3 本研究所使用之材料粒徑分布曲線與所使用材料之碎形維度關係 圖...103
VIII
圖7.4 No.4實驗之基質吸力與深度及體積含水率與深度之關係圖...103 圖7.5 No.5實驗之基質吸力與深度及體積含水率與深度之關係圖...103 圖7.6 利用滲透曲線碎形模式得到不同碎形維度之相關水力滲透係數與 水力勢能及正規化體積含水率關係圖...105
1
第一章
緒論
1.1 研究目的
台灣位於板塊交界帶,地質構造破碎、地震頻繁。此外降雨多集 中於颱風季,且以暴雨型態出現,因此地滑、崩坍、土石流等坡地災 害頻傳。近十年來因土石流造成的經濟損失及人員傷亡尤其慘重。如 1996 年 7 月賀伯颱風,造成南投縣陳有蘭溪流域發生嚴重土石流災 害;1999 年 921 集集大地震,使得土石流危險溪流增至 721 條;又 2001 年 7 月桃芝颱風,土石流危險溪流暴增至 1420 條;2004 年 7 月 敏督利颱風,造成中部及南部地區重大損失。根據統計,台灣地區從 民國70 年到 89 年,因土石流災害造成的死亡人數約 200 人,在鄰國 日本1967 到 1987 年有 1257 人死於土石流災害。在土石流發生的區 域,常見房屋淹沒、橋樑沖毀與交通中斷,造成人民生命財產嚴重的 損失。因此預測土石流發生之可能性,可以進一步預防土石流的發 生,提供給下游居民預先警告,以避免及減少生命財產的損失。 早期對發生土石流預測,往往只利用降雨與破壞紀錄,如Canon & Ellen (1985)以美國舊金山灣區 6 次的降雨量與災害紀錄,建立該地 區土石流發生的臨界降雨基準,並作為防災警報發布的依據。謝正倫 與張東炯(1996)、范正成與姚正松(1997)等亦提出以降雨強度與2 累積雨量為指標之方法。而近年來相關的研究,則將地文因子加入考 慮,亦即對分析對象進行分類,除提升分析的精度外,並藉以對災害 可能發生的時間與地點一併作較佳的預測。 土石流發生破壞機制之研究方面,目前國內外學者多考慮水可位 於土層內任何位置引發土石流之可能性;也有考慮土體內部滲流力作 用及超額孔隙水壓引發土石流之可能性(周憲德與廖偉民,1998),雖 所考慮之機制已相當完善,但有關入滲效應所產生之孔隙水壓變化對 土石流發生機制影響,尚缺深入探討。Tarantino & Bosco (2000)曾提 出入滲對非飽和邊坡穩定之影響,但卻未對濕潤線(wet front)以上 孔隙水壓變化進行討論;在土石流發生臨界雨量線方面,均以統計方 法分析土石流發生之降雨特性參數,如降雨強度、降雨延時、累積雨 量及臨前降雨量等四項 (詹錢登,2000),這些方法多不考慮現場地 文因子條件。因此,如何同時考慮地文及水文因子為本研究主要研究 內容。
1.2 研究方法
有鑒於近年來土石流災害頻傳,嚴重影響人民生命安全、財產, 因此瞭解土石流發生機制、土石流危害度分級與土石流災害防制乃為 刻不容緩的工作。而土石流發生時,往往具有流速快、衝擊力大、發3 生時間短等特性,而且常伴隨著豪雨、崩塌材料而發生,在此情況下, 欲瞭解土石流之整體特性確有其困難,所以學者紛紛著手研究土石流 發生的力學機制與影響因素,其最終目的旨在對土石流災害的成因及 行為能夠有清楚的認識,並進而建立有效的防治對策。 研究土壤水分的傳輸現象是十分複雜的,若要探討土體內孔隙水 壓,必須先從其傳遞機制著手。而孔隙之多寡與分佈形態,直接控制 其滲透行為,且孔隙之型態或連通方式又受顆粒粒徑分佈範圍與顆粒 排列方式之影響。其中,土壤內部在平衡時所保存的含水量與基質吸 力大小有關,其關係曲線即所謂的土壤保水曲線(soil moisture retention curve),以下簡稱『保水曲線』。而不同土質或不同初始孔隙 率狀況下,保水曲線各不相同。甚至具同一孔隙率之土體亦可能有多 種不同之顆粒排列方式,若僅依孔隙率將無法反應真實情形。因此建 立保水曲線之預測模式,將可進一步了解土壤顆粒內水流現象,並對 了解土石流材料性質有所助益。 另一方面,土壤之篩分析可以提供一些土壤基本物性及級配情況 之資料,對於土石流水分傳輸現象之完整描述亦具關鍵性的影響。然 而,土石流材料粒徑差異非常大。加上土石流發生地一般位於山區, 不易到達,不論是採樣或是現地直接進行篩分析試驗均非常困難。故 本研究嘗試以數位影像紀錄現地土石流材料粒徑分佈情況,冀望能解
4 決目前所遭遇之困難。
5
第二章
文獻回顧
本章對本研究相關部份作介紹,主要針對土石流發生影響因 素,作一概略性描述,再就國內外各學者之相關研究文獻回顧整理, 以瞭解前人研究之成果,並作為本研究之理論基礎。2.1 土石流發生機制
Takahashi (1991)認為土石流發生之原因有三種:(1)滑動塊體之液 化。(2)崩塌形成之土石壩潰決。(3)溝床材料之起動或因逕流引起沉積顆粒流動之現象。Anderson & Sitar (1995)則認為有另一種發生之情
況,為地下水壓過大,造成土體內孔隙水壓變大,使土壤有效應力降 低引起之破壞。綜合上述,將土石流發生機制分為下列四種現象: (引 自陳榮河(1999)土石流之發生機制)。 a. 滑動型土石流 土塊滑動初期,除了滑動面附近之位移量較大以外,滑動體內 之變形量並不大。然若持續滑動,塊體內之變形量逐漸增大,使得整 體結構逐漸瓦解,若有足夠的水量進入孔隙中,使土體飽和並發生液 化現象,將形成土石流。Takahashi(1978)由理論推測,若邊坡坡度大 於20 度,且滑動距離超過塊體厚度數十倍之距離,並有充足之水量,
6 通常塊體會發生液化之現象。 b. 潰壩型土石流 由滑動造成之天然土石壩,其潰決常急劇發生於急短時間內, 這種破壞與某些土石流發生前,其溝谷流量突然停止或減少之現象吻 合。由高橋保(1997)試驗之觀察分為: (1)溢流沖蝕型(圖 2.1(a))-發生於土體透水性低,但剪力強度高 且有大水量供應之情況,因土體之透水性低,故上游水位面上升之速 度遠比入滲於土體內之速度快,且因土體強度高不易發生破壞,故快 速上升之水位,終使水流溢流,並沖蝕壩體表面。 (2)驟然崩潰型(圖 2.1(b))-發生於土體滲透性較上述類型大,且 剪力強度較弱之情況,因透水性較佳,故壩體內及上游之水位面均同 時上升,直至某一臨界水位時,因剪力強度不足,土體內部發生滑動。 滑動初始,崩落之土石因尚未完全飽和,故無法形成土石流,而堆積 於壩趾附近,然而隨即因滑動促成上游之水體突然受到釋放,而快速 沖過崩塌後之壩表面及僅短暫堆積於壩趾附近之崩落土石,而形成土 石流。 (3)漸進式破壞型(圖 2.1(c))-發生於土體滲透性高之情況,其特 徵為在上游水位面尚不高時,即可見滲流出現於壩下游面之低處。若 土體強度不高,將於滲流出處產生破壞,並向上漸次延伸,直至接近
7 上游水位面時,因水之推移力,使未滑動之剩餘部分發生破壞引發土 石流。 圖2.1 天然土石壩潰決的三種型態(高橋保,1997) c. 沖蝕型土石流 當土壤的剪力強度低於逕流產生之剪應力時,整個堆積層為不 穩定之狀態,若有急速大量的地表逕流沖蝕,將使堆積層不穩定引發 土石流動。
2.2 孔隙水壓引發型土石流
8
Sitar(1990), Ala & Mathewson (1990)等,均認為土體內之孔隙水壓 突然上升,是導致土體不穩定,而轉變成土石流的重要因素。傳統上, 一般假設雨水的垂直入滲(infiltration)是使土體飽和、孔隙水壓上升之 最重要的途徑,然而對於堆積於溪谷或邊坡上的崩積土層而言,雖有 構造疏鬆易於透水的內部,但其表面可能因風化或植物之生長而不易 透水,再者,Sidle & Milner (1990) 在調查過美國 Oregon 到
Southeastern Alaska 的海岸山脈可能發生土石流的區域之後,從孔隙
水壓記錄與雨量資料來研判, Silde 認為孔隙水壓激發上升的速度相
當快,很難完全以雨水垂直入滲來解釋。
除此之外,Johnson & Sitar (1990) 於 San Francisco Bay 地區 (1982,1983 及 1986 年)及 Ala & Mathewson (1990)等人於美國 Wasatch Front 地區(1983 及 1984 年)調查土石流災害,發現在災 後許多產生土石流之處,仍有相當多的水從出露的破碎底岩中流出, 這些滲流水有時持續好幾天甚至好幾個星期。且從其水壓紀錄發現 到,土石流發生時有正的超額孔隙水壓被激發,如圖2.2,初始時, 在深度為120 公分時的水壓均低於在深度 30 公分與 75 公分的水壓 值,而當降雨強度達某一強度後,會有兩處正的超額孔隙水壓尖峰值 出現於深度120 公分處,即水壓有突然被激發的現象,隨即有土石 流發生。且周憲德、廖偉民(1998)亦提出孔隙水壓會使土石流發生的
9
臨界坡度降低,造成土石流機率大增,亦即超額孔隙水壓使土石流在 緩坡時就可能會發生。
圖2.2 雨量記錄及在 nest3 處所量測孔隙水壓紀錄(Sitar et al., 1992)
2.3 降雨對土石流發生機制之影響
降雨會影響土體的含水量和入滲情形,進而影響土石流發生的 時間和規模大小。降雨特性與土石流發生關係之模式,其中較常用的 降雨特性參數為降雨強度I、降雨延時 T、累積雨量 R、以及前期降 雨量P 等。降雨量、降雨強度與降雨延時到底要多大才會發生土石 流?青木佑九(1980)研究日本地區 23 場降雨所造成的 46 場土石流 災害,提出日本地區土石流的降雨特性為: (1) 不考慮前期降雨的情況下,強度大的降雨延時數小時或普通降雨 延時12 小時以上,再加上持續 3~6 小時強度約 30~40 mm/hr 的10 降雨,即會發生土石流災害(此時累積雨量達 100 mm 以上)。 (2) 若累積降雨量在 150~200 mm 以上,即使小於 30~40 mm/hr 的降 雨也會發生土石流災害。 (3) 累積降雨達 400 mm 土石流災害以上,一定會發生。 降雨延時、降雨強度與累積雨量一直被視為是激發土石流的主要 條件,主要因為土石流發生地附近設有雨量站或氣象觀測站,降雨量 資料便能夠容易的獲得,所以國內外的學者大都以這幾個因子來劃定 土石流的臨界降雨線或設立土石流預警之基準。
2.4 水流入滲之效應
水由重力及毛細力而進入土中,因此當降雨直接、間接落於地 面,不論雨量是否足以聚集成地表逕流,皆會影響土壤水份的狀況。 若降雨強度小於入滲率,則地面無逕流發生,雨量多轉變為地下水 流;若降雨強度大於或等於入滲率,地面即生逕流;又降雨初期的入 滲率最大,其值會在數小時內快速降低而達到一平衡值。 雨水會由土壤間的孔隙向下入滲,因為土體之結構不同,會有 不同之入滲量,使得地下水位的變化情形亦不同,本節將土層內之孔 隙水壓、滲透情形及地下水流變化狀況作一分析探討。 前期降雨能使邊坡表層充滿水分,使得水能在邊坡中更容易流11
動。而邊坡破壞所需要的前期降雨量是依據土壤表層的覆蓋、土壤之 水力傳導度、植生的蒸散和邊坡的水文情況來決定的。在非飽和的土 壤中,水分被土壤的吸力或稱負的孔隙水壓力留在土壤孔隙中,而前 期的降雨關係到土壤的水分含量和孔隙中水的張力(Johnson & Sitar, 1990)。根據 Johnson & Sitar (1990),暴雨發生在濕的情況下,比暴雨 發生在乾的情況下,更能產生正的孔隙水壓力。而隨著孔隙水壓力之 上升,將逐漸降低土壤之有效應力,進而導致邊坡發生不穩定之現象。
Fiorillo & Wilson(2004) 以ㄧ厚 1 m 的土壤為例,其孔隙率 n =
0.57,說明田間含水量(field capacity)θmax= 0.51 時累積變化的情形,
如圖2. 3 所示。在體積含水量 θ 上升至 51%前,水分不斷累積,且
保持在毛細現象所維持的孔隙及吸附水中。當雨量累積與蒸散作用達 到平衡時,此時的土壤含水量稱為田間含水量(field capacity)。當降 雨入滲超過排水速率時,正的孔隙水壓便會產生。
Atkinson & Farrar (1985)於調查英國高速公路路堤之淺層破壞, 並埋設水壓計,以便觀察孔隙水壓,其結果發現路堤邊坡滑動係因孔 隙水壓激發造成。Tarantino & Bosco (2000)舉出地滑造成之土石流常 發生於短暫延時之高降雨強度下、延時較長之小雨,甚至發生於降雨 停止後數小時內。會造成此種現象可能原因為降雨入滲造成之濕潤面 影響土層穩定。
12 前期降雨對邊坡穩定的影響已經被研究了很多年,Lumb (1975) 發現前期降雨對邊坡破壞的影響,特別是,他發現如果前期降雨量較 高的話,將會伴隨著較多的邊坡破壞事件。 根據所獲得的降雨資料,Lumb 訂定了不同等級事件的範圍,並 以15 天的前期降雨和日雨量來說明。最嚴重的事件發生在日雨量超 過100 mm,其 15 天的前期降雨量超過 350 mm;而嚴重的事件發生 在日雨量超過100 mm,其前期降雨量達到 200 mm。台灣地區夏季炎 熱且多雨,若遇上颱風或豪大雨時,累積雨量甚至可高達1000 mm 以上,故造成許多災害。
13
第三章
碎形維度理論與相關推導
3.1 碎形理論
碎形的概念可以從兩方面建立起來:首先,畫一個線段、正方形 和立方體,邊長都是 1。將它們的邊長二等分,此時,原圖的線段縮 小爲原來的 1/2,而將原圖等分爲若干個相似的圖形。其線段、正方 形、立方體分別被等分爲21、22和23個相似的子圖形,其中的指數1、 2、3,正好等於與圖形相應的維度。一般說來,如果某圖形是由原圖 以1/a 的比例縮小,且其總數為 b,則: aD=b, (3.1.1) D=logb/loga (3.1.2) 指數D 稱爲相似性維度,亦即碎形維度。如上例 a = 2 ,D = 1、2、 3。其中 D 可以是整數,也可以是分數。 碎形為具有擴展對稱性的幾何對象。擴展對稱性又稱為自我相似 性,它指的是:對一複雜的幾何對象,適當地取出一部份,並加以放 大,觀察者所見之結果與整體對象完全相同,此種具有擴展對稱性的 對象,即使在尺度的變化下亦是不變的。14
3.2 方格覆蓋維度
在盒子計數法中,還有一種應用更為廣泛之測量碎形維數的方 法,並特別適用於計算機的模擬。首先把空間的邊長分割成t 等份, 則產生許多大小相同的小立方體,然後計算覆蓋某一形狀所需小立方 體的數目。如果所研究的碎形位於平面上,則把平面分割成邊長為t 的許多小方格,依照上述方法計算覆蓋某一形狀所需方格的數目。 方格覆蓋維度(box-counting dimension)由於概念與計算並不複 雜,因此常被應用於物理學與地理學之中。就平面上的碎形而言,只 需將碎形圖像放在適當大小的方格中,再計算碎形圖像佔據了多少小 方格即可。 如圖3.1 所示,對同一圓形顆粒樣本分別利用 6 種不同尺規 (scale)的網格覆蓋。再以 6 種尺規為橫座標,以各尺規所對應覆蓋 格數為縱座標,將之繪成雙對數圖形並利用最小平方法回歸得一直 線,該直線斜率的絕對值即為此圓形顆粒之方格覆蓋維度DB。 同理,在篩分析試驗,一系列由大至小變化排列之篩網可比擬作 圖3.1 之 6 種不同尺規方格,所不同的是原本受測顆粒母體經過篩分 析後已被分為多組小母體於各篩網上(如圖3.2(a)、(b)、(c)、(d))。 換言之,各篩網上顆粒之小母體皆不同,而只是原本顆粒母體的一部15 份。故若欲將方格覆蓋維度的計算觀念應用於解釋篩分析上(如圖 3.2 所示),即必須要使用同一組顆粒母體作分析計算,而只能改變篩 網網格孔徑(方格尺規)大小。但是做篩分析試驗時,同一組粒料母 體顆粒必會分別停留於各號篩網上,所以每號篩網上顆粒小母體都應 具有各自不同篩網尺寸之覆蓋格數。 圖3.1 在不同網格尺規下同一顆粒之覆蓋方格數示意圖
16 圖 3.2 方格覆蓋法邏輯與篩分析過程之比對(楊振榮, 2001) 圖 3.3 篩分析試驗各號篩停篩結果示意圖(卓佳良,1999) 1/4 1/8 1/16 (a) 方格覆蓋法 (b) 篩分析
17 以圖3.3 之篩分析結果為例:若令四個孔徑 d1、d2、d3及d4篩網 上顆粒的獨自覆蓋格數分別為
( )
N1 a、( )
N2 b、( )
N3 c及( )
N4 d,則經轉換 並累加相當覆蓋格數後,可得原顆粒母體於篩網尺寸d1網格上之覆 蓋格數應為:( ) ( ) ( ) ( )
1 1 a 1 b 1 c 1 d N = N + N + N + N (3.2.1) 其中,( )
Ni b、( )
Ni c及( )
Ni d等三數(i = 1)分別為三篩上所停留顆粒, 轉換對應於網格尺寸d1時之覆蓋格數。經此種轉換並累加格數,才 符合方格覆蓋維度之計算中保持使用同一組粒料母體的概念。此外, 其他每個篩號之篩網都必需重複這個種轉換並累加格數的動作,如此 同一組粒料才會歷經不同尺寸篩網之檢測以得所對應覆蓋格數。同 理,可對N2、N3及N4作同樣處理方式。這些累加後的總格數 N1、
N2、
N3及N4,即表示同一組粒料母體於各篩號篩網上所覆蓋之總格 數。即原本顆粒母體於篩網尺寸d2上之覆蓋格數應為:( ) ( ) ( ) ( )
2 2 a 2 b 2 c 2 d N = N + N + N + N (3.2.2) 而於篩網尺寸d3上之覆蓋格數應為:( ) ( ) ( ) ( )
3 3 a 3 b 3 c 3 d N = N + N + N + N (3.2.3) 而於篩網尺寸d4上之覆蓋格數應為:( ) ( ) ( ) ( )
4 4 a 4 b 4 c 4 d N = N + N + N + N (3.2.4) 式(3.2.2)~式(3.2.4)中,( )
Ni a、( )
Ni b、( )
Ni c及( )
Ni d(i = 2、3、4)18 意義類似於式(3.2.1)。最後將累加的總格數Ni(i = 1、2、3、4)與 其對應篩網尺寸大小di(i = 1、2、3、4)繪成雙對數圖,並以最小 平方法求得一迴歸直線,則此直線斜率的絕對值即為所求的方格覆蓋 維度DB(如圖3.4 所示)。 圖3.4 方格覆蓋法計算過程示意圖 (卓佳良,1999)
19
3.3 現地與篩分析試驗之碎形維度差異
篩分析試驗是屬於空間中三維度之試驗,因此其所求出之碎形維 度會介於2~3 之間;而現地照片是屬於平面上二維度之量測,其所求 出之碎形維度會介於1~2 之間,理應兩者無從比較。但根據二維與三 維自我相似理論的結論,如果碎形母體具自我相似特性,則二維碎形 維度會與三維碎形維度僅差1。因此若現地顆粒級配具自我相似特 性,即可將現地照片之二維碎形維度加1,以便與篩分析試驗結果之 碎形維度作比較。3.4 土壤保水曲線碎形模式
Arya & Paris (1981)以土粒之粒徑分佈、乾密度、統體密度等已 知參數,由已知土粒粒徑分佈轉換成孔隙大小,提出一保水曲線模 型:由土粒體積與孔隙體積關係,可知第i 級土粒粒子所對應之孔隙 體積VV(i),如下: ( ) ( ) ( ) p i V i i p m V e ρ = (3.4.1) 其中,mp(i) = 第 i 級土粒 p 的質量、ρp = 土粒 p 的密度、e(i ) = 第 i 級 孔隙比。若設ρb = 土壤統體密度,則土體總體積為 Vb ,土粒體積為
20 Vp。 ( ) ( ) ( ) ( ) / / / b p i b p i p i p p i p V V m m e V m ρ ρ ρ − − = = (3.4.2)
若孔隙內不含水時,m(i) = mP(i),則可進一步導得此時 e(i) = (ρp-ρb)
/ ρb。由於土粒的密度 ρp可假設為一定值,若設土體統體密度ρb固定, 則各級孔隙比e(i)亦等於整體孔隙比 e。 再將土粒粒徑大小分佈分成i 個等級、設每等級的顆粒數有 N(i) 個,且平均半徑為R(i),則顆粒體積為 ) ( 3 ) ( ) ( 3 4 i i i p R N V =
π
(3.4.3) 進一步假設每個等級顆粒間之孔隙,僅形成一條毛細管通過,且毛細管半徑為r(i)、長度為h(i),則其孔隙體積VV(i)可計算為:
) ( 2 ) ( ) (i i i v
r
h
V
=
π
×
(3.4.4) 其中,若在每個等級內土粒球形粒徑(R(i))大小相等,則毛細高度h(i)等於2 × R(i) × N(i);但一般土粒常不為球形,則毛細管真實高度或
長度可修正為: α ) ( ) ( * ) (i
2
R
iN
ih
=
(3.4.5) 式中,修正因子α 經驗值介於 1 至 2。該式意義即土粒粒徑減小則孔 隙管徑急速減小,毛細高度將以指數倍數升高。 再將式(3.4.5)中 h(i)* 代入式(3.4.4)中可得:21 2 ( ) ( )
2
( ) ( ) v i i i iV
=
π
r
i
R N
α (3.4.6a) 又因VV(i) = e Vp(i),故由式(3.4.3)可得 3 ( ) ( ) ( ) 4 3 v i i i V = ieπ
R N (3.4.6b) 令式(3.4.6a)=(3.4.6b)可整理得毛細管徑為 2 1 1 ) ( ) ( ) ( 3 2 ⎥⎦ ⎤ ⎢⎣ ⎡ = i i−α i N e R r (3.4.7) 而毛細水頭高度公式為: ( ) ( ) 2 cos i i T h rβ
= i (3.4.8) 其中,r(i)為孔隙之半徑,T 為表面張力,β為接觸角,故將(3.4.7) 代入式(3.4.8)可得各等級土粒間之毛細水頭 h(i): 1 2 1 ( ) ( ) ( ) 2 cos 2 3 i i i T e h N R αβ
− − ⎡ ⎤ = ⎢ ⎥ ⎣ ⎦ i (3.4.9) 式(3.4.9)中 α 的意義,直到 Tyler & Wheatcraft(1989)引進碎形理 論觀點(如圖 3.5 所示),認為 α 可能是管道不規則度(degree of irregularity),並證實可用描述粒徑分佈之碎形維度 D 替代之,即 α = D,如圖 3.4 所示。其中 D 值之推求可將土體分成數等級,再以土壤 重量與粒徑推估土壤顆粒數,再依碎形理論求知碎形維度D(其中 1≦D<2)。對土壤而言,當 D = 1.5 時,管道路徑已屬非常曲折;而 D 越大,他們表示含有越多細粒黏土質土壤,去水速度較慢;例如粗22 粒土壤孔隙之D = 1.46,細粒土壤孔隙之 D = 1.89。 因此以 α = D 代入毛細水頭公式(3.4.9)得: 1 2 1 ( ) ( ) ( ) 2 cos 2 3 D i i i T h e N R
β
− − ⎡ ⎤ = ⎢ ⎥ ⎣ ⎦ i (3.4.10) 由於水對較小孔隙之吸附能力較強,因此通常是由較小孔隙填 滿後再填較大孔隙,故其含水量亦可分為n+1 級,而第 i 級之體積含水量(volumetric water content)便是水由第 i 級之孔隙填滿至第 n+1 級所需之體積除以土體之總體積而來,其可表為: 1 ( ) ( ) n V i i i b V V ω + =
∑
(3.4.11) 由式(3.4.10)及式(3.4.11)即可求得毛細水頭。其中,由土 體組構關係可知(Arya & Paris ,1981):3 p 1 v( i ) b( i ) ( pi ) ( i ) ( i ) b V V V N R ρ ρ ⎛ ⎞ = − = ⎜ − ⎟ ⎝ ⎠ (3.4.12) 設 e = e(i)由式(3.4.3)、式(3.4.12)及原始孔隙比(e)定義知: ( ) 3 3 1 3 1 4 4 3 p i i p b vi i pi b i i N R V e e V R N
ρ
ρ
ρ
π ρ
π
⎛ ⎞ − ⎜ ⎟ ⎛ ⎞ ⎝ ⎠ = = = = ⎜ − ⎟ ⎝ ⎠ (3.4.13) 故僅須知道土粒粒徑分佈特性D、土粒密度 ρp、土體密度 ρb等參數, 便可預測保水曲線。23
3.5 案例驗證
Turcotte(1986)曾經證實土粒的分裂是一種碎形特性,土粒粒 徑愈小,其數量愈多,即 N(i)×R(i)D為一個定值。故各等級土粒間具有 下列關係: ( ) ( )* ( 1) ( 1)* D D i i i i N i R = N + i R + (3.5.1) 其中,D*為土壤顆粒粒徑分佈的碎形維度;土粒的 D*值常接近 於3;N(i) = 粒徑大於 R(i)之顆粒總數(累積停篩顆粒總數)。若兩邊 取對數可得: ⎟⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + + ( 1) ) ( * ) 1 ( ) ( log log i i i i R R D N N (3.5.2a) 即−
D
*=
∆
∆
(log N )
(log R )
(3.5.2b) 類比方格覆蓋法,可由傳統進行篩分析所得到的粒徑分佈曲線 半對數關係:將縱軸由(% finer)改為 log (% finer)重新以雙對數關 係表示粒徑分佈曲線,便可由圖形中之斜率得知該顆粒大小分布之碎形維度,如圖3.6 所示。
Tyler & Wheatcraft(1989)指出碎形增額 Di(即碎形維數與拓
24 度空間顆粒群體,故拓璞維數DT = 3,所以土粒的 Di: Di = D* −3 (3.5.3a) 孔隙管道是通過此一土粒骨架之通道,故具相同的碎形增額亦 是Di,而且孔隙管道之DT = 1,所以孔隙管道之碎形維數 D: Di = D−1 (3.5.3b) 由式(3.5.3a) = 式(3.5.3b)得孔隙管道之碎形維數 D,與土粒粒徑分佈 的碎形維度D*,具下列關係: D = D* −2 (3.5.4) 即土粒粒徑分布之碎形維度,與描述管道不規則度之碎形維數 值差2。因此,藉由篩分析試驗資料,可先得知土粒粒徑分布之碎形 維度D*,再以式(3.5.4)關係,獲知土壤內部毛細管孔隙管道分佈之特 性D 值,其值反應毛細管道彎曲的程度 D 值愈大,毛細管道愈不規 則,其意義與前述α 一致。 本研究摘錄Bousnina(1984)在距地表 0-0.2 m 處所採得含黏 土質砂土之現地保水試驗資料,以驗證本文模式之正確性。其土粒粒 徑分佈資料如圖3.7 之篩分析粒徑分佈曲線;現地保水試驗資料則繪 於圖3.9 之保水曲線。利用土壤保水曲線碎形模式便可得到繪製保水 曲線之參數,如表3.1 所示。其中對碎形維度 D 值之求得,本研究採 用兩種方法,分述如下。
25 (1)改良類比方格覆蓋法 本研究中由於考慮各孔徑篩上之土體重量與篩孔直徑之關 係,故將原類比方格覆蓋法中之縱軸改以累積停篩百分比(log(% accumulation mass)),橫軸仍維持原狀作圖,重新將粒徑分佈曲線以 雙對數關係表示,結果如圖3.7 所示,可由圖形中之斜率得知顆粒大 小分布之碎形維度 (Dr = 3 - 斜率,Dr 為改良類比方格覆蓋法所 得之碎形維度値)。以此法將 Rieu & Sposito(1991)土壤粒徑與累積
停篩重量百分比作圖得斜率為0.61,依前述理論得知描述土粒的碎形 維度3-Dr=0.61,即 Dr=2.39;故表達分佈於土粒間之毛細管道彎曲 程度碎形值為D=0.61,可直接以原累積停篩粒徑分佈曲線雙對數圖 迴歸直線之斜率代替之,所以毛細管道的碎形維度D 即為等於該圖 的斜率。 (2)管道不規則度估算法 於篩分析中,可由兩篩篩網網徑之平均數當作顆粒粒徑,再將 粒徑除以2 得其各部分之 R(i)。由停留在此篩徑上之土重,除以土粒 之單位重得其體積,再除以顆粒體積((4/3)×π×R(i)3)而得各篩上之 土粒數目N(i),將各篩之顆粒數加起來便得到N 值。將篩孔直徑-土 壤顆粒數取雙對數作圖即可得一斜率,如圖3.8 所示。由圖形中迴歸 直線之斜率得知該顆粒大小分布關係之碎形維度D*=2.66,利用 D*
26 值可以反應孔隙大小分佈性質或毛細管道彎曲的程度。以兩種方法得 到Dr=2.39 及 D*=2.66,再由公式(3.4.10)換算得毛細管道之 D=0.66。 將上述二法所得之D 值分別代入,得到兩種保水曲線預測值, 並與現場保水試驗作比較如圖3.9 所示。由圖 3.9 的結果中,發現此 二預測曲線趨勢與原始現場實驗所得之保水曲線趨勢相符。且以改良 類比方格覆蓋法所得斜率做為毛細管道D 值預估時,曲線前半部原 實驗值有低估的情形,在曲線後半部與原實驗值相較有高估的情形; 而以管道孔隙不規則度估算法所換算得之D*值所得結果趨勢與真實 情形較為類似。故由比較得知,若將保水曲線分為二個部分,前面部 份以由管道孔隙不規則度估算法所得D 值推估與現場實驗所得之保 水曲線趨勢較為相符,後面部份D 值以改良類比覆蓋法推估與現場 實驗所得之保水曲線趨勢較為相符。但整體而言,孔隙不規則度估算 法所得知結果較佳,故本研究擬採用此一模式,以預測保水曲線。
27 圖 3.6 雙對數關係表示粒徑分佈曲線示意圖(楊振榮, 2001) 0.1 1 10 10 100 mean size (mm) A ccu mu la ti on m as s (% ) 圖3.7 篩分析曲線雙對數圖(Slope=0.61)
28
圖3.8 改良篩分析曲線對數圖(Slope=2.66)
29 表 3.1 估算保水曲線之基本參數 累積停篩 百分比 (%) 顆粒數 (E+15) 各層孔隙 體積Vv(i) 各層含水量 w(i) 水頭高Hi (D* = 2.66) (cm) 孔隙比e 水頭高Hi (D* = 2.39) (cm) 13.864 0.034 21154 0.296 2.202 0.203 1.01 34.735 0.105 31845 0.446 2.306 0.203 1.029 47.168 0.127 18971 0.266 2.825 0.203 1.254 57.503 0.210 15771 0.221 3.26 0.203 1.429 66.846 0.373 14254 0.2 3.705 0.203 1.601 73.176 0.513 9660 0.135 4.442 0.203 1.905 82.121 2.001 13648 0.191 4.944 0.203 2.049 89.276 6.48 10918 0.153 6.452 0.203 2.597 96.899 38.94 11631 0.163 8.468 0.203 3.258
3.6 滲透曲線碎形模式
3.6.1 理論推導 根據碎形維度理論,土壤孔隙分布數目可寫為下式: ( ) D N r =Cr− (3.6.1) 其中N 為在半徑 r 範圍內的孔隙數目,C 是常數,D 為孔隙之碎形維 度值。因此,土壤孔隙體積分佈關係式可寫為: 3 ( ) D f r =Ar − (3.6.2)30
其中A = 4πC/3。將半徑 r → r+dr 範圍之孔隙填滿所需含水量可寫為:
( )
dΛ =df r (3.6.3)
式中Λ 為相關體積含水率(relative volumetric water content),Λ = θ-θr,
θ 與 θr為體積及殘餘含水率,再將(3.6.3)式積分為下式: 3 (3 ) D A r D − Λ = − (3.6.3a) 而相關體積含水率在飽和時可將(3.6.3a)式改寫為下式: 3 (3 ) D S A R D − Λ = − (3.6.3b) 式中ΛS為飽和時之相關體積含水率,R 為孔隙之最大半徑值。 此時水之表面曲率為 r,因為在相關體積含水率 Λ 時,孔隙半徑 小於r 者皆被水填滿。而水力勢能 ψ 與孔隙半徑之關係可改寫自(3.4.8) 式 2 cosT r β ψ = (3.6.4a) 2 cos e T R β ψ = (3.6.4b) 式中ψ 為水力勢能,ψe為空氣進氣値(air-entry value),T 為表面張力, β 為接觸角。將(3.6.4a)、(3.6.4b)式與(3.6.3a)及(3.6.3b)式結合便是體 積含水率與水力勢能之關係式, e δ ψ ψ ⎛ ⎞ Θ = ⎜ ⎟ ⎝ ⎠ (3.6.5) 其中Θ 為正規化之體積含水率値( Θ = Λ/ΛS ),δ 是一常數,其與碎形 維度之關係為:
31 3 D δ = − (3.6.6) 對照(3.5.3a)式此處 D 值即為土粒粒徑分布之碎形維度 D*值。 Burdine(1953)發展出一相關水力滲透方程式(水力滲透係數與飽 和度之比例式) 2 0 2 1 2 0 r d k d ψ ψ Λ Λ = Θ
∫
Λ∫
(3.6.7) 其中kr為相關水力滲透係數,kr = ks/kw; kw與 ks為在任何飽和度下之 滲透係數與飽和狀態下之滲透係數。 將(3.6.5)式代入(3.6.7)式,並且考慮 dΘ = dΛ,則相關水力滲透方 程式可得到下式: r k = Θλ (3.6.8) r e k η ψ ψ ⎛ ⎞ = ⎜ ⎟ ⎝ ⎠ (3.6.9) 3 11 3 D D λ = − − (3.6.10) 3D 11 η = − (3.6.11) (3.6.8)式顯示了相關水力滲透係數與正規化之體積含水率之關係,而 (3.6.9)式顯示了相關水力滲透係數與水力勢能之關係。 由土壤顆粒之孔隙分布進而推導出非飽和土壤之水力滲透係數 之重要參數可由(3.6.10)式及(3.6.11)式決定。32 3.6.2 滲透曲線碎形模式探討 本研究以 Bousnina(1984)在距地表 0-0.2 m 處所採得含黏土質 砂土之現地資料,討論在土壤水分分布情形於不同顆粒排列下,其相 關水力滲透係數與正規化之體積含水率及水力勢能比之關係。此處之 原始水分分布情形如圖3.9 及表 3.1 所示。其空氣進氣値為 2.0 kPa, 在假設土粒分布分別維度值為D = 2.2、2.5、2.7 三種情形下之水分變 化情形。 圖 3.10 表示相關水力滲透係數(kr = ks/kw; kw與 ks為在任何飽和 度下之滲透係數與飽和狀態下之滲透係數)與水力勢能値之關係圖。 如圖所示,隨著水力勢能値增加,其相關水力滲透係數隨之減少,其 變化幅度在維度值較低時(D=2.2)較為明顯。但在水力勢能値較小 時,相關水力滲透係數幾乎重合,這表示在低水力勢能値時,土壤顆 粒的變化並不會對相關水力滲透係數造成太大的影響。 圖 3.11 表示相關水力滲透係數(kr = ks/kw)與正規化之體積含水 率値(Θ = Λ/ΛS,Λ 及 ΛS分別為任何飽和度下之相關體積含水率與在 飽和狀況下之相關體積含水率)關係圖。由關係圖得知在高相關體積 含水率情形時,可得到高相關水力滲透係數值;此外,其顆粒變化時 (D = 2.2、2.5、2.7 三種情形下)相關水力滲透係數變化明顯,尤
33 其在維度值較高時(D = 2.7)最為明顯。此圖表示體積含水率與滲透 係數之關係受到顆粒分布影響十分顯著,完全沒有重合的情形發生。 0.1 1 10 1E-007 1E-006 1E-005 0.0001 0.001 0.01 0.1 1 D=2.2 D=2.5 D=2.7
Water potential (kPa)
R e la ti v e w a te r p e rm e a b ili ty 圖3.10 不同顆粒分布下相關水力滲透係數與水力勢能之關係圖 0.1 1 1E-008 1E-007 1E-006 1E-005 0.0001 0.001 0.01 0.1 D=2.2 D=2.5 D=2.7
Normalize volumetric water content
R el a ti v e w a te r p e rm e a b ili ty 圖3.11 不同顆粒分布相關水力滲透係數與正規化體積含水率關係圖
34
第四章
影像於篩分析之應用
4.1 二維影像分析之研究
大部分的研究中,多以實驗室所得篩分析曲線描述粒徑分佈,僅 有少部分採用現地照片,進行土石流材料完整顆粒大小之判識。 (Bonnet-Staub, 1999, Ter-Stepanian, 2002, Varnes, 1978, Znamensky& Gramani, 2000)。 Mancini et al.(1987)提出影像分析之步驟如下: (1)在土石流發生區內取樣,且影像取景必須垂直於取樣平面,同 時照片中需含比例尺 。 (2)無論使用影像軟體或人工判讀照片,必須同時紀錄其顆粒之周 長、長軸長度、短軸長度、等值圓半徑、及顆粒距離影像圖形 之中心軸位置及密度等資料。 (3)顆粒大小的量度,必須以圖中比例尺為衡量標準。 (4)為避免大顆粒遮蔽效應,經驗上以分佈係數作為計算之比例。 (5)根據上述資料所得之體積百分比,需經由經驗修正才可轉換至 重量百分比。
Baroni et al.(2003)於義大利 Ravaneti 發生土石流之上游區域實
35 上述Mancini et al.(1987)提出之影像分析方法計算粒徑分佈曲線, 結果發現粒徑在10mm~100mm 範圍內,以傳統篩分析及影像分析所 得粒徑分佈曲線斜率相同。Baroni et al.認為可將兩種結果疊合,成為 完整粒徑分佈曲線,如圖4.3 所示。 圖 4.1 義大利 Ravaneto 區域影像拍攝位置之斷面全景 (Baroni et al.,2003) 圖4.2 影像分析照片與再製圖(Baroni et al.,2003)
36 圖4.3 影像分析與篩分析比對(Baroni et al.,2003)
4.2 數位影像理論
本研究擬以土石材料現地之數位影像,發展數位影像篩分析法, 此套簡易方式主要功用是繪製無法採樣部分之粒徑分佈曲線。可採樣 部分仍採用傳統篩分析法,而顆粒大小超過標準篩號尺寸之部分則以 二維影像分析取代篩分析。即嘗試以兩種方法之結合彌補傳統篩分析 之不足,以繪製完整粒徑分佈曲線。本研究擬進行二維影像分析,所 謂二維影像分析是指以數位相片紀錄現地土石之分佈情況。而影像分 析前,必須先了解數位影像的一些基本特性。 4.2.1 影像分割 影像分割(image segmentation)的主要目的是從影像中,將有興趣37 的目標、物體、物件分割出來,再分析所得的有效資訊,經由資料數 據的判別,將不同意義給予不同的指令,並做不同的處理動作,此為 影像分割的功用。 灰階影像是一個二維亮度函數f(x,y),而彩色影像是三個二維亮 度函數R、G、B 組成。茲定義點的空間坐標為 (x,y),該點影像的亮 度或灰度為f(x,y)。影像分割演算法通常以影像亮度值的兩個基本特 性之一為基礎:相似性和不連續性。在相似特性中,因為原屬同一個 物件,其影像亮度值相似,所以根據影像亮度值的相似性,從影像中 找出物件。而在不連續特性中,利用位於物件邊緣的影像亮度值產生 明顯的變化(即影像亮度值的不連續性),找出物件邊緣的封閉曲線, 進而找出物件。一般來說,一個物件可以使用封閉的曲線來表示,同 理,一個封閉的曲線可以描述一個物件,所以無論是找出物件的封閉 邊緣曲線或是直接找到物件,都是達到影像分割的目的。 影像分割主要問題之一是影像資料含糊不清的特性。由於影像常 受到外來雜訊的影響,使得影像呈現出亮度不均,或物件與背景間亮 度值產生低對比的情況。影像中的物件往往受到光影上的影響,造成 原本應該一致的顏色可能變為漸層色,或是有些地方變為陰影,有些 地方則變為反光面,造成這些區域原本應該被分類為同一物件,卻因 為亮度不均而被分類為另一物件。
38 對於影像分割來說,除了使用影像亮度值的基本特性外,還需要 一些額外的資訊來輔助,以便達到分割物件的目的。影像分割的方法, 根據使用的特性分析可以分成四大類: (1)門檻值法(Thresholding):利用統計整張影像或部份影像的灰階統 計長條圖 (histogram) 來分析。 (2)以邊緣為基礎的分割法(Edge-based Segmentation):利用物件和背 景間,影像亮度值產生明顯變化的特性找出物件邊緣。 (3) 以區域為基礎的分割法(Region-based Segmentation):利用同一物 件中,其影像亮度值相似的特性,從影像中找出物件。 (4)分水嶺分割法(Watershed Segmentation):使用前面三種方法混合的 特性分割影像。 在這一節中,將分別描述本研究使用的二種方法:分水嶺分割法 與門檻值法中的二值化法。 4.2.2 分水嶺分割法 (1)形態學梯度(morphological gradient) 在邊緣檢測影像處理中有多種梯度計算方式,其基本原理大多基 於下面考量:如果影像中某一像素點位置的梯度值大,則表示在該點 位置影像有快速明暗變化,意味該點可能為邊緣通過。這些梯度計算
39
多半以數學微分形式求得(Michael , 1999)。在數學形態學影像處理 中也提出了幾種梯度計算方法,此處僅討論本研究所用的形態學梯 度。形態學梯度運算能夠針對影像中物件邊緣創造出一邊緣加強的影
像,此影像的產生過程如下:複製原始影像A 成為 A’,對 A 進行侵
蝕,對A’進行膨脹。接著從膨脹運算後的 A’減去侵蝕運算的 A。由
於被膨脹過影像物件的大小會輕微地增加,而被侵蝕過影像物件的大 小會輕微地減少,兩者之間的差異便能標示出物件的邊緣(Baxes , 1994)。圖 4.4(a)~(d) 展示計算一灰階影像的形態學梯度過程。(a)為 輸入影像,(b)為膨脹後影像,圖(c)為侵蝕後影像,計算圖(b)與圖(c) 的差異即產生梯度影像圖(d)。 (a)輸入影像 (b)膨脹後影像 (c)為侵蝕後影像 (d)梯度影像 圖4.4 一灰階影像計算形態學梯度影像過程 影像將只會呈現出物件的輪廓並使得輸入影像中的尖銳灰階變 化過程變亮了。形態梯度運算在一維灰階處理過程展示則如圖4.5 所 示。
40 圖4.5 形態學梯度運算處理過程示意圖(Baxes , 1994) (2)分水嶺演算法 (Watershed algorithm) 分水嶺演算是一種流域方法概念。將流域方法應用在待分割影像 的形態學梯度,可以得到非常有效的灰階值分割方法。假設待分割的 影像由不連接的物體組成,且物體具有較暗的灰階值。影像中具有均 勻低灰階值得區域將其稱為極小區域。影像中有三種空間點:(a)屬 於極小區域的點;(b)若一滴水珠落於影像上某點處,水珠必定滾入
41 極小區域的點;(c)水珠由該點滾入一個以上極小區域的可能性相同 之點。對於一個給定的極小區域,水珠會滾入該區域的所有點所構成 的集合,稱為該極小區域的流域或集水域。水珠從脊線滾入一個以上 極小區域的可能性均等的點所構成的集合,稱為分水嶺線(Baxes , 1994)。此概念可參見圖 4.6。 圖4.6 流域、極小區域、分水嶺線概念圖(Baxes , 1994) 依上所述,影像分割問題可歸類成一種求影像流域分界線的問 題。現在假設水珠不是從流域滾下,而是從水溢出的角度來看,假設 影像中極小區域部分被打了小孔,水從這些小孔中以一定速度溢出, 當不同流域中的水面不斷升高到要匯合在一起時,築起一道堤防防止 不同流域中的水匯合。水面持續上漲直到所有的分水嶺線被找出
42 (Baxes , 1994)。此概念可參見圖 4.7。 圖4.7 流域築堤概念圖(Baxes , 1994) 從地形學(topographic)的角度來看,把一張影像的灰階值看成是 地形上高低起伏的表面高度時,運用地形學的概念,影像就猶如有高 低起伏的地理景觀。因此,根據這些灰階值,將每個局部區域高度最 小值當成是地表一個盆地(basin)的最低點,然後想像用水逐步浸入, 模擬水從不同盆地的最低點開始慢慢的上漲,直到來自兩個不同盆地 區域的水將要混合時,便建立一分隔水壩將彼此隔離,而此分隔水壩 就等於所要找的分水嶺線 (watershed line)。直到水漲到最高點後,所 有的分水嶺線就都能被找出來,而每一個分水嶺區間被給予一個單一
43 編號,此時就完成了分水嶺演算法分割過程。這樣的方法稱為是以浸 泡(immersion)的模擬方式求得分水嶺演算法,此方法快速、簡單、準 確且計算效率高,本研究所採用的分水嶺演算法即以此為基礎。一個 浸泡式分水嶺演算一維灰階處理過程展示見圖4.7。 在整個處理過程中,通常是以水漲的高度來決定影像中的像素是 否要被處理,也就是以像素的灰階值大小來與水上漲的高度來做比 較,當水漲高度等於某一灰階值H 時,所有小於 H 的灰階值是已被 處理完畢並給一適當標籤(label)值。此方法是以每個像素為中心,依 其周圍8 個鄰近點的標籤值來決定此像素的標籤值(Vincent , 1991)。 4.2.3 門檻值法(Thresholding) (1)二值化(Binarization) 影像二值化處理作業是將每個像素的灰階值,降低至特定範圍的 亮度值。每個二值化的作業,將把整張影像轉變成兩個灰階值:0, 以及影像中最大的灰階值(若是 8 位元的影像,則最大值即為 255)。 在進行二值化作業時,門檻值的選取是相當重要的。若門檻值設 定不理想,經過二值化的影像會遺失許多原有的細節,將無法找出待 測物體的正確邊界與某些特定輪廓特徵。
44 (2)二值化基本概念 二值化的運算方法,是將整張影像中,每個像素之灰階值與預設 之門檻值(Threshold)逐一比較,若高於門檻值,則該像素的灰階值 將被設定為最大值;反之,則設定為0。 以圖 4.8 為例,假設(a)中所示的灰階直方圖對應到由陰暗背景上 的明亮物體所組成一幅影像f(x , y)上,使得物體和背景的像素被分類 為兩個主要群體的灰階。若要從背景中抽取出物體的一個明顯方式即 選擇分開這些群體的門檻值T。於是,對 f(x , y)>T 的任一點(x , y)稱
為物體點(Object Point),其餘則為背景點(Background Point)。
圖 4.8(b)為一般較常見的情況,其中在此影像直方圖上有三個主
要群體的特徵(例如,在陰暗背景上的兩種明亮物體)。此處,如果
T1(x , y)≤T2,則多階門檻值法(Multilevel Thresholding)將點(x , y)歸
屬於一個物體點;如果 f(x , y)>T2,則歸屬於另一個物體點;如果f(x ,
45 (a)單一門檻值 (b)多門檻值 圖4.8 門檻值法之灰階直方圖 根據以上所述,門檻值法可以看成是一種依據函數T 檢測的運 算,其中T 的形式為:
( ) ( )
, , , , , = ⎡⎣ ⎤⎦ T T x y p x y f x y (4.2.1)其中 f(x , y)是點(x , y)的灰階,p(x , y)代表了這一點的區域性質—
例如,以(x , y)為中心其鄰域的平均灰階值。一個門檻值化的影像 g(x , y)定義為:
( )
, 1 ( , ) 0 ( , ) f x y T g x y f x y T > ⎧ = ⎨ ≤ ⎩ 若 若 (4.2.2) 因此,這些標記為1(或其他自訂的灰階值)的像素對應於物體, 標記為0(或其他沒有指定給物體的灰階值)的像素對應到背景。 當 T 只跟 f(x , y)有關時,這個門檻值稱為整體(global)門檻值。46 如果T 與 f(x , y)和 P(x , y)有關,這個門檻值稱為區域(local)門檻 值。 (3)Otsu Method Otsu(1979)提出運用組間(between-class)差距越大越好的概念來 進行二值化分割。若灰階影像之長條圖具有雙峰特性,則可以利用統 計學的原理,來找出最佳的門檻值,以分割灰階影像中的二個群集的 像素。設一已知影像像素的灰階分佈為[ 1 , 2 , … L ],灰度值為 i 的 像素個數ni,且像素之總和N=n1+n2+n3+…+ni。則機率分佈為:[ 1 , 2 , …, K ] 之 C0像素群;C1則表示灰度值範圍為[K+1, …L]之像素群。 此時二個群集之機率分佈分別為: 0 r 0 1 1 r 1 1 P ( ) ( ) P ( ) 1 ( ) k i i L i i k C P k C P k ω ω ω ω = = + = = = = = = −
∑
∑
(4.2.3) 和平均數( ) ( )
( )
( )
0 r 0 0 1 1 1 r 1 1 1 1 P ( ) P / / P ( ) P / [ ]/[1 ] k k i i i L L i r i K i k i i C i k k i i C i k k µ ω µ ω µ ω µ µ ω = = = + = + = = = = = = − −∑
∑
∑
∑
(4.2.4) 其中 ω = =∑
1 ( ) k i i k P ,µ = =∑
1 ( ) k i i k iP 分別為長條圖之中前K 個值的零階和 一階動量(Moment),且 1 ( ) L r i i L iP µ µ = = =∑
(4.2.5) 為原始影像的總平均數,我們可以任一值來得到下列之關係:47 0 0 1 1 0 1 2 2 2 1 1 1 1 1 1 1 ; 1 ( ) ( ) ( ) / r L L r i i k i k i P i C i P
ω µ ω µ
µ ω ω
σ
µ
µ
ω
= + = + + = + = =∑
− =∑
− (4.2.6) 為了求取更佳的門檻值(第k 個灰度值),可以分辨分析中之分辨標準量測(Discriminant Criterion Measures)來討論。其中
η = σ2 σ2 / B r (4.2.7a) 2 2 2 0 0 1 1 2 0 1 1 0 ( ) ( ) ( ) B r r σ ω µ µ ω µ µ ω ω µ µ = − + − = − (4.2.7b) σ µ = =
∑
− 2 2 1 ( ) L r r i i i P (4.2.7c) 分別是組間(between-class)變異數和總變異數。由於 Otsu 法的概念 是組間差距越大越好(組間變異數越大越好),即 η 越大越好。因此 求η 值是找出閥值 K 最簡單途徑。亦即 η σ σ σ µ ω µ ω ω = = − − 2 2 2 2 ( ) ( ) / ( ) [ ( ) ( )] /[ ( )(1 ( ))] B r B r k k k k k k k (4.2.8) 此方法中,若兩個高峰的高度差異不大,則 Otsu 法可以找到不 錯的閥值。但是當灰階長條圖中,兩個高峰的高度大小差異過大,此 時利用Otsu 法找到的值將會過度偏向總個數比較大的高峰,如此便 無法取得理想的閥值。本研究中大顆粒材料灰階差異不大,因此可以 Otsu 此法獲得不錯之結果。48
4.3 二維影像分析模型試驗
本研究以影像處理分析程序,改良人工描繪之缺點,最終目標冀 望以程式自動搜尋土石邊緣,以利後續估算土石之面積、體積、重量, 並將資料轉換至篩分析曲線。 以一般數位相機拍攝之影像,先將影像轉灰階處理,如圖4.9 所 示。再者,以Otsu 法直接對影像進行二值化處理,結果如圖 4.10 所 示。比較圖4.9 與 4.10,雖然二值化過程明顯判別出物體與背景值差 異,但許多物體均呈相連狀態,對於後續分析十分不利。因為此過程 會低估物體數量及高估物體尺寸,甚至也未達完全分離(segmentation) 的效果。所以在二值化處理前必需先配合其他影像處理方法。 本試驗先以分水嶺演算法將灰階圖像進行分析,結果如圖4.11 所示。將所得結果利用Otsu 法對影像進行二值化處理,結果如圖 4.12 所示。即可獲得邊界較清楚及個別顆粒較明顯之二值化圖形。將此二值化圖形以軟體IPTK(Image Processing Tool Kit, V5.0)進行後處理分
析。IPTK 先將個別顆粒判識、圈繪並標註號碼,如圖 4.13。利用 IPTK
強大之後處理運算技術,將顆粒進行周長、長短軸長、圓度 (Roundness)、中心點位置…等分析,如圖 4.14。
49
圖 4.9 彩色影像轉灰階影像
50
圖4.11 分水嶺演算法後之結果
51
圖 4.13 標註個別顆粒物體
52 經過校正單位後可獲得一些基本資料(如周長、長短軸長、圓度、 中心點位置…等),利用以下幾種體積之估算模式計算體積: (1) Frey (2003)法 2 6 V =
π α
d D (4.3.1)其中,d為中值尺寸(Medium dimensions);D 為最大尺寸(Maximum
dimensions);α 為形狀係數(Shape coefficient)。 (2) Illerström (1998)法 圖4.15 體積幾何示意圖(Illerström,1998) V = ⋅ ⋅ = ⋅ (4.3.2) d d D A d 其中,d為中值尺寸;D 為最大尺寸;d 為最小尺寸;A 為面積。 (3)以球體模擬法 圓球體: 4 3 3 V =
π
r (4.3.3) 其中,r 為等面積圓之半徑;53 橢圓體: 4 2 3 V =
π
ab (4.3.4) 其中,a, b 代表橢圓體最大或最小軸長。若 a > b 為橄欖球狀,a < b 為鐵餅狀。 將經由IPTK 統計資料,分別以上述各種方法計算,並與物理 量測(阿基米得法)所得結果比較,整理如下表所示: 表4.1 各種體積估算結果比較 Frey 法 Illerström 法 圓球體法 橢圓體法 修正圓球 體法 物理量測法 體積 (cm3) 7824.7 6981.7 5676.3 42135.1 4164.9 4080 誤差 (%) 91.8 71.1 39.1 932.7 2.0 0 由表 4.1 之前四欄得知,圓球體模擬結果較佳,因此將圓球體另 乘ㄧ修正係數,結果更接近物理量測 。本研究選定之修正係數為 Compactness factor(緊密係數),所謂 Compactness factor 為量測物體圓 度的指標。計算方法為:Compactness factor=等面積圓直徑/長軸長度。係數介於0 之 1 間,越接近 1 代表物體形狀越接近圓形。反之,
54 利用上述方法進行其他大小石頭影像量測,誤差範圍約介於 2~35.3%之間。推論誤差可能原因來自於(1)石頭排列方式:若遮蔽 越多,則誤差越大(2)前景與背景色差:色差越大,較易將石頭獨 立分離出來;因此石頭與土壤顏色差別越大,分析結果越好。(3)形 狀規則度:形狀較規則、較無稜角之石頭,以公式估算體積較接近真 實大小。
55
第五章
現地採樣及粒徑影像分析
5.1 研究區域地理位置
本研究區域選定位於北部橫貫公路(省道台七線)復興至巴陵段(參 考圖 5.1 ),二度 TM 座標左上角為 280000,2746000;右下角為 2730000,292000。本區域屬於桃園縣石門水庫集水區內,以大漢溪 流域為區域中心。由於自民國53 年起北橫公路通車多年,雖縮短桃 園、宜蘭縣間的距離,但沿途崩坍及土石流常造成交通中斷的情形, 故本研究將以此路段之材料進行物理模型試驗。 北橫公路由復興(17K)至巴陵橋(48.2K),全長約 31.2 公里, 主要沿大漢溪左岸修築,為區域內主要對外交通通道。區域內交通除 了北橫公路外,僅在羅浮附近有桃-188 號公路通往新竹縣關西鎮, 此外也僅有少數產業道路連接聚落與北橫公路。56
57
5.2 研究地點地形概況
從巴陵以下到石門之間為大漢溪的中游,是水庫、攔砂壩等水利 工程集中的河段。在此,大漢溪主流斜切地層走向,向北流經蘇樂、 高義、 榮華、高坡、義興、合流, 並且在合流附近以直角轉向西流, 直抵石門,參照圖5.2。在蘇樂到高坡之間,大漢溪橫切過插天山背 斜構造,而在背斜位置上露出來的岩層,正好是經過輕度變質的堅硬 板岩層。岩性堅硬,維持陡坡而不易崩壞,使得板岩為峽谷的形成提 供了有利的條件。這段河谷深窄,兩岸裸露岩石絕壁,是標準的峽谷 地形,被稱為高坡峽谷。特別是匹亞外到榮華攔砂壩之間的一段,大 漢溪橫切插天山背斜軸部,使得河谷更具峽谷外形。高義與高坡之間 常有堅硬的砂岩構成的山脊突出,形成顯著的交切山腳,位於北橫公 30.5~31 公里、41~45 公里間者最為顯著。至於峽谷的南、北兩端, 巴陵到蘇樂之間以及高坡到羅浮之間,露出來的地層是由砂、頁岩互 層所構成,使得大漢溪在這兩段的河谷,擁有較不穩定的邊坡,河谷 因而比較開闊,並且出現河階地形,例如巴陵、高坡、義興、合流等 地。合流又稱羅浮或拉號,是霞雲坪以南面積最大的河階。從巴陵到 合流的大漢溪河谷, 明顯表現出岩性與河流地形之間的密切關係。58 採樣處
59 採樣地點如圖 5.1 與 5.2 中紅色星號所標示。土石流源頭地區, 如圖5.3 ~ 5.7 所示。經沿線調查發現,附近岩性十分複雜,沿途出現 頁岩、砂岩、板岩等。甚至於區域一內有二處主要源頭區,由肉眼觀 察顏色的顯著差異,也可判斷為分屬二種不同土石材料,如圖5.4 所 示。本研究採樣一地點為圖5.4 中左側。此區近照如圖 5.5 所示。其 餘之區域二及區域三照片分別為圖5.6 及 5.7 所示,此三區域皆位於 北橫公路沿線。 圖5.3 研究區域一之土石流上游全貌
60
圖 5.4 研究區域一之二處主要源頭之近照
61
圖 5.6 研究區域二採樣區域照片
62
5.3 研究地點地質概況
本研究區屬於雪山山脈地質區的北部,在何春蓀(1997)的台灣 地質圖概論中提到,雪山山脈地質區屬於中央山脈地質區的西面,通 稱雪山山脈帶,參照圖5.8。雪山是台灣第二高山,約位於本帶的中 央,本帶即因其而得名。本帶長約200 公里,平均寬約 20~25 公里。 東北起自北海岸的福隆,向南延經烏來、雪山、埔里和日月潭地區, 到達玉山山脈南邊的荖濃溪的上游為止,全省最高的玉山山嶺也包括 在本地質亞區之內。雪山山脈帶在西邊以屈尺斷層和西部麓山地質區 分隔,在東邊以梨山斷層和脊樑山脈帶相隔,參照圖5.9。 圖 5.8 台灣地質分區圖(何春蓀 , 1974)63 圖5.9 雪山山脈帶與西部麓山帶和脊樑山脈帶分界圖 (何春蓀 , 1974) 本研究區域的地層為中新世至漸新世的地層,參照圖5.10,其敘 述如下: 【木山層】(MF) 屬中新世地層,本層以厚層之白灰色、黃灰色及白色,細粒至粗 粒砂岩及砂泥岩為主,夾有灰色頁岩,或細粒砂岩及頁岩之互層,偶 夾有薄煤層。砂岩含少量黏土礦物。部分地區有玄武質之火山凝灰岩。 【大寮公館層】(TK) 屬中新世地層,本層以深灰色至黑色頁岩為主,中段夾有灰色細 粒泥質砂岩,部分地區砂岩含有鈣質,常造成懸崖或山脊,部分地區
64 夾凝灰岩。頁岩中偶夾有薄層細粒砂岩,頁岩風化後呈洋蔥狀之構造。 【媽岡層】(MK) 屬中新世地層,本層應歸為澳底層之上段,由深灰色頁岩夾有灰 色細粒砂岩薄互層組成。這些岩石只受到輕微的變硬作用,葉理結構 並不明顯。 【石底層】(ST) 屬中新世,以厚層或中厚層,淺灰色或黃灰色細粒至中細粒砂岩 為主,夾有薄層頁岩及砂頁岩互層,岩質堅硬緻密,層理發達。 【大桶山層】(TTS) 屬漸新世地層,本層以黑色硬頁岩為主,夾有灰至灰黑色細粒泥 質變質砂岩,所夾之變質砂岩較乾溝層為多,堅硬緻密的泥質粉砂岩 抗蝕力強,常沿河床形成陡壁。此層在許多地方含有狹小的玄武質火 山碎屑岩,或玄武岩流,通常成為不規則體。 【粗窟砂岩】(TSK) 在雪山山脈帶北部,大桶山層的下部岩層中有一厚砂岩段,砂岩 是暗灰色、泥質、細粒,並含有少許硬頁岩的夾層,此砂岩段後來被 命名為粗窟砂岩,屬漸新世,可以成為大桶山層和其下乾溝層分界的 依據。
65 【乾溝層】(KK) 屬漸新世地層,本層主要是由黑色硬頁岩夾有灰黑色泥質緻密細 粒變質砂岩組成,硬頁岩一般呈厚層狀。本層層理較不清楚。以粗窟 砂岩段和大桶山層分界。 大致而言,霞雲坪至高坡及蘇樂至巴陵間為砂頁岩互層區。而大 灣至高義間及巴陵附近為黑色頁岩與板岩狀黑色頁岩。 圖 5.10 北橫公路復興至巴陵段研究區域地質圖(蔡允瀚,2002)
66