Douglas-Peukcer抽稀算法在水文中的應(yīng)用及編程實(shí)現(xiàn)與優(yōu)化

代 斌,施 藝

(四川省成都水文水資源勘測(cè)中心,成都 611130)

0 引言

水文學(xué)是研究水資源、水循環(huán)和水環(huán)境的科學(xué),涉及大量的水文數(shù)據(jù)。隨著信息技術(shù)的快速發(fā)展和數(shù)據(jù)采集技術(shù)的成熟,水文學(xué)領(lǐng)域也迎來了大數(shù)據(jù)時(shí)代。在目前的大數(shù)據(jù)環(huán)境下,水文數(shù)據(jù)的數(shù)據(jù)量呈現(xiàn)爆炸式增長(zhǎng),傳統(tǒng)的數(shù)據(jù)處理方式已經(jīng)無法滿足需求,數(shù)據(jù)抽稀成為了處理數(shù)據(jù)的一種關(guān)鍵技術(shù),也是現(xiàn)在水文數(shù)據(jù)管理和分析的必要手段。

水文數(shù)據(jù)抽稀是指在保持水文數(shù)據(jù)重要特征的同時(shí),從原始水文數(shù)據(jù)中有選擇地減少數(shù)據(jù)點(diǎn)的數(shù)量,以達(dá)到減少數(shù)據(jù)量、提高計(jì)算效率、降低存儲(chǔ)成本和提高數(shù)據(jù)質(zhì)量的目的。

1 抽稀算法的選擇

1.1 常用算法

目前來說抽稀較為常用的算法有:等步長(zhǎng)法、線段過濾法、垂距限值法及Douglas-Peuker算法。

水文數(shù)據(jù)是按時(shí)間順序排列的數(shù)據(jù)點(diǎn)集合。等步長(zhǎng)法是按照一定的間隔距離對(duì)點(diǎn)集進(jìn)行抽取,其余點(diǎn)全部壓縮掉,在相鄰抽取點(diǎn)間用直線連續(xù);線段過濾法是從點(diǎn)集的第一點(diǎn)開始,若連續(xù)兩個(gè)點(diǎn)組成的線段長(zhǎng)度比給定的閾值小時(shí),刪除這兩點(diǎn),并將這兩點(diǎn)的中點(diǎn)插入兩點(diǎn)原來的位置;垂距限值法是從點(diǎn)集的第一點(diǎn)開始,根據(jù)中間點(diǎn)到前、后兩個(gè)相鄰點(diǎn)連線的距離,來確定是否保留該點(diǎn)的一種抽稀算法,當(dāng)距離大于給定的閾值時(shí)保留該點(diǎn),否則刪除該點(diǎn);Douglas-Peukcer算法是計(jì)算點(diǎn)集中除頭尾兩個(gè)端點(diǎn)外的每個(gè)點(diǎn)與頭尾兩點(diǎn)直線之間的垂直距離,若這些距離的最大值小于給定的閾值就刪除這一組點(diǎn),若大于給定的閾值,則從最大距離處將點(diǎn)集分為兩部分,繼續(xù)以上步驟,該算法以遞歸方式應(yīng)用于這兩個(gè)部分,直到達(dá)到所需的簡(jiǎn)化級(jí)別[1]。

1.2 算法選擇

將水文數(shù)據(jù)的點(diǎn)集看作是以時(shí)間為x軸、數(shù)據(jù)值為y軸的一條連續(xù)曲線,在以上四種抽稀方法中,等步長(zhǎng)法在這四種方法中顯然是相對(duì)較差的一種方法,因?yàn)樗耆魂P(guān)心原始曲線的數(shù)據(jù)是什么,只是按照給定的步長(zhǎng)去抽取數(shù)據(jù)。若步長(zhǎng)較小,抽稀效果不明顯;若步長(zhǎng)較大,則原始曲線嚴(yán)重失真[1]。

另外三種方法與等步長(zhǎng)法有明顯的不同,因?yàn)樗鼈兊某橄《紩?huì)關(guān)系原始曲線中數(shù)據(jù)之間的聯(lián)系。而從它們觀察的對(duì)象來看,線段過濾法觀察的是兩個(gè)點(diǎn)之間的距離,垂距限值法觀察的是三個(gè)點(diǎn)之間的關(guān)系,而Douglas-Peukcer算法觀察的是整條曲線中各個(gè)點(diǎn)的關(guān)系。從設(shè)計(jì)思路上來看,Douglas-Peukcer算法的抽稀效果是相對(duì)較好的,線段過濾法的抽稀效果是相對(duì)較差的[1]。

