顧及面積差的層內等效二段線替代法優化聲速剖面

肖元弼,彭認燦,董 箭,馬政偉,劉 聚1,

1. 海軍大連艦艇學院作戰軟件與仿真研究所, 遼寧 大連 116018; 2. 海軍大連艦艇學院戰術學博士后科研流動站, 遼寧 大連 116018; 3. 海軍大連艦艇學院軍事海洋與測繪系, 遼寧 大連 116018; 4. 91388部隊, 廣東 湛江 524000

對于多波束勘測這種靠聲波在水中傳播進行測距的測深系統,聲速作為最重要的計算參數,其精度的高低將直接影響測深結果。聲線跟蹤法在追蹤波束腳印時需要用到聲速剖面(sound velocity profile,SVP),而高采樣密度的聲速剖面數據會導致聲線跟蹤法在分層計算和逐層累加過程中耗時較長,降低計算效率。通過簡化聲速剖面來減少計算層數是提高聲線跟蹤計算效率的一種最直接方法,但盲目的簡化會忽略聲速剖面一些特征信息,從而降低測深精度。再加之,多波束測深系統測深覆蓋面廣、數據量大的特點,因此,需要探究一種在聲速剖面簡化過程中能兼顧數據精度的方法[1-16]。

道格拉斯-普克法(Douglas-Peucker,D-P)是一種在制圖系統中常用的地形特征線簡化算法,主要通過觀測點到曲線首尾連線的垂直距離作為閾值參考指標,利用幾何相似性來簡化復雜曲線。若應用于聲速剖面簡化,其存在的問題是聲速剖面橫縱坐標分別是速度與深度,其單位不統一將導致幾何垂直距離及閾值沒有實際物理意義。為解決聲速剖面簡化問題,國內外學者進行了大量研究。文獻[17]提出了一種基于聲速最大偏移量(maximum offset of sound velocity,MOV)的D-P改進算法。將聲速維度上最大距離作為閾值的參考指標,對采樣點進行取舍,相比較D-P算法,MOV算法在聲速剖面應用下更為合理,并且能夠保留D-P算法對于曲線幾何特征保留性強的優勢。但該算法以整體水深誤差1%作為極限指標,來尋找最優閾值,存在分層內聲速剖面等效性不完全合理的情況,尤其在特殊的聲速剖面情況下閾值收斂性較差。文獻[18]提出了一種基于面積差約束的聲速剖面自適應簡化方法,通過對所有采樣點進行逐點取舍判斷,在確保聲速剖面精度和基本空間變化結構情況下,簡化率和跟蹤精度較高,在聲速剖面起伏波動頻率較高的區域尤為明顯。但整體聲速剖面運用面積差約束尋找等效替代,會導致層內聲速剖面幾何特征保持性較差,并且整體計算過程較為復雜。

本文針對聲速剖面簡化后層間精度不高的情況,提出一種聲速剖面自適應分層優化方法——層內等效二段線替代法(in-layer equivalent two-segment line substitution method,IETL)。該方法在實測SVP數據的基礎上,對MOV算法進行改良,利用等效面積差作為約束條件,通過構建層內等效二段線,來對層間聲速剖面進行優化,以便提高層間聲速剖面精度。試驗表明,經過本文方法優化后的聲速剖面,在保證整體數據精度的基礎上,能夠進一步提高層間數據精度,尤其在地形波動較大的地區,其潛在的工程應用優勢尤為明顯。

1 方法與簡化模塊

1.1 顧及面積差的IETL法

1.1.1 聲速剖面簡化原則調整

聲速剖面簡化的主要目的是通過減少聲速剖面上的聲速點個數來減少聲速的分層數,從而提高測深點位置計算效率。但不同于簡單的幾何優化,聲速剖面簡化原則最重要的是要保證測深點位置精度,并在此基礎上提高聲線跟蹤效率。因此,傳統聲速剖面簡化原則可歸納為:①保證波束在海底總體測深數據精度;②聲速變化相對穩定的區域可合并為一個水層;③不破壞聲速剖面結構特征信息,保留特征性較強的采樣點。

但傳統簡化原則存在幾個問題,首先原則①由于只考慮整體聲速替代剖面簡化,而忽略了層間聲速剖面等效情況,對于海底地形變化幅度較大的區域(如斷崖等),勢必會影響測深精度;原則③對于聲速剖面變化復雜的水層,是否可以在保證測深精度的情況下,對特征采樣點進行改變,以達到簡化目的。本文針對問題情況,將簡化原則調整為:①保證波束在海底總體及每個水層的測深數據精度;②聲速變化相對穩定的區域可合并為一個水層;③不破壞聲速剖面結構特征信息,保留或適當改變特征性較強的采樣點。

