基于Python的多波束測深數據壓縮研究

張明希，許捍衛,于艷超，俞禮彬

（1. 河海大學 地球科學與工程學院，江蘇 南京 210098）

基于Python的多波束測深數據壓縮研究

張明希1，許捍衛1,于艷超1，俞禮彬1

（1. 河海大學 地球科學與工程學院，江蘇 南京 210098）

在曲線數據壓縮的垂距限差法基礎上，引入總體最小二乘算法對多波束測深的ping條帶數據進行壓縮，利用Python工具庫實現總體最小二乘壓縮算法并集成到ArcToolbox中。與傳統道格拉斯-普克壓縮算法具有接近的壓縮比，但精度更高，能夠在整體上更好地表示原始數據。

多波束數據;抽稀;總體最小二乘；Python

多波束測深系統是一種高效、高精度的水下地形測量系統，數據繁雜，即使單個ping也能產生幾兆數據，對多波束測深數據的壓縮處理非常重要。多波束測深系統沿著測線逐ping觀測，并以ping為順序存儲，可以從Ping數據的壓縮入手。關于數據壓縮問題，國內外學者發展了許多算法。道格拉斯-普克算法（RDP）從整體上對曲線段進行考慮，通過保留特征點、刪除冗余點來達到壓縮目的[1]。總體最小二乘（TLS）方法是近年發展起來的能同時顧及觀測值誤差和模型系數矩陣誤差的數學方法[2]。與RDP不同的是,TLS算法沒有簡單地刪除距離誤差小于限差的點, 而是略微放松對端點的精度要求, 通過TLS作直線擬合,從而提高曲線壓縮的精度，在整體上最優逼近原始曲線。盧銀宏等[3]綜合利用RDP算法和TLS算法來壓縮多波束測深數據，采用Matlab仿真二維ping數據進行壓縮處理，忽略了當曲線的曲率很大時，兩條擬合直線的交叉點可能會遠離原始曲線的情況。本文基于ArcPy站點包和開源Python工具庫，以Python作為腳本語言，開發了多波束測深三維ping數據點抽稀算法，具有易用性和重用性特點。

1 Python空間數據處理特點

Python語言作為穩定的動態腳本語言與ArcGIS高度集成，其靈活性和重用性可以方便實現地理空間數據的處理流程[4]。目前，ArcGIS10．0已集成基于Python的ArcPy站點包，可理解為ArcGIS的Python API。ArcPy提供了大量類和函數，可以輕松執行ArcGIS工具箱中的工具，并創建原生對象，如幾何、柵格、空間參考等[5]。

另一方面，Python具有高級的內建數據結構，優秀的科學計算擴展庫SciPy、NumPy（包括科學計算、統計分析、可視化等模塊），高效的空間分析開源庫GDAL、Shapely等。其中SciPy包含的Python模塊有最優化、線性代數、積分、插值、特殊函數、快速傅里葉變換、信號處理和圖像處理、常微分方程等函數[6]。Shapely庫是針對平面幾何對象的空間分析工具。基于點集理論，其將幾何對象分為點、線、面3種類型，每種數據類型又有3個屬性：內部集、邊界集、外部集。Shapely以Python風格的代碼形式訪問GEOS庫，將PostGIS的一般操作用Python語言加以實現，適合處理一些非常規的地理空間數據問題。

2 基于Python的總體最小二乘點抽稀

2.1 總體最小二乘直線擬合模型

傳統最小二乘法擬合直線時，一般令直線方程為：

其中，x，y是測點坐標；a，b為待估參數。最小二乘法認為測量誤差都存在于y中，而x沒有誤差或在進行平差時不予考慮。則誤差方程可以表示為：

按最小二乘準則，得出其最小二乘解為：

一般認為，觀測值y的獲取方法相同，因此P為單位陣。

最小二乘只認為觀測值y含有測量誤差，事實上x、y均是采用測量方法獲取的，按照總體最小二乘思想，應該對觀測值x和y都進行改正。考慮到觀測值x的誤差，則條件方程為：

寫成矩陣形式為：

其中， Vy為y的改正數；EA為系數陣A的改正矩陣。由于系數陣里含有固定列，運用混合總體最小二乘的迭代解法[7],求出參數a與b。

總體最小二乘直線擬合算法需要在SciPy的線性代數庫linalg基礎上，利用NumPy豐富的數組處理函數，按照總體最小二乘的迭代解法，重新編寫Python腳本計算。

2.2 基于總體最小二乘的多波束數據點抽稀算法

