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

王永東，韓尚煒，謝小菊

（1.武漢中迪聯(lián)創(chuàng)科技有限公司，湖北 武漢 430074；2.中國地質(zhì)大學(xué)（武漢） 地理與信息工程學(xué)院， 湖北 武漢 430074）

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

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

1 常用的等值線光滑算法

1.1 五點法

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

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

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

第三點處斜率為：

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

通過式（4）、式（5）可求解拋物線方程，進(jìn)而求解補(bǔ)點坐標(biāo)。通過補(bǔ)點可計算得到等值線首尾4個點處的一階導(dǎo)數(shù)值。

1.2 張力樣條函數(shù)法

1.2.1 張力樣條函數(shù)

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

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

求解微分方程可得：

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

1.2.2 節(jié)點關(guān)系式與端點條件

對式（9）進(jìn)行微分，再結(jié)合式（10），可得到節(jié)點關(guān)系式如式（11）所示。

節(jié)點關(guān)系式可簡記為：

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

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

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

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

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

1.2.3 張力系數(shù)選擇

在地形圖上采用等高線描繪地形時，特征點之間的曲線取最短為優(yōu)，為避免相鄰光滑曲線相交且要通過等值點，采用規(guī)范化的張力系數(shù)σ作為標(biāo)準(zhǔn)值，即

σ為規(guī)范化張力系數(shù)，需要預(yù)先進(jìn)行試驗，取最適合當(dāng)前地圖的張力系數(shù)，當(dāng)發(fā)生變化時，則由標(biāo)準(zhǔn)值重新計算張力系數(shù)[8-9]。

1.3 Bezier函數(shù)法

1.3.1 Bezier曲線插值

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

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

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

1.3.2 切線抹角法

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

圖1 切線抹角法原理圖

為避免出現(xiàn)多余拐點，切線P1P2應(yīng)垂直于∠A的平分線，設(shè)∠A的平分線與BC相交于D，，，進(jìn)而求得切線的斜率為：可計算得到

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

圖2 切線端點計算

設(shè)dt已知，則點E坐標(biāo)的計算公式為：

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

2 算法編程實現(xiàn)與結(jié)果分析

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

2.1 五點法測試

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

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

表1 五點法測試結(jié)果

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

2.2 張力樣條函數(shù)法測試

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

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

2.3 改進(jìn)后的Bezier函數(shù)法測試

采用改進(jìn)后的Bezier函數(shù)法生成光滑曲線需分段進(jìn)行，每次取相鄰3個點坐標(biāo)，首先計算切線點P1、P2的坐標(biāo)，然后利用3個等值點和2個切線點生成Bezier曲線。若生成曲線效果欠優(yōu)，可調(diào)整相關(guān)系數(shù)dt后重新進(jìn)行等值線光滑。

與張力樣條函數(shù)法類似，改進(jìn)后的Bezier函數(shù)法利用dt來控制曲線與原始折線的偏離程度，當(dāng)dt趨近于1時，曲線趨向原始折線；當(dāng)dt較小時，曲線會偏離原始折線。因此，當(dāng)?shù)雀呔€較密集時，應(yīng)適當(dāng)增大dt以避免等值線相交；當(dāng)?shù)雀呔€較稀疏時，應(yīng)適當(dāng)減小dt使地形特征更加明顯。與張力樣條函數(shù)法類似，改進(jìn)后的Bezier函數(shù)法也可通過調(diào)節(jié)系數(shù)值來使最終效果達(dá)到最優(yōu)。同時，相較于張力樣條函數(shù)法，Bezier曲線不會產(chǎn)生多余拐點，光滑效果更優(yōu)。

3 算法特點和適用性

五點法的計算速度快，但可能出現(xiàn)等值線相交等情況，適用于原始數(shù)據(jù)質(zhì)量高、數(shù)據(jù)分布平緩、等值線稀疏、結(jié)果要求不嚴(yán)格的場合，操作簡單，可一鍵生成；張力樣條函數(shù)法的等值線嚴(yán)格正確，但需手動調(diào)節(jié)張力系數(shù)，適用于結(jié)果要求嚴(yán)格、結(jié)果需逐步調(diào)整以達(dá)到最優(yōu)狀態(tài)的場合；改進(jìn)后的Bezier函數(shù)法的等值線嚴(yán)格正確、不產(chǎn)生多余拐點，但需手動調(diào)節(jié)相關(guān)參數(shù)，適用于結(jié)果要求嚴(yán)格、結(jié)果需逐步調(diào)整以達(dá)到最優(yōu)狀態(tài)、不產(chǎn)生多余拐點的場合。