垂距法較等步長(zhǎng)法與線段過濾法精度更高,同樣具有算法簡(jiǎn)單和易于實(shí)現(xiàn)等優(yōu)點(diǎn),但它是一個(gè)局部的算法,無法從整體上考慮線要素的變化形態(tài)。因此,當(dāng)各相鄰點(diǎn)起伏變化不大時(shí),即使給定的閾值較小,也可能會(huì)刪除掉構(gòu)成大彎曲形態(tài)上的所有中間點(diǎn),造成線要素形態(tài)失真[2]。

Douglas-Peukcer算法的優(yōu)點(diǎn)是它是一個(gè)整體算法,能做到在刪除與保留之間達(dá)到較好的平衡,即能充分減少點(diǎn)的數(shù)量,又能盡量保留特征點(diǎn)[2]。但Douglas-Peukcer算法較其余算法更復(fù)雜,實(shí)現(xiàn)時(shí)需要采用遞歸方法,有一定的難度。本文選取Douglas-Peukcer算法進(jìn)行水文數(shù)據(jù)的抽稀,并對(duì)其在編程中的實(shí)現(xiàn)進(jìn)行了分析探討。

2 在編程中的實(shí)現(xiàn)與優(yōu)化

2.1 算法實(shí)現(xiàn)

目前遙測(cè)技術(shù)已經(jīng)得到廣泛應(yīng)用和普及,動(dòng)輒5 min甚至1 min一個(gè)的數(shù)據(jù)雖然可以幫助我們更及時(shí)地掌握水文信息以及更精細(xì)化地進(jìn)行水資源管理,但在數(shù)據(jù)存儲(chǔ)的成本、數(shù)據(jù)的分析整編、計(jì)算所需的時(shí)間和資源、數(shù)據(jù)可視化效果上都帶來了不少的麻煩,需要先對(duì)數(shù)據(jù)進(jìn)行抽稀處理后再進(jìn)行應(yīng)用。因此,本文以遙測(cè)水位數(shù)據(jù)為例,對(duì)Douglas-Peukcer抽稀算法的實(shí)現(xiàn)進(jìn)行分析探討。

為方便普通職工在EXCEL中調(diào)用該過程,本文采用VBA語言進(jìn)行編譯,實(shí)現(xiàn)代碼如下:

Function Compr_DP(Arr_COMPR, Arr_TEMP, Start_Index, End_Index)

If Start_Index

Dim Start_X#, Start_Y@, End_X#, End_Y@

Start_X = Arr_COMPR(Start_Index, 1): Start_Y = Arr_COMPR(Start_Index, 2)

End_X = Arr_COMPR(End_Index, 1): End_Y = Arr_COMPR(End_Index, 2)

Dim i&, Index&, Dist@, Max_Dist@, Point_X#, Point_Y@

For i = Start_Index + 1 To End_Index - 1

Point_X = Arr_COMPR(i, 1): Point_Y = Arr_COMPR(i, 2)

Dist = Calc_DIST (Start_X, Start_Y, End_X, End_Y, Point_X, Point_Y)

If Dist >Max_Dist Then Index = i: Max_Dist = Dist

Next i

If Max_Dist >3 Then

Arr_TEMP(Index, 1) = Arr_COMPR(Index, 1)

Arr_TEMP(Index, 2) = Arr_COMPR(Index, 2)

Call Compr_DP(Arr_COMPR, Arr_TEMP, Start_Index, Index)

Call Compr_DP(Arr_COMPR, Arr_TEMP, Index, End_Index)

End If

End If

End Function

Function Calc_DIST#(Start_X, Start_Y, End_X, End_Y, Point_X, Point_Y)

Dim a, b, c, p, S

a = Abs(Sqr((Abs(Start_X - End_X) ^ 2) + (Abs(Start_Y - End_Y) ^ 2)))

b = Abs(Sqr((Abs(Start_X - Point_X) ^ 2) + (Abs(Start_Y - Point_Y) ^ 2)))