RDP算法在給定限差情況下的逼近精度不高。如圖1，當采用RDP算法對所示曲線數據壓縮時，會得到線ABC。在保留點AB之間，被刪除的點全部位于保留直線的一側；在保留點BC之間，被刪除點也全部位于保留直線的一側。顯然，用ABC表示原始數據具有明顯的偏差，且均方誤差較大。如圖1a，采用基于TLS的分段直線擬合算法，對AB和BC之間的原始點分別作直線擬合，得到兩條直線A′B′、B′C′，它們相交于B′。此時，利用A′-B′-C′代替A-B-C，可以提高對A、B和C點的壓縮精度。但在個別情況下，如圖1b，曲線的曲率很大時，兩條擬合直線的交叉點可能遠離原始曲線，明顯偏離了原曲線的形狀。

圖1 直線擬合和交叉點遠離原始曲線的處理

利用多波束測深的ping點數據，每個ping的相鄰數據點由相鄰波束腳印的中心確定，每個ping描述了某個橫向的地形剖面，其x、y、z坐標可映射到二維平面。基于道格拉斯-普克的總體最小二乘三維點抽稀算法的基本步驟如下：

Step1：將多波束測深ping點三維坐標（x，y，z）變換為二維曲線坐標（dx，dy）。計算ping條帶上某點O到空間直線AB（ping條帶上的首尾點）的垂足坐標N（xn，yn，zn），然后求得N到點A和C的距離dx、dy，其中A（x1，y1，z1）、B（x2，y2，z2）為：

Step2：利用RDP算法對曲線坐標作壓縮抽稀，選取抽稀后的特征點（Mx，My），并標記特征點所對應的索引值index。

Step3：搜索index和index+1中所包含的二維曲線點坐標對。考慮到TLS至少需要3個點坐標，當點坐標對數大于2時，利用TLS求出系數a、b。遍歷所有index，求出直線擬合的交點坐標集（Hi，Hj）（i，j=1，2，…），記錄其序號。

Step4：計算直線擬合的交點坐標（Hi，Hj）和特征點（Mx，My）間距離d。如果d大于限差ε，將直線擬合交點坐標替換成特征點，再運用基于Shapely庫函數中的曲線intersection方法求出擬合直線段和原始曲線的交點。若某段出現多個交點，排序得出距特征點最近的交點，并補充到直線擬合交點坐標集中；否則只保留直線擬合交點坐標。

Step5：用步驟4得到的二維點集，解算空間平面坐標方程，反算得到抽稀后的三維ping點數據。

算法用Python腳本實現后，添加到ArcToolBox中，如圖2。

圖2 總體最小二乘算法點抽稀工具界面

2.3 實 驗

選取某江水下地形多波束測深數據ping條帶數據做測試，驗證算法的有效性和正確性。該ping帶含有62個測深點。用ArcToolbox中的自定義點抽稀工具腳本進行實例計算。

根據國際海道測量規范s44 特等測量的精度要求,抽稀閾值設為0．26 m[8]。圖3顯示了垂距限差取0．26 m時，分別采用道格拉斯-普克算法和分段總體最小二乘的壓縮效果。其中藍色實線代表原始ping數據，黑色點表示RDP算法抽稀效果，紅色點表示TLS點抽稀效果。

圖3 點抽稀效果圖

為評定算法的抽稀精度，計算RDP算法和TLS算法的均方誤差，公式為其中di表示原始曲線上第i個點到壓縮后曲線的垂直距離，n表示總點數。垂距限差分別取0．26、0．5、0．1進行測試。TLS和RDP的壓縮比近似相等，均方根誤差為當垂距限差較大時,TLS點抽稀算法比RDP算法精度提高更明顯，能更好地逼近原始數據。

信息科學版,2013,38(7):850-856

[3] 盧銀宏,岳東杰，宋飛鳳．基于總體最小二乘的Douglas-Peucker算法在多波束測深數據抽稀中的應用[J]．水利與建筑工程學報,2012,10(2):4-5

[4] 彭海波,向洪普．基于Python的空間數據批量處理方法[J]．測繪與空間地理息,2011,34(4):81-82

[5] Dobesova Z．Programming Language Python for Data Processing[C]．International Conference on Electrical and Control Engineering (ICECE), 2011

[6] 張若愚．Python科學計算[M]． 北京：清華大學出版社，2012

[7] 魯鐵定,周世．總體最小二乘的迭代解法[J]．武漢大學學報：信息科學版,2010,35(11):1 351-1 354

[8] 夏偉,黃謨濤,劉雁春,等．Douglas-Peucker算法在多波束測深數據抽稀中的應用[J]．測繪科學,2009,34(3):159-160

[1] 楊云,孫群,朱長青．曲線數據壓縮的總體最小二乘算法[J]．西安電子科技大學學報：自然科學版,2008,35(5):946-950

[2] 王樂洋,許才軍． 總體最小二乘研究進展[J]．武漢大學學報：

P229.3

B

1672-4623（2014）04-0117-03

10.11709/j.issn.1672-4623.2014.04.041

張明希，碩士，主要研究方向為地理信息系統開發與應用。

2013-10-10。

項目來源：國家自然科學基金資助項目（41101374）。