鄭子璉 成功大學水利及海洋工程研究所 研究生 |
周乃昉 成功大學水利及海洋工程學系 副教授 |
關鍵詞:高度平衡多邊形、三角網、平均雨量、集水區 摘要 估算集水區平均降雨量常採用算數平均法、徐昇多邊形法與等雨量線法進行分析,惟高度平衡多邊形法以雨量站高程輔助分配其雨量站之權值,一般認為此法於地形崎嶇地區較徐昇多邊形法為佳,但因計算過程繁雜,以致分析不易,其應用亦相對減少。 本研究發展之電子計算機高度平衡多邊形法的幾何計算模式,支援中央大學太空遙測中心之數值高程資料檔及 ArcView 圖檔,可迅速依序求得雨量站之三角網、高程中點三角網及高度平衡多邊形網,求出各雨量測站在集水區內的控制面積加權因子,有助於現場即時估計平均降雨量。研究中建議了兩項高度平衡多邊形法補充修正的原則, (1) 宜採高程加權中點取代標高中點; (2) 宜以自相交點作對角平分線修正發生自重疊的高度平衡多邊形。 |
Tze-Lien Tseng | Frederick N.F. Chou |
Department of Hydraulics & Ocean
Engineering National Cheng Kung University Tainan, Taiwan 71048, R.O.C. |
Key words: height balance polygons, triangulated irregular networks, areal rainfall, watershed Abstract The most frequently adopted algorithms in estimating average rainfall within a certain watershed are Arithmetical Averaging Method, Thiessen Polygons Method, and Isohyeatal Method. Because of involvement of geographic elevation factors in the estimation, when evaluating the average rainfall in rugged and craggy area, Height Balance Polygons Method (HBPM) is more accurate than the former three estimation algorithms. However, complicated calculation and analysis procedures make itself less applicable in the routines of estimating average rainfall. In order to estimate real-time average rainfall in the field, a computerized HBPM model was developed. This numerical module quickly processes the triangulated irregular networks (TIN), center of elevation of TIN, height balance polygons and height balance weight factors with DTM data and ArcView format. Besides, two revisions were suggested in this study; (1) using weight central point instead of central point of elevation in the computation process. (2) Applying the angular bisector at the auto-intersection of the height balance polygons. |
在集水區中由點降雨推估全流域之平均降雨有很多方法,一般常用之方法有算數平均法、徐昇多邊形法 (Thiessen Polygons Method) 與等雨量線法等,在此多種分析方法外,另有高度平衡多邊形法 (Height Balance Polygons Method) 可考慮集水區地形高程變化,藉此推估各點雨量對集水區平均降雨量之影響,過去的研究中認為徐昇多邊形法 (Thiessen Polygons Method) 與等雨量線法計算結果相近,適用於坡度平緩地區,高度平衡多邊形法考慮集水區地形高程變化,適用於地勢崎嶇、山嶺險阻之區域,惟高度平衡多邊形網因計算過程繁雜,在多數研究中,較少運用本法直接進行集水區面積雨量估算,甚至進而探討在不同雨量站組合下的情形;而徐昇多邊形法計算過程雖然仍然繁雜,但較高度平衡多邊形法分析為易,且多有相關數值分析模式可直接加以運用,分析結果廣受相關專家學者接受。在高度平衡多邊形法中以高程加權雨量站所分配之控制權重,一般認為徐昇多邊形法僅考慮二維平面權值分配,若在地形崎嶇多變之地區,應用高度平衡多邊形法會較徐昇多邊形法為佳。但高度平衡多邊形法若以人工方式繪製,須耗費相當的時間及人力,且其中若有一雨量站移動、增設或裁撤,則整個高度平衡多邊形網便須重新繪製,若應用於即時防洪預警系統中,可能因電傳雨量站故障導致必須重新計算分析,若以人工作業進行本項分析,將無法配合即時防洪預警系統的時效。
隨著近年來地理資訊系統快速發展,高度平衡多邊形法所需之地理基本資料建構漸趨完整,除可利用衛星定位系統快速取得雨量站座標,並可結合中央大學太空遙測中心所提供的數值高程資料 (DTM) 檔取得集水區內各座標的高程資料,可以此基本資料建立高度平衡多邊形法計算機分析模式。因此在研究中利用電子計算機建立完整高度平衡多邊形法幾何計算分析模式,可迅速完成高度平衡多邊形網的繪製與計算,並計算出每一個雨量測站的高度平衡加權因子,有助於現場估計即時平均降雨量的計算應用。研究中以 Visual Basic 為發展工具,並直接取用中央大學太空遙測中心之數值高程二進位資料檔,在分析上包含雨量站三角網分析、三角網連線上兩雨量站高程中點分析、高程中點三角網分析、高度平衡多邊形網、交集分析與高度平衡加權因子計算。
高度平衡多邊形法與徐昇多邊形法在基礎三角網上分析是相同的,經過不同面積範圍計算後,再經相同的交集分析得到各雨量站的控制面積及加權因子。
國內採用徐昇多邊形法進行相關分析之研究相當多,而探討徐昇多邊形法於電子計算機上之應用分兩大研究方向,一是由計算幾何的基礎下,所發展出的向量分析方法,一是由計算圖學的基礎下,所發展出的網格分析方法。向量分析的優點是其計算結果為解析解,缺點為計算過程繁雜,不易撰寫分析模組;網格分析的優點是計算過程單純,易撰寫分析模組,缺點為計算結果為數值近似解,計算時間隨網格精度以平方倍增加。
向量分析透過三角網與徐昇多邊形網互為正交之定理,以標準作圖法由三角網推得徐昇多邊形網為依據,相關研究有周乃昉、鄭子璉[1-3]於民國 84 年搭配 AutoCAD 與 AutoLisp 的初步研究,民國 88 年改善計算演算法並擴展成具有獨立人機介面之分析模組,民國 89 年修正交集方法縮短計算時間,增補內點判別原則,並支援 ArcView 圖檔格式。網格分析以任意點至控制點最短距離之基本假設為依據,相關研究有謝錦志等[4]於民國 89 年證明網格分析之依據與徐昇多邊形網基本假設吻合,龔誠山等[5]於民國 89 年配合 ArcView 發展整合性的數值模式,另 ArcView 擴展套件 Spatial Analyst 中徐昇多邊形法 (Proximity) 的相關分析上亦是採用網格分析,唯其分析因透過內建之直譯器進行,以致執行效能不佳。
國內相關著作中,最早敘述說明到高度平衡多邊形法的文獻可能為徐世大等[6]於民國 59 年 6 月所出版之「實用水文學」乙書,較為農田水利從業人員廣為接觸乙書為王如意等[7]於民國 74 年所出版之「應用水文學」,此二份文件在高度平衡多邊形法分析程序及原理相當接近,均僅概念性介紹,於實際運用之說明較少。而在探討高度平衡多邊形法於電子計算機上之應用中,由於高度平衡多邊形法在作圖上需經過高程中點三角網分析,因此,網格分析法難以直接運算分析,僅向量分析法方能加以運算。相關研究有鄭子璉等[8]於民國 88 年配合數值高程資料所提出的數值模式,其中並提出兩項高度平衡多邊形的補充分析原則。
探討徐昇多邊形法與高度平衡多邊形法各雨量站的控制面積及加權因子的一般性比較,國內鄭子璉等[9]於民國 89 年討論比較。
高度平衡多邊形法為先標明各水文站之標高,然後連接相鄰兩站,取其標高之中點處,該標高之中點處不一定為該兩站連接線之中點。再連接各標高中點成甚多個三角形,再做各三角形之內角等分線,交於各三角形之內心。連接各內心可圍住一雨量站,此面積即為該雨量站之標高及所控制之面積,對全流域面積之比例即為高度平衡加權因子,詳細計算方式分述如後。
運用電子計算機進行計算幾何分析計算前,應先透過數位板或其他方式將流域資料進行數位化,以供電子計算機分析使用。流域及雨量站經數化處理後,在往後的各項應用中非常方便,如須計算其它組合的高度平衡多邊形網時,只須將已數化完畢之圖檔讀入,即可立刻進行分析使用,即使測站有增加、遷移或重新測量座標等情形,只須略加修改圖檔後,即可立刻重新分析。鑑於地理資訊系統日漸普及,欲計算之流域亦可能已有數化完成的圖檔,或僅須對不同系統撰寫資料轉換界面,或經由各地理資訊系統將圖檔轉換成可以讀取之格式存檔,再進行分析,因此僅需考量資料數值分析部分,資料數位化之建立、管理或編修可仰賴現有之軟體達成。
高度平衡多邊形三角網中,各三角形的形狀影響高度平衡多邊形法估算加權因子的均勻性及精確度,由不同組的三個雨量站所形成的三角形,以每個三角形都接近正三角形最為理想,但三角形受雨量站之分佈位置受地形及雨量站是否故障等情況影響,實難以達此理想。在本研究中三角網之決定,以採用任三個雨量站所決定的外接圓內不包含其它雨量站者,構成一個基本三角形,利用此法可決定出較為理想的三角形,並可儘量使三角形的各邊長最短,以此原則逐一分析所有測站組合後,即可得到較佳的雨量站三角網。
高度平衡多邊形法在計算分析上的核心即為高程中點的決定,在本研究中以兩雨量站間加權高程中點為平均高程中點,程式分析可從中央大學太空遙測中心所提供的數值高程二進位資料檔中直接擷取高程資料,及搜尋指定目錄下所有子目錄內所需之數值高程二進位資料檔,並由高於較低高程的雨量站至低於較高高程的雨量站間的剖面資料進行高程加權中點之分析。高程加權中點之計算式如方程式 (1) , (2) 。
![]() |
(1) |
![]() |
(2) |
其中,
![]() |
= 高程中點座標 |
![]() |
= 兩雨量站連線的內插高程座標 |
在求取高程中點時,先以雨量站座標位置所處的最小格網對雨量站進行高程檢定,以降低可能因文獻紀錄資料精度不足所造成的誤差,另外在兩雨量站連線的高程擷取上,只內插求取與 DTM 圖檔格網線上的交點座標,以降低空間內插所可能造成的誤差累積,如圖 1 ,在 x 方向與 y 方向交叉點為已知格網座標,內插與 x 方向及 y 方向格線相交點之高程座標為地形剖面上之點並進行加權高程中點計算。
由兩兩雨量站連線已分析完成的高程中點可對應至雨量站連線之三角網上,逐一連接三角網上的所有高程中點所形成新的三角網即為高程中點三角網。
由高程中點三角網內每一頂點作對角平分線可交於各高程中點三角形之內心,並反向延伸各對角平分線至計算邊界,則可形成包含於對應雨量站的閉合多邊形,此即為高度平衡多邊形,所有對應於各雨量站的高度平衡多邊形總稱高度平衡多邊形網。
當高度平衡多邊形網分析完成後,應分析各雨量站所控制之高度平衡多邊形依序與流域多邊形進行交集分析[2]並求出所佔面積,方能計算高度平衡加權因子,在兩非凸多邊形進行交集分析時,先以兩多邊形相交的交點座標,將兩多邊形切成多段之多重折線,並逐一判斷兩多邊形所切成多段之多重折線所有座標點是否在另一多邊形內,分別剔除所有未在另一多邊形內的多重折線,所得之所有座標點即為欲求交集多邊形邊界座標,依序逐一串接交集多邊形邊界座標及交點座標,即可得完整之交集多邊形。可得各雨量站在該流域中控制的閉合多邊形,其中交集分析採自行發展的交集演算法,以提高執行效率及計算精度。
在空間中之一任意點可位在多邊形內、在多邊形外及在多邊形上等三種情況,若在多邊形上又可細分成在線段上或在線段端點上。若將點與多邊形畫在同一圖紙上以肉眼判斷,可以很容易看出點與多邊形的關係,但在數學上要判斷如圖 2 所示之「任意點是否在多邊形內」,可考慮自該點做任一方向射線求取多邊形與射線的交點數,若為奇數則可得知點在多邊形內;若為偶數(包含 0 )則可知點在多邊形外,在一般計算幾何相關參考文獻[11]上多採用此法。但此種方法仍有例外情況存在,例如自點 Pa 的射線恰與多邊形相切於某一端點,或射線與多邊形之一邊恰巧重合,或如點 Pc 或點 Pd 恰巧落在多邊形上,因此前述計算方法多僅用在凸多邊形的問題上,而流域多邊形通常為不自相交之凹多邊形,故此法可能會因相切或重合而不適用。
研究中依幾何學的原理,採「多邊形各頂點依序對任意點求取角度差之總合」來判別任意點是否在多邊形內。計算如下式:
![]() |
(3) |
式中,
q | = 多邊形各頂點依序對任意點求取角度差之總合 |
n | = 多邊形總點數 |
![]() |
= 所需判斷之任意點之座標 |
![]() |
= 多邊形上第 i 點之座標 |
![]() |
=第 i 點與第 i+1 點角度差 |
在式 (3) 中,計算結果可能如式 (4) ,
![]() |
(4) |
研究中發現,當 q 計算結果為
0
時,任意點在多邊形上或多邊形外;當計算結果為
時,任意點在多邊形內;當計算結果為其它值時,任意點在多邊形的端點上。其分析原理可藉由圖
2 來協助說明之。
在圖 2 中,如任意點 Pa 在多邊形內時,多邊形各頂點與任意點的角度差會沿著多邊形旋轉一圈,累加為
,如任意點 Pb 在多邊形外,多邊形各頂點與點
Pb
的角度差累加為 0 ,且角度差和變化量由 0 增加再減少,至最後一多邊形頂點與起點之角度差變化量恰巧合為 0 ;如任意點 Pc
落在多邊形的端點上,由於同一點不能做夾角分析,因此最後計算結果會等於任意點前後兩邊的角度差值。
在這其中如任意點 Pd 位在多邊形的線段上時,其角度差和為 0 ,但在其計算過程中必有一個角度差值為 ,因此當多邊形各頂點依序對任意點求取角度差之總合為 0
時,又可細分為當角度差
計算過程中出現為
時,任意點在多邊形的線段上;否則任意點在多邊形外。
在式 (3) 中,若計算結果為
時,可以另外獲得多邊形座標點排列順序旋轉方向之資訊,當計算結果為正值時,可以推得多邊形座標點排列順序沿正方向旋轉,即為沿逆時鐘方向旋轉;反之,當計算結果為負值時,可推得多邊形座標點排列順序沿負方向旋轉,即為順時鐘方向旋轉,可分析欲求交集多邊形座標點排列順序之關係,在進行依序逐一串接交集多邊形邊界座標及交點座標時,方不致發生交集多邊形座標連結錯亂。
對於一流域多邊形來說,在分析前並無法得知該多邊形之旋轉方向,必須取得任意點在流域多邊形內來判別,對於凸多邊形而言,可採形狀重心為任意點,但流域多邊形屬於凹多邊形,其形狀重心並不一定保證在流域多邊形內,因此研究中由流域多邊形任意點起,依序取三個點求三角形重心,若重心在流域多邊形外,表示取得之第二點為凹多邊形內角大於 π 之頂點,並由此點再依序取三個點求三角形重心,直至重心在流域多邊形內,依幾何學原理,任意多邊形必有 3 個以上頂點之內角小於 π ,故此演算法必能找到任意點在多邊形內。
在研究中若流域多邊形內點計算結果為負值時,令該多邊形先進行反轉後才進行分析;而本研究所產生各多邊形均令其沿正方向旋轉,以要求所有各多邊形均沿正方向旋轉。
在進行兩多邊形的交點分析時,兩線段的交點可採解兩直線方程式的二元一次聯立方程式進行計算,若兩直線平行時為無解,若兩直線恰可交於一點時為唯一解,若兩直線重合則為無限多組解;其狀況可運用高斯消去法演算完成的最後一式判定。但若為唯一解時,並不代表兩線段一定有交點,應增加判別交點是否落在兩線段上而非線段外。
欲檢查任意點是否落在線段上可使用兩種方式且都成立:
1. 空間幾何關係
任意點若在線段上,則任意點座標應在線段兩端點座標之間,任意點之點座標必不同時大於或小於兩端點座標,其表示式如下式:
![]() |
(5) |
式中,
![]() |
= 所需判斷之任意點座標 |
![]() |
= 線段兩端點之座標 |
2. 直線方程式
任意點若在線段上,則任意點座標應在線段所形成之直線方程式上,其表示式如下式:
![]() |
(6) |
本文中兩線段的交點座標已採聯立解兩線段之直線方程式而得,所以交點座標必在兩線段上,因此僅需檢查交點座標與兩線段的空間幾何關係,若交點座標均在兩線段的端點座標內,則兩線段必有交點。
由前述交集分析相關理論及方法後,本研究採用之兩任意多邊形的交集分析步驟如下:
圖 3 為交集分析步驟示意圖,圖 3(a) 為多邊形交集前;圖 3(b) 為所有相交之交點,並沿多邊形各相交點,將多邊形從各相交點切成多段多重折線,虛線表示該段多重折線上之任意頂點經內點判別後得知在另一多邊形外;圖 3(c) 為保留交集多邊形之各段多重折線,並逐一串接;圖 3(d) 為完成交集之多邊形。
根據管晏如[10]提出之多邊形面積之計算方法,簡記如 (7) 式:
![]() |
(7) |
式中,
n | = 多邊形總點數 |
A | = 多邊形面積 |
![]() |
= 多邊形上第 i 點之座標 |
各高度平衡交集多邊形為各雨量站在該流域中控制的閉合多邊形,計算此一多邊形與整體流域多邊形的比例面積即為該雨量站的高度平衡加權因子:
![]() |
(8) |
式中,
![]() |
= 高度平衡加權因子,無因次; |
![]() |
= 雨量站的高度平衡多邊形在集水區內面積; |
A | = 集水區流域多邊形總面積。 |
本研究採曾文水庫集水區為計算案例[12]。曾文水庫集水區集水面積為 481.16 平方公里,集水區內外目前共有九個電傳雨量站,在即時防洪運轉過程中,由各電傳雨量站整點傳回之雨量資料可估計集水區內的平均降雨,並可推估未來數個時段內水庫可能的進水量,以此資料可進行水庫的防洪運轉操作。若在即時防洪運轉過程中發生電傳雨量站故障,應立即依可正常作業中之電傳雨量站重新分析各雨量站的高度平衡加權因子,以利推估平均降雨,若能應用電子計算機進行分析,將可大幅提高分析效率。
為求精確計算各雨量站所控制之面積,在研究中利用全球定位系統 (GPS) 至現地各雨量站進行座標校正,並以校正後之雨量站座標進行分析,如表 1。
5.1 分析結果
圖 4 (a) 為已分析完成的程式螢幕展示畫面,圖面展示了分析完成之三角網、高程中點三角網、高度平衡多邊形網及與集水區邊界多邊形交集完成的控制區域,以 Pentium 233 MMX 計算此案例約共需 10 秒 (雨量站 9 站,邊界多邊形 688 點) ,三角網高程中點三角網及高度平衡多邊形網分析部分約佔 2 秒,高度平衡多邊形網與流域多邊形交集約需 8 秒,圖 4 (b) 為高度平衡加權因子輸出表之程式螢幕展示畫面,圖 4 (b) 表格中高度平衡網交集面積欄之單位為平方公尺。
就模式計算之時效性來說,純就人工圖紙作業,約需三到七個工作天,若以本模式計算與操作,包含流域多邊形及雨量站位置數化與模式分析計算,整體作業可在 15 分鍾內操作完成。所運用的編譯器為 Visual Basic 5.0 中文專業版,分析模組可在一般慣用之 32 位元之視窗環境系統上執行 (Win32) 。
5.2 研究討論
茲就本研究進行時所遭遇疑義部分提出討論。
高程中點在原文獻[6,7]定義上是取兩點間高程恰為兩點高程平均值之點,惟在本研究中計算案例電傳雨量站因考量到無線電波傳送,一般均設置於地形稜線上,因此若以此定義之高程中點分割,將使高度平衡加權因子分配十分不平均。例如圖 5 (a) ,在表湖與三角南山連線的剖面上,兩點高程平均點與整個剖面長度的比例約為 0.5 % ,若以此定義進行高度平衡多邊形網計算,在研究中認為並不合適。且高程中點若以兩點高程平均值之點為中點,則亦有可能出現三個點以上反而造成應用上的困擾,如圖 5 (b) ,在大棟山與樂野連線的剖面上,兩點高程平均點共計 11 個點,在高程中點反而難以選取,故本研究中改以高程加權中點進行分析,除解決前述選點的困擾外,並可始做圖有一致標準。
高程三角網計算中,以相同之兩雨量站連線為邊之三角形有 1 ~ 2 個,因此高程中點計算會造成重複分析,研究中考慮計算時效,將所有兩雨量站組合情形納入分析,並先將各種組合計算完成存檔,在高程中點分析以查表方式直接取得計算結果,以加速分析速度,此法對於所有高度平衡多邊形組合分析有極大之加速效果。
靠近圖面邊緣的雨量站之高度平衡多邊形多無法閉合,而對角平分線反向將無限延長至無窮遠,因此不但不能計算高度平衡多邊形面積,也無法執行交集分析。為解決此種情況,本研究依據人工圖紙作業以有限圖紙分析高度平衡多邊形網之方式,由雨量站位置座標點與集水區邊界之資料點中找出最小及最大的 x 值與 y 值,並向圖面外延伸各軸邊長的 5% 做為分析範圍邊界的頂點,憑以繪一矩形視為圖紙邊緣並當做分析範圍邊界。再針對各高度平衡多邊形逐一檢查其是否閉合,非閉合的高度平衡多邊形則由多重折線的起點與終點處向外延伸至分析範圍邊界,依正方向(逆時針方向)繪成一閉合多邊形;若原為閉合的高度平衡多邊形則與分析範圍邊界進行交集分析,以交集完成之高度平衡多邊形取代超過分析範圍邊界的高度平衡多邊形。附帶的,加入分析範圍邊界亦可使輸出之圖形較美觀,可依圖面內容完全展示在螢幕或圖紙上。
前述高度平衡多邊形雖與分析範圍邊界進行交集分析,但實際運用分析時,有可能因雨量站座標位置發生高度平衡多邊形內自相交的情形,如圖 6 ,而此高度平衡多邊形相鄰兩多邊形亦必發生重疊,過去參考文獻中並無討論到此情形之處理方式,在本研究中沿用高度平衡多邊形內心分割原則,令從相交點作一對角平分線延伸至分析範圍邊界,憑此線段修正共用這兩相交線段之三個多邊形。
若雨量站網由單一雨量站構成,則無三角網、高程中點三角網及高度平衡多邊形網,令高度平衡加權因子為 1 ,若雨量站網由兩雨量站構成,則三角網為單一直線,無高程中點三角網,高度平衡多邊形網由與兩雨量站連線垂直並通過高程中點之直線與分析範圍邊界構成。
內點判別在向量計算佔相當高比例之計算,若能改善內點判別之演算效率,對整體分析效率應能有相當之提昇,原內點判別演算法為每一小角度差均需經兩次反三角函數演算,若以欲判別之內點為相對座標之原點,則在此相對座標內之任一象限內之角度和應為象限內多重折線頭尾兩點之角度差,以此觀念可降低反三角函數之計算量,如圖 7 中,虛線視為直線兩端所得之角度差與視為多重折線逐一累加相同,而不跨越相對座標之任一象限多重折線,必不影響及包含內點座標位置,故內點判別可不受影響。
然上述簡化必須經象限判斷,若此處之分析演算不能快於反三角函數運算,則無加速之能力,在 32 位元視窗系統 (Win32) 中,以四個位元組 (byte) 之長整數運算效能最高,且邏輯運算元運算速度遠高於反三角函數運算,研究中參考陳雪美[13]譯著,直線的二度空間剪取演算法之邏輯,修改並提出以內點及等於內點座標所形成之空間為相對原點,逐一運算多邊形頂點對內點位元值可得描述特定空間之長整數,長整數共計有 32 位元 (bit) ,將 16 個高位元劃作等號成立,將 16 個低位元劃作大於判斷成立,可處理 16 個座標軸之分析,其計算流程如下:
For t =1 To t_Count Select Case P(t) Case Is > P0(t) cs = cs Or 2t Case Is = P0(t) cs = cs Or 2t+16 End Select Next |
流程中,
cs | = 特定空間之長整數 |
t_Count | = 座標軸總數 |
t | = 正在分析之座標軸 |
P(t) | = 正在分析 t 軸之值 |
P0(t) | = 正在分析 t 軸之內點座標值 |
而上述流程 cs 為利用長整數進行位元運算,其值不需經任何運算即可直接得到長整數,長整數電腦運算方式如下:
![]() |
(9) |
式中,
bi | = 該位元之係數值, 1 或 0 ; |
研究中僅運用到 2 維空間判別,故其運算結果必為圖 8 中 9 種情況所得長整數值之一。透過前述判別原則,可得一與多邊形各頂點對應之特定空間長整數陣列,在判別各象限連續之多重折線可利用前後兩點特定空間長整數值是否相同來判別。在分析中,若有任意端點其值為 196608 ,則代表該任意點與多邊形頂點重合,若前後兩端點未直接跨越對角象限(I 、 III 象限或 II 、 IV 象限),則可省略是否在多邊形上之判斷。
經修改加速之演算法,以獨立測試環境比較與原演算法之計算效率差,利用 4 組不同之起始亂數種子與分析範圍邊界產生各 1,000 點隨機亂數內點,與流域多邊形 688 點進行內點判別運算,研究中以 Pentium 233 MMX 將計算所需時間列於表 2 中,平均執行效能提高約 3 倍,可有效減少計算所需時間。
在現有運用地理資訊系統分析中,常會因圖檔格式眾多,且不同軟體彼此之間圖檔格式互不相容,而造成應用上的困擾及資料轉換或重新數化上時間的浪費,因此本研究中所發展的分析模式除可使用自行定義的 ASCII 純文字格式,亦支援 ArcView 的 Shape File 格式[14],幾何計算模式可直接匯入 ArcView 的圖檔,如圖 9 ,計算完成之圖形可匯出為 ArcView 的圖檔,供其他分析使用,如圖 10 ,經此處理後,可使分析結果之應用更加多樣化,此外,本研究中所發展的分析模式亦可直接匯出視窗加強型向量檔 (Enhanced Metafile, EMF) ,此部份較 ArcView 匯出之 EMF 檔案技術穩定,適合配合文章及簡報檔之製作,不似 ArcView 會有部分圖檔邊緣顯示不完全的情況。
在套裝軟體 ArcView 的交集分析功能中,會在交集邊界上形成階梯狀,格網方式分析對多邊形的邊界點數增加很快,在運用 ArcView 進行圖元分析原數化座標點為 688 點之多邊形,經過數次交、聯、差集分析,多邊形已迅速增為 2,203 點,又在 ArcView 中需設定座標解析度,此一設定亦限制其在交集分析上的交點座標必位於格點上,雖然其座標可以以倍精度儲存,但仍無法解析比格網還精確的座標,計算求得之交點座標必會經過適當的進位至格網座標精度方能順利分析,在本研究所發展的幾何計算模式僅增加多邊形交點座標,其餘座標均為原多邊形數化座標,交點座標完全使用倍精度可用的最多位數表示,在分析及應用上座標值失真較小。
徐昇多邊形法僅考慮二維平面的控制面積,在做圖處理亦以徐昇多邊形法較為單純,且目前此法已有多套的的現成計算機模組,因此徐昇多邊形法較常被採用於面積雨量的推算,在本研究中亦比較徐昇多邊形法與高度平衡多邊形法的分析成果。
徐昇多邊形法與高度平衡多邊形法計算分析圖如圖 11 及圖 12 ,加權因子如表 3 ,高度平衡多邊形法需參考雨量站連線剖面之資料,再據以分析,在作圖結果中可以看出,高度平衡多邊形網受集水區邊界高程影響,在集水區邊界上的雨量站有較大的控制範圍,徐昇多邊形網則無此影響。由於計算方法不同,對於 9 站電傳雨量站均選取分析的狀況來說,徐昇加權因子與高度平衡加權因子或有增減,分析結果趨近一致。
對於接近線性分布的雨量站網,如圖 12 及表 3 中案例 1 、 2 中曾文新村、樂野、馬頭山及龍美所構成之雨量站網,徐昇多邊形法分析結果為接近平行分布之分割,高度平衡多邊形法分析結果則為接近T 形分布之分割,徐昇加權因子尚維持均勻分配,高度平衡加權因子則在線性站網兩端端點雨量站佔大部分之權重,若雨量站網接近ㄩ或ㄇ字型分布,如圖 12 及表 3 中案例 3 、 4 ,則亦會使高度平衡加權因子在站網兩端端點雨量站佔大部分之權重,由圖 12 中,站網兩端端點雨量站在高度平衡多邊形法僅會有一條分割線經過,徐昇多邊形法則會由數條分割線經過,徐昇多邊形法在加權因子分配較合理,如表 3各計算案例中,曾文新村徐昇加權因子介於 2 % ~ 7 % ,曾文新村高度平衡加權因子介於 37 % ~ 52 %,曾文新村加權因子相差達 5 ~ 21 倍。
若將ㄩ或ㄇ字型分布的雨量站網視為凹形分布的雨量站網的正反集合,線性分布的雨量站網視為凹形分布的退化雨量站網,則可推知,高度平衡多邊形法受雨量站網分布位置影響頗大,在運用此法分析流域平均降雨量時應盡量納入所有可分析之雨量站進行分析,以提高平均降雨量計算之精度,當雨量站網成凹形分布的狀況下,應避免以高度平衡多邊形法計算平均雨量。
在 Visual Basic 中,預設變數資料型態為 Variant 型態,其變數精確度在實數方面可比照倍精度變數,因此可知至少具有 15
位以上的精確度,而在判斷二線段是否重合、相切,或點是否相同時,有時會因受到計算機變數的誤差限制,在數學上或實際上是同一點的情況中,有可能在計算機中會有極小不同的誤差,例如將點座標代入直線方程式時,可能因為點座標或直線方程式之係數為無理數,而無法以有限位數表達,因此本研究中選用相對誤差做判斷標準,並取
10-10 為容許誤差的門檻值,藉此判別是否為同一點或某點是否在目標直線方程式上。
台灣地區二度分帶座標值位數龐大,但實際流域邊界內座標變動位數小,分析計算案例曾文水庫集水區座標範圍約在 (200000, 2568000) ~
(232000, 2601000) 之間,以 y 方向之座標軸來說,小數點以上即有 7 位數,若採用容許誤差為絕對精度 10-10
,則容許誤差應用於本研究案例,會比計算機可使用 15 位有效位數還小 2 個值級,形成容許誤差完全失效,因此必須改採用容許誤差為相對精度 10-10
,並於計算分析時改採相對座標,將座標值未變動之值級部分進行座標平移,座標平移所採用之相對座標原點如下式:
![]() |
(10) |
式中,
dt | = t 軸方向上之差量位數所得主間距量; |
mt | = t 軸方向上之相對座標原點偏移量; |
轉換分析計算案例曾文水庫集水區座標範圍變更為 (0, 68000) ~ (32000, 101000) 之間,而在實際對外部進行輸入輸出時才採用絕對座標,提高電子計算機有效位數的運用。
在檢查任意點是否在線段上所使用兩種方式均有容許精度的問題,在空間幾何關係分析上,判別兩點是否為同一點時,可能因有效位數的影響,造成同一點可能座標值略有不同,這對串接多重折線成一多邊形時可能造成困擾,在代入直線方程式分析上,可能無法很精準的恰巧等於 0 ,因此配合採用容許誤差為相對精度 10-10 ,使得電子計算機因有效位數對計算所造成的誤差能消除。
以高斯消去法解得之兩線段交點會因兩直線方程式的求取過程而有可能稍有差異,如直線方程式在計算常係數時代入不同的線段端點座標,則有可能發生常係數在有效位數最後 1 ~ 2 位數不完全相同的情況,又兩直線方程式在求解的係數矩陣順序上的不同,亦有可能解出在有效位數最後 1 ~ 3 位數略有差異的交點座標,在高斯消去法分析之前可先對輸入之係數矩陣依序排列,以使同一組直線方程式組必解出相同交點座標,另需對所產生之直線方程式係數約分,如此可使兩直線方程式之係數在同一值級上,可避免在消除係數的步驟時,誤差位數過度放大而影響正確的計算。
經此誤差修正後,多數計算誤差在轉回原二度分帶座標後,均被擠出有效位數之外,在 9 站電傳雨量站均選取之情況下,誤差為有效位數不足所造成,如圖 4 (b) ,加權因子累積誤差小於 10-14 ,實際面積計算誤差僅為 3.80 平方公釐,約為總面積 10 兆分之一,遠小於座標平移前之誤差。
在流域邊界數化一般可分人工數化及自動數化,人工數化可能因流域邊界的線寬或流域邊界本身的精度及端點選點的數量,而導致產生數化誤差,自動數化則受限於演算法及資料本身的精度與解析度而導致產生數化誤差,高精密度的流域邊界數化受限於原始資料,最大可能產生約為百分之一的面積誤差,研究中曾文水庫集水區數化面積為 482.67 平方公里,如圖 4 (b) ,略大於文獻 481.16 平方公里,數話誤差為 0.31 %,在研究中認為此誤差應較利用求積儀所描繪之面積誤差為小,並不影響高度平衡加權因子之計算。
在雨量站座標數化來說,多由水文年報或文獻等記載雨量站座標經轉換後點入,亦可由相片基本圖上判讀出來,過去文獻上由於測量儀器簡陋,導致文獻記載雨量站座標精度不足,多有座標為 0 分或 0 秒的紀錄,部分座標經轉換後應在流域內的雨量站反而位於流域外,誤差較大。在曾文水庫集水區中,電傳雨量站經全球定位系統 (GPS) 校正後,與文獻資料最大相差距離達 2,544 公尺,詳如表 1 ,全球定位系統一般可達到 15 公尺以下的精確度,有助於快速檢核雨量站座標是否正確。雨量站座標數化精度對高度平衡加權因子計算影響視週遭各雨量站位置而定,在曾文水庫集水區中,加權因子最大產生約為 0.1 的控制面積誤差。
另在本研究中並修正前研究之交集分析方法,有效利用計算機變數可用的位數,使計算誤差均降至小於相對精度 10-14 ,使計算結果更為可靠。相較利用 ArcView 進行交集分析,可能受 ArcView 格網分析使得各計算區多邊形邊界形成階梯狀邊界,而使計算誤差大於相對精度 10-4 ,在本研究中所運用的交集分析較 ArcView 可靠,並可維持原數化精度。
前述所提之數化誤差遠高於模式計算誤差 11 個值級以上,因此模式之計算誤差均可忽略不計,研究中認為特別應提高雨量站座標數化精度,以提升高度平衡加權因子計算精度。
在本研究期間,承蒙曾文水庫管理中心洪燈河先生協助,以全球定位系統針對計算案例之雨量站進行座標位置校正,中央大學太空及遙測中心饒見有研究員提供 DTM 圖檔格式資料,仲琦科技總公司林開輝先生及高雄分公司彭美香小姐、蔡培元先生提供 ArcView 圖檔格式資料,謹致謝意。
為尊重智慧財產權,特將本文所提及之各項軟體、商標及所屬公司名稱列出,以示尊重。列出說明如下:
本研究成果所發展之幾何計算模式可於網際網路上自由下載,網頁位置為:
http://feitsui.hyd.ncku.edu.tw/TLCheng/Thiessen/
表 1 曾文水庫集水區雨量站位置座標 (東經 120 度、北緯 23 度,表內單位:分 - 秒)
站名 | 相差距離 (公尺) | 文獻位置 | 修正位置 | ||
東經 | 北緯 | 東經 | 北緯 | ||
曾文新村 | 1,092 | 30-00 | 13-20 | 29-21 | 13-21 |
水山 | 779 | 48-45 | 28-28 | 49-15 | 28-17 |
樂野 | 1,481 | 43-00 | 28-17 | 43-30 | 27-48 |
里佳 | 1,163 | 43-23 | 23-07 | 42-43 | 23-17 |
表湖 | 392 | 39-09 | 15-43 | 39-23 | 15-43 |
馬頭山 | 564 | 35-23 | 20-19 | 35-43 | 20-17 |
龍美 | 2,477 | 38-00 | 24-00 | 39-15 | 24-42 |
三角南山 | 1,946 | 34-00 | 13-00 | 34-48 | 13-45 |
大棟山 | 2,544 | 30-00 | 18-00 | 31-22 | 18-35 |
組別 | 修改前 (秒) | 修改後 (秒) |
1 | 25.65039 | 8.738281 |
2 | 25.92969 | 8.730469 |
3 | 25.64844 | 8.730469 |
4 | 25.65039 | 8.728516 |
平均 | 25.71973 | 8.731934 |
站名 | 曾文新村 | 水山 | 樂野 | 里佳 | 表湖 | 馬頭山 | 龍美 | 三角南山 | 大棟山 | |
徐昇加權因子 | 9 站 | 0.0054 | 0.0992 | 0.1352 | 0.1813 | 0.1353 | 0.1780 | 0.1043 | 0.0809 | 0.0805 |
a-1 | 0.0726 | - | 0.2882 | - | - | 0.4167 | 0.2224 | - | - | |
a-2 | 0.0726 | - | 0.4050 | - | - | 0.5223 | - | - | - | |
a-3 | 0.0240 | 0.3607 | - | - | 0.0000 | - | - | 0.1691 | - | |
a-4 | 0.0359 | 0.1004 | 0.3046 | - | - | 0.4468 | - | - | 0.1123 | |
高度平衡加權因子 | 9 站 | 0.0123 | 0.1081 | 0.1316 | 0.1157 | 0.1775 | 0.1841 | 0.1071 | 0.0731 | 0.0905 |
b-1 | 0.4636 | - | 0.3910 | - | - | 0.0222 | 0.1232 | - | - | |
b-2 | 0.3788 | - | 0.3798 | - | - | 0.2413 | - | - | - | |
b-3 | 0.5232 | 0.2717 | - | - | 0.1817 | - | - | 0.0234 | - | |
b-4 | 0.4284 | 0.149 | 0.1824 | - | - | 0.1512 | - | - | 0.0891 |
圖 1 高程剖面取點示意圖
圖 2 任意點與多邊形的關係
圖 3 任意多邊形交集分析示意圖
圖 4 高度平衡多邊形法計算結果畫面展示
圖 5 兩點高程平均中點
圖 6 自相交之高度平衡多邊形
圖 7 任意點與多邊形的相對座標
圖 8 各特定空間之長整數值與內點座標 (x, y) 之關係
圖 9 匯入 ArcView 圖檔
圖 10 匯出後利用 ArcView 展示畫面
圖 11 徐昇多邊形法計算結果展示
圖 12 徐昇多邊形法與高度平衡多邊形法比較