c = Abs(Sqr((Abs(Point_X - End_X) ^ 2) + (Abs(Point_Y - End_Y) ^ 2)))

p = (a + b + c) / 2

S = Sqr(Abs(p * (p - a) * (p - b) * (p - c)))

Calc_DIST = S * 2 / a

End Function

其中Function Compr_DP是主過程,閾值設(shè)定為3(cm),根據(jù)曲線上所有點(diǎn)的坐標(biāo)依次計(jì)算各點(diǎn)到Start、End兩點(diǎn)所在直線的垂直距離(通過調(diào)用Function Calc_DIST),得到最大距離Max_Dist,若Max_Dist小于閾值3,則將這條曲線上的中間點(diǎn)全部舍去,以直線段作為曲線的近似,該段曲線處理完畢,若Max_Dist小于閾值3,保留Max_Dist對(duì)應(yīng)的坐標(biāo)點(diǎn),并以該點(diǎn)為界,把曲線分為兩部分,對(duì)這兩部分重復(fù)使用該方法,直到所有Max_Dist均小于閾值3,即完成對(duì)曲線的抽稀;Function Calc_DIST是子過程,通過3個(gè)點(diǎn)的坐標(biāo)求得各點(diǎn)到Start、End兩點(diǎn)所在直線的垂直距離。

圖1 Douglas-Peukcer算法示意

2.2 算法優(yōu)化

上一章節(jié)中,我們計(jì)算得到的距離是各目標(biāo)點(diǎn)到Start、End兩點(diǎn)所在直線的垂直距離,即目標(biāo)點(diǎn)與Start、End兩點(diǎn)形成的三角形的高,而與大部分其他類型的數(shù)據(jù)集是以空間為坐標(biāo)進(jìn)行分布不同,水文數(shù)據(jù)基本上都是以時(shí)間為序列進(jìn)行分布,我們可以將各目標(biāo)點(diǎn)到Start、End兩點(diǎn)所在直線在同一坐標(biāo)軸(即以同一時(shí)間為參考)的距離代替垂直距離,這樣既能減輕計(jì)算量(將三角運(yùn)算變?yōu)榧訙p乘除),又更符合水文數(shù)據(jù)以時(shí)間為序列進(jìn)行分布的特征。由于在VBA中默認(rèn)的時(shí)間運(yùn)算以d為單位,而水文計(jì)算中水位默認(rèn)以m為單位,我們以更精細(xì)的刻度來保證計(jì)算的準(zhǔn)確性,以 min為時(shí)間單位、以cm為水位單位。通過修改子過程Function Calc_DIST的代碼完成以上內(nèi)容,實(shí)現(xiàn)代碼如下:

Function Calc_DIST@(Start_X, Start_Y, End_X, End_Y, Point_X, Point_Y)

Dim a1&, a2&, b1&, b2&

a1 = (End_X - Start_X) * 1440

If a1 = 0 Then Calc_DIST = 0: Exit Function

a2 = (Point_X - Start_X) * 1440

b1 = (End_Y - Start_Y) * 100

b2 = (Point_Y - Start_Y) * 100

Calc_DIST = Abs(b2 - b1 * a2 / a1)

End Function

圖2 Douglas-Peukcer算法(改進(jìn)后)示意

2.3 實(shí)例計(jì)算與分析

為驗(yàn)證以上代碼的實(shí)際運(yùn)行效果,本次分析以石堤堰水文站實(shí)測(cè)水位資料為實(shí)例樣本,對(duì)Douglas-Peukcer算法改進(jìn)前后的效果,以及跟四川省水旱災(zāi)害防御決策支持系統(tǒng)的抽稀處理進(jìn)行對(duì)比(閾值均設(shè)置為3 cm)。