1.1.2 MOV算法

D-P算法是一種簡化多段線的常用方法[19],如圖1(a)所示。其基本思想是先確定閾值Q,對復雜線段的首尾端點進行連線,再依次計算出原始折線段間其余各點P到該線段的距離Hi,找到其中最大的垂直距離Himax與閾值Q進行比較,計算公式為[19]

(1)

圖1 D-P算法與MOV算法Fig.1 D-P algorithm and MOV algorithm

式中,x表示聲速;y表示深度;P1表示原始聲速剖面第一個點;Pn表示原始聲速剖面最后一個點;Pi表示替代聲速剖面點。

若Himax≤Q,則除首尾兩點,其余各點都刪除,算法結束;若Himax>Q,則保留該點Pi,并從該點處分別向首尾兩端P1、Pn進行連線,形成兩條新的線段P1Pi、PiPn。重復迭代上述步驟,找到各自線段的所有點到新形成線段的最大垂直距離,繼續與閾值進行比較,直到所有點到各自線段垂直距離都小于閾值,即完成簡化操作。但由于聲速剖面縱坐標為深度(單位為m),橫坐標為聲速(單位為m/s),兩者單位不統一使得上述方法的幾何連線垂直距離無實際物理意義,導致D-P算法無法應用于SVP簡化。

針對D-P算法在SVP方面的應用問題,文獻[17]在D-P算法基礎上,提出了改進的MOV算法,由點到線段的垂直距離改為聲速的最大偏移量,同閾值進行比較。如圖1(b)所示,對首尾端點進行連線,依次計算出原始折線段間其余各點到該線段的橫坐標偏移量,即聲速偏移量Ci,找到最大距離與閾值進行比較,計算公式為[20]

(2)

而對于點的判斷保留步驟與D-P算法相同。因此,MOV算法能夠很好地解決D-P算法在簡化復雜聲速剖面時,閾值與垂直距離沒有實際物理意義的問題,使得此簡化模型能夠應用于聲速剖面。

等效聲速剖面表示若采用某一曲線作為簡化聲速剖面進行計算得到的測深點位與原始聲速剖面得到的測深點位一致,即認為該曲線為等效聲速剖面[20]。雖然MOV算法作為D-P模型的改進算法能夠對曲線進行最優簡化,在設定合適的閾值后,能夠將整體測深數據誤差控制在1%以內,滿足國際海道測量組織(International Hydrographic Organization,IHO)對水深測量的標準,但不同于誤差傳遞,聲線跟蹤法計算深度是一個累加的過程,可能會存在某一層聲速整體偏低,另外一層又整體偏高,與替代聲速剖面間偏移量抵消導致測深結果影響不大的情況。因此,整體測深數據誤差低并不能說明聲速剖面每一層測深誤差同樣低。對聲速剖面而言,每一層的簡化曲線并不一定等效于原始聲速剖面,在地形變化較大的區域,尤其是聲速變化較大的淺水區,簡化曲線會對測深數據產生較大的影響。如圖1(b)所示,明顯其中第三層替代聲速剖面與原始聲速剖面存在面積差,導致該層部分的聲速剖面等效性較差?？梢?雖然整體測深數據滿足精度要求,但每一層的替代聲速剖面并不一定是等效的聲速剖面,因此需要對該層進行等效替代。

1.1.3 最大等效面積差提取方法

圖2為某一水層聲速剖面與等效聲速剖面面積差的仿真關系,圖中紅虛線為簡化聲速剖面,實線表示原始聲速剖面及等效聲速剖面面積差的曲線。由圖2可以看出,除了水層的兩端點外,極值也產生在a、b、c3點,這3點都對應的是替代聲速剖面與原始聲速剖面的交點處,因此尋找替代聲速剖面與原始聲速剖面面積差的最值,只需在端點及兩條線交點處比對即可。記錄下水層中聲速剖面面積差的最大值,若最大值符合精度要求,即整個水層的替代聲速剖面都符合精度要求,如圖2中c點為該水層等效聲速剖面面積差的最大值。

圖2 聲速剖面與等效聲速剖面面積差關系Fig.2 Area difference between SVP and equivalent SVP

面積差與聲線跟蹤精度的關系模型如下

(3)

