等值線光滑算法及其適用性分析

王永東，韓尚煒，謝小菊

（1.武漢中迪聯創科技有限公司，湖北 武漢 430074；2.中國地質大學（武漢） 地理與信息工程學院， 湖北 武漢 430074）

等值線是由某一數量指標值相等的各點連成的平滑曲線，能反映該指標在一定空間或時間范圍內的分布規律，被廣泛應用于各種圖件中。等值線的繪制一般分為離散數據格網化、等值線追蹤和等值線光滑 3個步驟，其中離散數據格網化是將數據以格網的形式重新組織，為等值線追蹤做準備，常見的格網類型包括不規則三角網和規則格網；等值線追蹤是針對某一指標，利用格網找到若干等值點，并根據格網拓撲關系確定點的連接順序，從而生成折線型等值線；等值線光滑是對折線型等值線進行光滑處理，使等值線更加符合實際、更加美觀。目前，等值線光滑方法主要包括手工勾繪和計算機自動生成兩類，手工勾繪方法的圖形美觀、精度較高，但費時費力，一般多用于復雜地形圖的出版；計算機自動生成主要依靠光滑曲線生成算法，通過在已知等值點間添加曲線點的方式實現等值線的光滑。目前的光滑曲線生成算法可以分為兩類：①過點的光滑算法，即生成的曲線嚴格經過等值點，以張力樣條函數法為代表[1]，由于曲線必須過 等值點，當等值點較密集時，易產生震蕩現象，與實際地形特征不符合；②不過點的光滑算法，即光滑過程中曲線會經過等值線附近，但不一定嚴格經過等值點[2]，以Bezier函數法為代表[3-4]。

針對Bezier函數法存在的問題，本文采用在拐角處以等值點和切線點共同作為控制點的方法使光滑曲線經過所有等值點[5]，并對比分析了五點法、張力樣條函數法、改進后的Bezier函數法在不同類型地形圖中的適用性。

1 常用的等值線光滑算法

1.1 五點法

五點法的基本原理為：在兩個相鄰數據點之間建立一個三次多項式曲線方程，且要求整條曲線具有連續的一階導數以保證曲線光滑，每個點的一階導數由該點及其相鄰兩側各兩點（共5個點）來確定。對于離散點序列中相鄰兩點(xi,yi)、(xi+1,yi+1)，經過這兩點的光滑曲線可表示為：

為保證曲線光滑，要求從(xi-1,yi-1)到(xi,yi)以及從(xi,yi)到(xi+1,yi+1)這兩個區間上的曲線在(xi,yi)點處具有相等的一階導數。假定已知(xi,yi)和(xi+1,yi+1)處的導數，結合已知點(xi,yi)、(xi+1,yi+1)，代入式（1）可求解參數c0、c1、c2、c3，進而得到從離散點(xi,yi)到(xi+1,yi+1)的曲線方程。

若已知曲線上5個點以及曲線經過各點的順序，可由Akima五點法估算曲線在第三點處的導數值[6]。設平面上5個點為Pi(xi,yi)，i=1,2,3,4,5，則每相鄰兩點間折線段的斜率為：

第三點處斜率為：

若等值線為開曲線，則首尾處各有兩點無法找齊計算所需的5個點，此時采用補點的方法進行處理。具體方法為：對于等值線末端最后3個點(x1,y1)、(x2,y2)、(x3,y3)，補充點(x4,y4)、(x5,y5)滿足式（4），且5個點同屬于式（5）所表示的拋物線上。

通過式（4）、式（5）可求解拋物線方程，進而求解補點坐標。通過補點可計算得到等值線首尾4個點處的一階導數值。

1.2 張力樣條函數法

1.2.1 張力樣條函數

設平面上一組已知點為(x1,y1)、(x2,y2)、…、(xn,yn)，x1＜x2＜…＜xn，常數σ≠0為張力系數，求解一個二階導數連續的三次樣條函數應滿足：

式（7）是一個二階非齊次的常系數線性微分方程，令

求解微分方程可得：

式（9）為通過所有等值點的張力樣條函數，式中 未知項為各等值點處的二階導數值，因此若能確定 各等值點處的二階導數值，則可確定張力函數。

1.2.2 節點關系式與端點條件

對式（9）進行微分，再結合式（10），可得到節點關系式如式（11）所示。

節點關系式可簡記為：

式（12）是含有n個未知量，i=2,3,…,n?1的n-2個方程的線性方程組，要唯一定解，需附加兩個方程式，由于等值線會出現開曲線和閉曲線兩種情況，需要分開討論。

對于開曲線，首尾端點處則有：

由式（12）和式（14）組成三對角線性非奇異方程組，可求得唯一解，再代入式（9）可求得張力樣條 函數方程。

對于閉曲線，首尾端點處則有：

可得到一個三對角線并帶有右上角和左下角系數矩陣的線性方程組，它也是非奇異的，求得唯一解后代入式（9）求得張力樣條函數[7]。

1.2.3 張力系數選擇

在地形圖上采用等高線描繪地形時，特征點之間的曲線取最短為優，為避免相鄰光滑曲線相交且要通過等值點，采用規范化的張力系數σ作為標準值，即

σ為規范化張力系數，需要預先進行試驗，取最適合當前地圖的張力系數，當發生變化時，則由標準值重新計算張力系數[8-9]。

1.3 Bezier函數法

1.3.1 Bezier曲線插值