抽稀算法的優(yōu)劣,一般通過數(shù)據(jù)處理時(shí)間、數(shù)據(jù)抽稀率、數(shù)據(jù)偏差來評(píng)價(jià)。數(shù)據(jù)處理的速度(時(shí)間復(fù)雜度)主要以運(yùn)行時(shí)間來評(píng)價(jià);數(shù)據(jù)抽稀率(空間復(fù)雜度)是抽稀掉的數(shù)據(jù)點(diǎn)數(shù)與原始數(shù)據(jù)點(diǎn)數(shù)之比。數(shù)據(jù)偏差依據(jù)《水文資料整編規(guī)范》(SL/T 247-2020)要求,處理后的水位變化過程完整,經(jīng)處理以后計(jì)算的日平均水位與采用所有數(shù)據(jù)計(jì)算的日平均水位相差不宜超過2 cm[3]。我們分別對(duì)抽稀前后的水位數(shù)據(jù)進(jìn)行整編,通過對(duì)比抽稀前后得到的日平均水位及月最大、最小、平均水位進(jìn)行數(shù)據(jù)偏差評(píng)價(jià)。

對(duì)石堤堰水文站2022年全年共98 488個(gè)水位數(shù)據(jù)進(jìn)行抽稀處理(見表1),四川省水旱災(zāi)害防御決策支持系統(tǒng)的抽稀用時(shí)14.9 s,抽稀數(shù)據(jù)81927條,抽稀率83.2 %;Douglas-Peukcer算法(改進(jìn)前)用時(shí)12.7 s,抽稀數(shù)據(jù)87108條,抽稀率88.4 %;Douglas-Peukcer算法(改進(jìn)后)用時(shí)3.47 s,抽稀數(shù)據(jù)86707條,抽稀率88.0 %。通過對(duì)比,四川省水旱災(zāi)害防御決策支持系統(tǒng)抽稀用時(shí)最長(zhǎng),抽稀率最低,改進(jìn)前的Douglas-Peukcer算法抽稀率最高,但數(shù)據(jù)處理的用時(shí)較長(zhǎng),改進(jìn)后的Douglas-Peukcer算法數(shù)據(jù)處理的速度明顯更快,且抽稀率與改進(jìn)前相差不大,在時(shí)間復(fù)雜度與空間復(fù)雜度上均能達(dá)到比較滿意的效果,效率最佳。

表1 抽稀算法實(shí)際運(yùn)行效率比較

由于全年遙測(cè)數(shù)據(jù)量較為龐大,對(duì)未抽稀的原始數(shù)據(jù)進(jìn)行整編難度較大,因此本文采用單月數(shù)據(jù)進(jìn)行數(shù)據(jù)偏差比較,以石堤堰水文站2023年5月水位資料為樣本,各算法與未抽稀前的整編成果對(duì)比見表2。

表2 抽稀算法實(shí)測(cè)數(shù)據(jù)偏差比較

通過與原始數(shù)據(jù)的整編成果進(jìn)行對(duì)比,各算法的成果均能滿足《水文資料整編規(guī)范》(SL/T 247-2020)要求,Douglas-Peukcer算法較四川省水旱災(zāi)害防御決策支持系統(tǒng)稍好,算法改進(jìn)后效果與改進(jìn)前基本一致。

將抽稀前后的水位過程線進(jìn)行對(duì)比,選取5月1日8時(shí)至5月2日8時(shí)的時(shí)間段,可以看出,在大部分時(shí)間,Douglas-Peukcer算法的水位變化過程都更貼合原始數(shù)據(jù),算法改進(jìn)后效果與改進(jìn)前基本一致。

圖3 水位過程線對(duì)比

3 結(jié)語

本文對(duì)Douglas-Peukcer抽稀算法及其與其他抽稀算法的優(yōu)劣進(jìn)行了簡(jiǎn)要的分析,通過編程實(shí)現(xiàn)了Douglas-Peukcer抽稀算法在水文中的應(yīng)用,并根據(jù)水文數(shù)據(jù)的特性對(duì)算法進(jìn)行了優(yōu)化。本文根據(jù)石堤堰水文站實(shí)測(cè)資料為實(shí)例樣本對(duì)算法進(jìn)行了運(yùn)行效率測(cè)試與數(shù)據(jù)偏差分析,Douglas-Peukcer抽稀算法在各方面均領(lǐng)先于四川省水旱災(zāi)害防御決策支持系統(tǒng)的抽稀處理算法,改進(jìn)后的Douglas-Peukcer抽稀算法在處理效果上與改進(jìn)前基本一致,但在運(yùn)行效率上較改進(jìn)前有較大提升,實(shí)現(xiàn)了高效率、高抽稀率與保真效果之間的平衡。