式中,εx、εy分別表示水平位移和深度的相對誤差,由原始聲速剖面聲線跟蹤得到的測深點坐標(x,y)及替代聲速剖面得到的坐標(xt,yt)作參考得到;gi、gt分別表示原始聲速剖面的梯度和等效聲速剖面的梯度;ΔS為面積差;c0為初始聲速;θ0為入射角。式(3)表明在入射角θ0與初始聲速c0已知的情況下,相對誤差εx、εy只與面積差ΔS有關?？赏ㄟ^比較水層內最大面積差ΔSmax得到相對誤差εmax與限差ε0的關系來判斷替代聲速剖面水層是否需要進行進一步優化。若εmax≤ε0,說明層內替代聲速剖面優化程度滿足精度指標,保留層內替代聲速剖面點不做處理;若εmax>ε0,則說明層內替代聲速剖面優化精度不夠,需要進一步優化。

1.1.4 層內聲速剖面等效替代二段線模型建立

已知,i層間原始聲速剖面與y軸(深度軸)圍成的面積是固定的Vi,因此其層內聲速積分也是一個定值

0.5(y-ya)(x+xa)+0.5(ya+Δy-y)·

(x+xa+Δx)=Vi

(4)

由式(4)可知,深度y相同的情況下存在唯一一點(x,y)分別與i層內兩端點(xa,ya)、(xa+Δx,ya+Δy)進行連線得到兩條線段(又稱二段折線段),并與y軸(深度軸)圍成的面積為Vi,即與原始聲速剖面所得的面積差值為0,認為該二段折線段為該層等效聲速剖面。

如圖3(a)所示,對于某一水層,根據MOV算法原始聲速剖面的聲速偏移量都沒有達到閾值k的界限,因此不需要進一步分層,但由圖3中可明顯看出原始聲速剖面與替代聲速剖面EG不等效。而對于聲速最大偏移量F處所在的深度m線上,存在一點F2分別到E點與G點所形成的線段,與y軸圍成的面積與該層原始聲速剖面圍成的面積等效(圖3(b))。

圖3 顧及面積差的等效二段線替代法Fig.3 In-layer equivalent two-segment line substitution method considering area difference

式(4)轉化為式(5)

(5)

由式(5)可得存在一條直線,該直線上任意一點到兩端點(xa,ya)、(xa+Δx,ya+Δy)的連線所形成二段折線段為等效聲速剖面,并且該直線的斜率與兩端連線的直線斜率相同,如圖3(c)中l線所示。

在等效聲速剖面點集合所形成的直線l上,存在一點(xt,yt)使得等效聲速剖面與原始聲速剖面間偏移量差的平方和為最小值(也可理解為凈偏移量最小值,為了方便計算采用忽略正負號影響,這里采用平方和),即該點為最優解,是所有替代的等效聲速剖面中最接近原始聲速剖面的,如圖3(c)中Fa所示。為方便尋找最優解,式(6)構造拉格朗日函數求條件極值

F(xt,yt)=f(xt,yt)+λφ(xt,yt)

(6)

式(7)為約束方程由式(4)轉化而來

φ(xt,yt)=0.5(yt-ya)(xt+xa)+0.5(ya+
Δy-yt)(xt+xa+Δx)-Vi

(7)

式(8)—式(9)分別表示以最優解(xt,yt)作為新增點,分別與上下兩端點連線所形成的兩部分替代聲速剖面表達式

(8)

(9)

式中,x=g(y)表示原始聲速剖面,是個已知的定值。

式(10)為偏移量差的平方和

(10)

將式(7)和式(10)代入式(6),分別對xt、yt、λ求偏導(式(11)—式(13))

(11)

(12)

(13)

可得到最優解(xt,yt),即等效面積差為0作約束條件的最優替代聲速剖面二段線拐點。分別計算測點與等效二段線做偏移量差的平方和,計算最小值。

1.2 顧及面積差的等效二段線替代法優化流程

1.2.1IETL方法

直接在層內增加等效聲速剖面二段折線,雖然會提高層內數據精度,但同時也會增加整個聲速剖面的水層數,降低簡化率。若在水層數保持不變的情況下,盡可能對已簡化的采樣點進行優化,即可在保證簡化率的情況下提高替代聲速剖面。先從等效聲速剖面面積差最大的水層進行分析,如圖4(a)所示,該圖為經過MOV優化的聲速剖面,其中i層等效聲速剖面面積差最大,在不增加水層的前提下先分析是否能夠對i層兩端點Di、Di+1進行優化。對點Di進行分析,把Di-1DiDi+1所在的i-1和i層合并為一層,按照1.1.4節算法,得到最優點DDi代替原Di作為新聲速剖面簡化點,并記錄下DDi與Di+1形成的新水層內等效聲速面積差(圖4(b));同理,對原i層另外一端點Di+1的分析,把DiDi+1Di+2所在的i和i+1層合并為一層,得到最優點DDi+1,記錄下DDi+1與Di形成的新水層內等效聲速面積差(圖4(c))。對比兩個等效聲速面積差,若最小值滿足精度要求,則將該替代點作為新的保留,繼續優化。否則說明該層兩端點難以進行優化,需增加新的等效聲速剖面點。