Bezier曲線是一條需要在起點和終點之間構建插值多項式的混合函數的不規則曲線，可擬合任意數目的控制點，其曲線方程次數比控制點個數少1。通過n+1個頂點定義一個n次多項式，假設n+1個點的位置矢量為Pi，i=0,1n，則Bezier曲線上各點的坐標插值為[10-11]：

式中，Pi為原特征多邊形的等值點；Bi,n(t)為n次Bezier 多項式。

Bezier函數法的缺陷是光滑后的曲線不再通過所有等值點，為解決該問題，本文首先對曲線進行分段擬合，每段取相鄰3個等值點，擬合得到的曲線將不再通過中間的等值點，再利用切線抹角法進行處理，使光滑后的曲線通過所有等值點。

1.3.2 切線抹角法

如圖1a所示，對于折線B-A-C，利用Bezier函數法生成的光滑曲線不再經過A點，為解決該問題，先過點A作線段P1P2，A為P1P2中點（圖1b），再以P1、A、P2共同作為切線抹角法的控制點，即可得到經過 點A的光滑曲線，如圖1c所示。

圖1 切線抹角法原理圖

為避免出現多余拐點，切線P1P2應垂直于∠A的平分線，設∠A的平分線與BC相交于D，，，進而求得切線的斜率為：可計算得到

曲線的松緊程度由P1P2兩點間的距離決定，如 圖2所示，延長DA到E點，取AB與AC的中點F、G，連接FE、GE與切點交于P1、P2，建立系數控制曲線松緊程度。

圖2 切線端點計算

設dt已知，則點E坐標的計算公式為：

基于點F、G、E的坐標以及直線P1P2方程，通過直線求交可計算得到點P1與P2的坐標值[5,12]。

2 算法編程實現與結果分析

本文主要針對地形圖中生成的等值線（即等高線）結果進行分析，測試數據來源于對一些具有代表性地形的實際測量數據。

2.1 五點法測試

為實現五點法，首先在程序中獲取等值線的條數以及每條等值線的點數，然后根據曲線的開閉決定是否需要在首尾各插補兩個點，最后循環建立所有相鄰兩點的三次曲線方程。

五點法測試結果如表1所示，可以看出，五點法的光滑效果主要受等值點密度和等高線密度的影響，當等高距較小、等高線較密集且等值點較稀疏時，光滑處理后的等高線易出現相交的情況。其原因在于如果相鄰兩條等高線走向相背，進行光滑處理時，生成的曲線會向二者中間凸起，若等值點分布稀疏，則會導致凸起部分相交（圖3a），中部的兩條等高線向二者中間凸起，若繼續縮小等高距，則會導致兩條等高線相交（圖3b）。

表1 五點法測試結果

圖3 不同等高距下五點法結果（局部放大圖）

2.2 張力樣條函數法測試

張力樣條函數法首先采用經驗值作為張力系數，再根據等值線的開閉，調用不同的計算函數生成等值線。若光滑效果不佳，則調整張力系數重新進行等值線光滑。在計算等值線的過程中，涉及帶有右上角和左下角的三對角線系數矩陣的線性方程組求解時，可采用追趕法（Doolittle分解）將矩陣分解為一個特殊的下三角矩陣和一個特殊的上三角矩陣的內積，再進行求解。

張力樣條函數法的光滑效果可通過張力系數進行調節。結果表明，當張力系數增大時，曲線會更加貼近原始折線；當張力系數減小時，曲線彎曲程度會增大。考慮到實際應用場景，當等高距較小、等值點密度不夠時，應適當增大張力系數以避免等值線相交，當等高距較大、等值點密度有余時，應適當減小張力系數使地形特征更加明顯。張力樣條函數法也可通過改變張力系數來調整等值線，有效解決了五點法中等值線相交的問題，但該方法首先需要提供經驗值作為張力系數，且需要多次調節張力系數使效果達到最優，因此操作復雜、使用門檻較高。

2.3 改進后的Bezier函數法測試

采用改進后的Bezier函數法生成光滑曲線需分段進行，每次取相鄰3個點坐標，首先計算切線點P1、P2的坐標，然后利用3個等值點和2個切線點生成Bezier曲線。若生成曲線效果欠優，可調整相關系數dt后重新進行等值線光滑。

與張力樣條函數法類似，改進后的Bezier函數法利用dt來控制曲線與原始折線的偏離程度，當dt趨近于1時，曲線趨向原始折線；當dt較小時，曲線會偏離原始折線。因此，當等高線較密集時，應適當增大dt以避免等值線相交；當等高線較稀疏時，應適當減小dt使地形特征更加明顯。與張力樣條函數法類似，改進后的Bezier函數法也可通過調節系數值來使最終效果達到最優。同時，相較于張力樣條函數法，Bezier曲線不會產生多余拐點，光滑效果更優。

3 算法特點和適用性

五點法的計算速度快，但可能出現等值線相交等情況，適用于原始數據質量高、數據分布平緩、等值線稀疏、結果要求不嚴格的場合，操作簡單，可一鍵生成；張力樣條函數法的等值線嚴格正確，但需手動調節張力系數，適用于結果要求嚴格、結果需逐步調整以達到最優狀態的場合；改進后的Bezier函數法的等值線嚴格正確、不產生多余拐點，但需手動調節相關參數，適用于結果要求嚴格、結果需逐步調整以達到最優狀態、不產生多余拐點的場合。