圖4 聲速剖面層內端點優化替代方法Fig.4 An alternative approach to SVP in-layer endpoint optimization

1.2.2 方法簡化模塊

根據前文的問題分析,本文方法是顧及面積差的等效二段線聲速剖面替代優化法,可歸納為兩大模塊分別為:MOV算法簡化模塊和替代優化模塊,其中替代優化模塊又可分為整體優化模塊及層內端點優化模塊,如圖5所示。

圖5 優化模塊流程Fig.5 Flowchart of optimization module

MOV算法簡化模塊用于聲速剖面的簡化,其核心是1.1.2節所述的MOV算法,主要實現步驟如下。

(1) 數據輸入。輸入原始聲速剖面點集,設定閾值參數T。

(2) 提取聲速最大偏移。計算各點的偏移量,并從中提取最大偏移量Cmax。

(3) 剖面線段切割。判斷Cmax和T的大小,若Cmax≤T,轉至步驟(4);若Cmax>T,保存該點,并將剖面線段在該點處分為兩段,并返回步驟(2)。

(4) 數據輸出。將原始聲速剖面的首尾端點及所有后形成的點記錄下,輸出為替代聲速剖面。

整體優化模塊是在MOV算法簡化的替代聲速剖面基礎上,對現有的特征點進行修正或者增加,以便提高聲速剖面精度。

(1) 數據輸入。輸入MOV算法簡化模塊所得的替代聲速剖面點集及原始聲速剖面點集。

(2) 提取最大等效面積差的水層。利用1.1.3節方法找到等效面積差最大的所在層,并對該層進行優化分析。

(3) 精度評估A。利用1.1.3節算法,通過最大面積差ΔSmax得到相對誤差εmax,若εmax≤ε0,轉至步驟(6);若εmax>ε0,則轉至步驟(4)。

(4) 對層內端點進行優化。對最大面積差所在的i層端點進行優化,詳細情況在層內端點優化模塊中說明。

(5) 精度評估B。評估端點進行優化后的i層,利用1.1.3節算法,若εmax≤ε0,則接受新端點替代原簡化點,轉至步驟(2);若εmax>ε0,則拒絕新端點,并在i層內添加新的等效點,形成新的替代聲速剖面,轉至步驟(2)。

(6) 輸出數據。輸出新的聲速替代剖面。

層內端點優化模塊的步驟如下。

(1) 數據輸入。輸入最大等效面積差所在i層的替代聲速剖面點集和原始聲速剖面點集。

(2) 對i層內兩端點進行替代分析。分析上端點,將i-1層與i層合并為一層,建立等效面積線;分析下端點,將i層與i+1層合并為一層,建立等效面積線。

(3) 得到最接近原始聲速剖面的點。在替代聲速點所得等效二段線使新替代聲速剖面與原始聲速剖面間的凈面積差為最小值時,計算新的i層上端點等效面積差ΔSCA和下端點的等效面積差ΔSCB。

(4) 保留最優替代端點。提取兩端點所得的數據進行判斷,若ΔSCA≤ΔSCB,用新替代的上端點DDi替代原始替代聲速剖面上端點Di;若SCA>SCB,用新替代的下端點DDi+1替代原始替代聲速剖面下端點Di+1。

(5) 輸出數據。輸出優化后新i層的聲速替代剖面。

2 試驗及分析

2.1 閾值的選取

閾值是簡化模型中重要的參數指標,為確定閾值的影響,對一段聲速剖面分別設定閾值為0.1、0.2、0.3、0.5、1、2 m/s,其替代優化分層效果如圖6所示,圖7記錄不同閾值影響下的精簡率和運算時間關系曲線。由圖6、圖7可以看出,當閾值設定為0.1和0.2 m/s時,兩種方法的精簡率是相同的,但本文方法會對一些聲速剖面簡化特征點進行改動;當閾值設定為0.3和0.5 m/s時,本文方法相較于MOV算法不僅對一些特征點進行改動,還增加了一些特征點;當閾值設定為1和2 m/s時,MOV算法對聲速剖面的精簡替代表現較差,雖然精簡率較高但是不能很好地表現出聲速剖面的特征,本文方法在這方面表現較好。

圖6 設定不同閾值的簡化聲速剖面效果圖對比Fig.6 Comparison of SVP renderings with different thresholds

圖7 不同閾值影響下的精簡率和運算時間關系曲線Fig.7 Relationship curves of reduction rate and operation time under different thresholds

由圖7可以明顯看出,本文方法相較于MOV算法計算時間較長,這是由于兩種方法從計算模型角度出發,MOV算法只需提取每個采樣點最大偏移量進行分析,而本文方法中多重循環嵌套明顯在計算量方面要比MOV算法大得多,因此在前期簡化過程中盡量使用MOV算法。MOV算法設定的閾值大小會直接影響整體改進方法的計算效率,閾值過小會導致簡化率降低,使得多波束測深系統計算深度時效率降低;閾值過大會使得MOV算法過早停止簡化,增加改進方法簡化計算,降低整體簡化效率。為在分層精度為主的基礎上兼顧計算效率,閾值范圍選取在0.2～0.4 m/s這個范圍最優,根據聲速剖面波動程度可做適當調整。

2.2 實測數據應用分析

為驗證本文方法的有效性,通過對東海某測量船作業的實測數據進行驗證分析。首先提取該區域的聲速剖面數據,并對其進行不同方法的簡化,如圖8所示。設定閾值為0.2 m/s,本文IETL方法能夠將300個聲速剖面采樣點簡化為15個,MOV算法能將其簡化為13個。在簡化率上MOV算法略勝一籌,然而由圖8可以看出,其中一部分原始聲速剖面和MOV算法簡化的替代聲速剖面間存在較明顯的面積差。圖8中A點為原始聲速剖面和替代聲速剖面的交點,根據1.1.3節推導,該點處二者面積差最大,并且使用MOV算法簡化的聲速剖面計算的水深數據與原始測深數據誤差也較大。

圖8 實測數據聲速剖面不同簡化效果對比Fig.8 Comparison of different simplification effects of measured data SVP

圖9為使用不同替代聲速剖面計算后的海底地形圖,其中圖9(a)、(b)分別為原始海底地形圖和使用MOV算法替代聲速剖面計算后的海底地形圖,對比兩圖可以發現,在圖9(b)中H、J、K3處深度差較為明顯,其原因為在該深度附近的替代聲速剖面等效面積差較大,導致計算后的深度存在誤差。圖9(c)為運用本文方法的替代聲速剖面計算后的海底深度,與原始海底地形深度差較小。對比海底地形直觀圖不明顯,為便于對比分析,對兩種替代方法的誤差進行統計如表1所示,圖10為深度差與深度關系的散點分布圖。由圖10可知,MOV算法在深度-55 m附近與原始聲速剖面間的深度差比較大,隨著深度的增加,替代聲速剖面與原始聲速剖面間面積差逐漸縮小,深度差也隨之縮小。本文IETL方法在該深度區間內,與原始聲速剖面間面積差較小,因此應用效果要優于MOV算法。由表1可知,MOV算法雖然平均相對誤差在1%以內,但是部分數據誤差大于1%,而IETL方法能夠將整體和部分測深數據誤差都能控制在1%以內,滿足國際海道測量組織對水深測量的標準,并且精度更高。

表1 不同SVP簡化方法的誤差統計

圖9 應用不同簡化聲速剖面的實測海底地形Fig.9 Measured seabed topography using different simplified SVP

圖10 應用不同簡化聲速剖面的深度差與深度關系散點分布Fig.10 Scatterplot of depth difference and depth relationship with different simplified SVP

3 結 論

本文提出了一種顧及面積差的聲速剖面層內等效二段線替代法。該方法在MOV算法的基礎上,針對其層間替代聲速剖面不等效和精度不高的情況,通過構造最優的等效二段線來優化替代聲速剖面。在原始聲速剖面與替代聲速剖面總體等效保證最低水深精度的同時,會出現由于部分層間聲速剖面等效面積差較大,導致該水層替代聲速剖面計算出深度精度較低。本文提出的IETL方法,適用于地形起伏較大的海底,此法在有效保留聲速剖面彎曲特征點的同時,明顯提高了層間聲速數據的精度。

本文方法也存在一定問題,例如,為保證層間數據精度而在客觀上犧牲了少量的精簡率,并且該方法模型的計算要比MOV算法復雜,聲速剖面簡化提取特征點的計算時間也較長。因此,在地形起伏不大的區域,其應用優勢不明顯,并且在后續工作中需要對地形不同起伏程度的適應性,以及不同聲速剖面等效面積差的影響效果進行進一步分析。