本發(fā)明屬于測量淤泥厚度技術(shù)領(lǐng)域，具體涉及一種基于克里金插值的水岸帶淤泥厚度計(jì)算及出圖方法。

背景技術(shù)：

湖濱水岸帶是介于陸生生態(tài)系統(tǒng)與水生生態(tài)系統(tǒng)之間的生態(tài)交錯類型之一，生產(chǎn)力高，生態(tài)邊緣效應(yīng)顯著。湖濱水岸帶河底土壤長期處于厭氧與好氧交替的環(huán)境中，強(qiáng)烈積累有機(jī)物質(zhì)并且分解緩慢，故有機(jī)質(zhì)豐富。對水岸帶河底淤泥的厚度測定對研究水岸帶土壤特性及活性有機(jī)碳采樣提取有著重要的參考意義。原始的淤泥厚度測量方法，一般情況下操作人員必須下水，直接面對河水及雜草、河底淤泥，工作環(huán)境比較惡劣，承載的測量船只也有傾覆的危險(xiǎn)。

近年來，測繪無人船技術(shù)得到了較大發(fā)展，無人船不僅能夠?qū)Υ钶d水質(zhì)監(jiān)測儀，對水質(zhì)進(jìn)行實(shí)時(shí)監(jiān)測，也能夠?qū)拥走M(jìn)行深度探測。而利用測繪無人船對于淤泥厚度的研究，由于對設(shè)備要求較高，目前工作中還缺乏實(shí)際的操作經(jīng)驗(yàn)。

技術(shù)實(shí)現(xiàn)要素：

發(fā)明目的：針對現(xiàn)有技術(shù)中存在的不足，本發(fā)明的目的是提供一種基于克里金插值的水岸帶淤泥厚度計(jì)算及出圖方法，提高測繪的工作效率和工作準(zhǔn)確度。

技術(shù)方案：為了實(shí)現(xiàn)上述發(fā)明目的，本發(fā)明采用的技術(shù)方案為：

一種基于克里金插值的水岸帶淤泥厚度計(jì)算及出圖方法：首先使用船載RTK和聲吶進(jìn)行點(diǎn)數(shù)據(jù)采集；然后，分別對高頻和低頻點(diǎn)數(shù)據(jù)進(jìn)行普通克里金格網(wǎng)插值，形成兩個頻段的規(guī)則格網(wǎng)數(shù)據(jù)；最后，結(jié)合桿高等固定參數(shù)設(shè)置，計(jì)算淤泥頂和淤泥底數(shù)據(jù)，得到淤泥厚度，并使用ArcGIS進(jìn)行柵格計(jì)算，得到空間淤泥厚度圖。

所述的基于克里金插值的水岸帶淤泥厚度計(jì)算及出圖方法，包括以下步驟：

1)通過岸端服務(wù)器進(jìn)行測繪無人船進(jìn)行路線及目標(biāo)點(diǎn)設(shè)置；確定河道的起始點(diǎn)和終止點(diǎn)，然后設(shè)置無人船測量的時(shí)間間隔為1次/秒，采樣路線采用Z字形采樣路線；將預(yù)先設(shè)定好的路線輸入控制系統(tǒng)，讓無人船能夠進(jìn)行自動巡航和測量數(shù)據(jù)，所測量得到的數(shù)據(jù)是以CAD軟件的DXF格式保存；

2)無人船將所測量的CAD的DXF點(diǎn)狀數(shù)據(jù)，通過無線傳輸?shù)姆绞?，傳輸?shù)桨抖朔?wù)器控制電腦中，利用ArcGIS軟件的格式轉(zhuǎn)換工具，轉(zhuǎn)為shp格式；

3)設(shè)置克里金參數(shù)并插值，生成格網(wǎng)數(shù)據(jù)；

4)利用克里金插值計(jì)算得到的高頻和低頻信號的柵格數(shù)據(jù)，結(jié)合桿高、水面高、無人船吃水深度固定參數(shù)設(shè)置，計(jì)算淤泥頂部和底部數(shù)據(jù)；

5)在獲得淤泥底高和頂高的基礎(chǔ)上，將頂高減去底高，就得到了淤泥的厚度圖；

6)最后，將實(shí)際人工測量數(shù)據(jù)與所得到的淤泥厚度進(jìn)行驗(yàn)證，并使用ArcGIS軟件進(jìn)行出圖。

步驟3)中，使用克里金插值將測量點(diǎn)構(gòu)建成測量面；將所要測試的內(nèi)河河道段區(qū)域設(shè)定為B，將需要構(gòu)面的河道淤泥層測量值設(shè)定為Z(x)，因此有{Z(x)∈B}，Z(x)是一個二階平穩(wěn)的隨機(jī)函數(shù)，它在空間中的取值設(shè)定為Z(x₁),Z(x₂),...,Z(x_n)，其中x表示淤泥測量點(diǎn)的空間位置；按照普通克里金算法的原理，非測量點(diǎn)x₀的高頻和低頻值的Z(x₀)估計(jì)值是多個已知測量點(diǎn)的加權(quán)求和，如下式所示：

其中，Z(x_i)(i＝0,1,2,...,n)為測量點(diǎn)x_i的值，x₀為非測量點(diǎn)，也就是需要插值得到的淤泥層值，其他為已知點(diǎn)；λ_i(i＝1,2,...,n)為權(quán)值；權(quán)值并非僅由距離決定，是需要在最小方差及無偏特征的假設(shè)條件下，由變異函數(shù)計(jì)算數(shù)值決定；

普通克里金方程組：

按照普通克里金估計(jì)方差最小，可以變形成為如下所示：

其中，C表示協(xié)方差函數(shù)，E{}為數(shù)學(xué)期望值。

步驟4)中，使用的克里金插值算法進(jìn)行淤泥測量點(diǎn)計(jì)算，步驟具體如下：

1)對于觀測數(shù)據(jù)，兩兩計(jì)算距離與半方差；

2)尋找一個擬合曲線擬合距離與半方差的關(guān)系，從而能根據(jù)任意距離計(jì)算出相應(yīng)的半方差；

3)計(jì)算出所有已知點(diǎn)之間的半方差r_ij；

4)對于未知點(diǎn)z_o，計(jì)算它到所有已知點(diǎn)z_i的半方差r_io；

5)求解步驟4)的方程組，得到最優(yōu)系數(shù)λi；

6)使用最優(yōu)系數(shù)對已知點(diǎn)的屬性值進(jìn)行加權(quán)求和，得到未知點(diǎn)z_o的估計(jì)值。

步驟5)中，具體公式如下所示：

H_淤泥頂高＝(H_水面高程-H_高頻-(H_桿高+H_吃水高度))

H_淤泥底高＝(H_水面高程-H_低頻-(H_桿高+H_吃水高度))

H_淤泥厚度＝H_淤泥頂高-H_淤泥底高。

有益效果：與現(xiàn)有技術(shù)相比，本發(fā)明的基于克里金插值的水岸帶淤泥厚度計(jì)算及出圖方法，將無人船應(yīng)用于河道的淤泥厚度測繪中，能夠提高測繪的工作效率和工作準(zhǔn)確度。通過對南京水西門河段進(jìn)行實(shí)地測試和出圖驗(yàn)證，本方法具有較高的精度和效率，能夠較好的完成典型水岸帶河底淤泥的精確測定，能夠?yàn)樗稁寥姥芯刻峁┘夹g(shù)支持。

附圖說明

圖1是基于克里金插值的水岸帶淤泥厚度計(jì)算及出圖方法的流程圖；

圖2是克里金插值結(jié)果圖；其中，a為高頻面結(jié)果圖，b為低頻面結(jié)果圖；亮度越暗，表示深度越深；

圖3是河道淤泥厚度圖。

具體實(shí)施方式

下面結(jié)合具體實(shí)施例對本發(fā)明做進(jìn)一步的說明。

以下實(shí)施例中所使用的測繪無人船，長約1.5米，寬0.7米，船身安裝了RTK裝置及雙頻聲吶裝置。

RTK(Real Time Kinematic，載波相位差分)定位技術(shù)是基于載波相位觀測值的實(shí)時(shí)動態(tài)定位技術(shù)，能夠?qū)崟r(shí)地提供測站點(diǎn)在指定坐標(biāo)系中的三維定位結(jié)果，并達(dá)到厘米級精度。精度高，作業(yè)方便。RTK作業(yè)不受通視條件限制，無需做控制，基準(zhǔn)站設(shè)置好，進(jìn)行點(diǎn)檢核后，即可開測，如用虛擬基站則更簡便。由于精度達(dá)到厘米級，所以可以用于河道淤泥的測量點(diǎn)的精確定位。

雙頻聲吶，指的是能夠在同一時(shí)段同一地點(diǎn)發(fā)射2個頻段的聲吶信號，利用聲波的頻率和強(qiáng)度存在的差異，測定所要求的被測物體的形態(tài)和距離等參數(shù)特征。由于聲波信號的頻率不同，其在水下的穿透力和強(qiáng)度也存在差異。高頻聲吶信號設(shè)定為180MHz，穿透力弱，信號強(qiáng)，能夠穿透內(nèi)陸河道中的浮游植物和懸浮垃圾等，但是無法穿透淤泥頂部，可以用來測定聲吶距離淤泥頂部的長度；低頻聲吶信號設(shè)定為50MHz，穿透力較強(qiáng)，信號偏弱，它不僅能夠輕易穿透浮游植物和懸浮垃圾，也能夠穿透淤泥層，直接射到淤泥底部的河床硬質(zhì)層。由于河床硬質(zhì)層的密度大，低頻聲吶信號無法穿透，因此可以用來測定聲吶距離河床硬質(zhì)層的長度，也就是距離淤泥底部的長度。

實(shí)施例1

一種基于克里金插值的水岸帶淤泥厚度計(jì)算及出圖方法，流程圖如圖1的所示，具體步驟如下：

1)測量淤泥樣本點(diǎn)數(shù)據(jù)，生成CAD數(shù)據(jù)。

通過岸端服務(wù)器進(jìn)行測繪無人船進(jìn)行路線及目標(biāo)點(diǎn)設(shè)置。確定河道的起始點(diǎn)和終止點(diǎn)，然后設(shè)置無人船測量的時(shí)間間隔為1次/秒，采樣路線采用Z字形采樣路線。將預(yù)先設(shè)定好的路線輸入控制系統(tǒng)，讓無人船能夠進(jìn)行自動巡航和測量數(shù)據(jù)，所測量得到的數(shù)據(jù)是以CAD軟件的DXF格式保存。

2)將CAD數(shù)據(jù)轉(zhuǎn)換成shp點(diǎn)狀格式數(shù)據(jù)。

shp格式是shapefile格式的簡寫，該格式是由ESRI公司開發(fā)，一個shp文件包括一個主文件，一個索引文件，和一個dBASE表。之所以要將原始格式DXF轉(zhuǎn)換成shp格式，是因?yàn)锳rcGIS軟件具有強(qiáng)大的出圖功能，最后出圖將會使用ArcGIS軟件進(jìn)行操作，而ArcGIS軟件無法直接操作DXF文件，但是能夠直接操作shp格式數(shù)據(jù)。無人船將所測量的CAD的DXF點(diǎn)狀數(shù)據(jù)，通過無線傳輸?shù)姆绞?，傳輸?shù)桨抖朔?wù)器控制電腦中，利用ArcGIS軟件的格式轉(zhuǎn)換工具，轉(zhuǎn)為shp格式。

3)設(shè)置克里金參數(shù)并插值，生成格網(wǎng)數(shù)據(jù)。

由于克里金插值能夠?qū)Ⅻc(diǎn)狀數(shù)據(jù)生成柵格格網(wǎng)數(shù)據(jù)，并且能夠?qū)γ總€數(shù)據(jù)區(qū)域的精度進(jìn)行求解，所以采用克里金插值將點(diǎn)生成面。所使用的軟件工具是ArcGIS軟件，計(jì)算得到生成高頻和低頻規(guī)則格網(wǎng)的柵格數(shù)據(jù)，用于后期的求取層差，從而獲得淤泥厚度。

4)基于格網(wǎng)數(shù)據(jù)，計(jì)算淤泥頂高和淤泥底高的格網(wǎng)數(shù)據(jù)。

利用克里金插值計(jì)算得到的高頻和低頻信號的柵格數(shù)據(jù)，結(jié)合桿高、水面高、無人船吃水深度等固定參數(shù)設(shè)置，計(jì)算淤泥頂部和底部數(shù)據(jù)。

5)基于淤泥底高和頂高的格網(wǎng)數(shù)據(jù)，進(jìn)行差值計(jì)算，獲得淤泥厚度圖。

在獲得淤泥底高和頂高的基礎(chǔ)上，將頂高減去底高，就得到了淤泥的厚度圖。采用的具體操作是實(shí)用柵格相減。

6)進(jìn)行結(jié)果驗(yàn)證，利用ArcGIS軟件進(jìn)行出圖。

最后，將實(shí)際人工測量數(shù)據(jù)與所得到的淤泥厚度進(jìn)行驗(yàn)證，并使用ArcGIS軟件進(jìn)行出圖。

本專利方法使用克里金(Kriging)插值將測量點(diǎn)構(gòu)建成測量面?？死锝鸩逯凳紫瓤紤]的是空間屬性在二維空間中的分布，然后用此范圍內(nèi)的采樣點(diǎn)來估計(jì)待插點(diǎn)的屬性值。將所要測試的內(nèi)河河道段區(qū)域設(shè)定為B，將需要構(gòu)面的河道淤泥層測量值設(shè)定為Z(x)，因此有{Z(x)∈B}，Z(x)是一個二階平穩(wěn)的隨機(jī)函數(shù)，它在空間中的取值設(shè)定為Z(x₁),Z(x₂),...,Z(x_n)，其中x表示淤泥測量點(diǎn)的空間位置。按照普通克里金算法的原理，非測量點(diǎn)x₀的高頻和低頻值的Z(x₀)估計(jì)值是多個已知測量點(diǎn)的加權(quán)求和，如下式所示：

其中，Z(x_i)(i＝0,1,2,...,n)為測量點(diǎn)x_i的值，x₀為非測量點(diǎn)，也就是我們需要插值得到的淤泥層值，其他為已知點(diǎn)。λ_i(i＝1,2,...,n)為權(quán)值。權(quán)值并非僅由距離決定，是需要在最小方差及無偏特征的假設(shè)條件下，由變異函數(shù)計(jì)算數(shù)值決定。

詳細(xì)推導(dǎo)請參見文獻(xiàn)[靳國棟,劉衍聰,牛文杰.距離加權(quán)反比插值法和克里金插值法的比較[J].長春工業(yè)大學(xué)學(xué)報(bào),2003,24(3):53-57]。以下列出普通克里金方程組：

按照普通克里金估計(jì)方差最小，可以變形成為如下所示：

其中，C表示協(xié)方差函數(shù)，E{}為數(shù)學(xué)期望值。

本實(shí)施例使用的普通克里金插值算法進(jìn)行淤泥測量點(diǎn)計(jì)算的基本步驟，具體如下：

1)對于觀測數(shù)據(jù)，兩兩計(jì)算距離與半方差；

2)尋找一個擬合曲線擬合距離與半方差的關(guān)系，從而能根據(jù)任意距離計(jì)算出相應(yīng)的半方差；

3)計(jì)算出所有已知點(diǎn)之間的半方差r_ij；

4)對于未知點(diǎn)z_o，計(jì)算它到所有已知點(diǎn)z_i的半方差r_io；

5)求解步驟4)的方程組，得到最優(yōu)系數(shù)λi；

6)使用最優(yōu)系數(shù)對已知點(diǎn)的屬性值進(jìn)行加權(quán)求和，得到未知點(diǎn)z_o的估計(jì)值。

本專利將游離于河床表面，一定密度內(nèi)的底質(zhì)層稱之為河道淤泥層。通常其組成為粘性細(xì)顆粒泥沙。本實(shí)施例使用的雙頻聲吶采用的工作原理是：聲吶裝置同時(shí)使用2個工作頻率，180kHz和50kHz。其中，低頻由于其波長長，穿透力強(qiáng)，能夠到達(dá)河床底；高頻的穿透能力較弱，只能到達(dá)淤泥頂部。然后，利用兩者的高程差，就能夠計(jì)算出淤泥的厚度。具體公式如下所示：

H_淤泥頂高＝(H_水面高程-H_高頻-(H_桿高+H_吃水高度))

H_淤泥底高＝(H_水面高程-H_低頻-(H_桿高+H_吃水高度))

H_淤泥厚度＝H_淤泥頂高-H_淤泥底高

實(shí)施例2

以江蘇省南京市水西門河段為例，采用實(shí)施例1的方法進(jìn)行實(shí)地測量，本實(shí)施例使用的計(jì)算機(jī)硬件主要為：CPU是i5處理器，內(nèi)存是8GB，硬盤為固態(tài)硬盤。使用的軟件為：操作系統(tǒng)是windows 7專業(yè)版系統(tǒng)，ArcGIS軟件是9.3.1版本。

(1)首先是對采集到的點(diǎn)狀數(shù)據(jù)進(jìn)行格式轉(zhuǎn)化。為了使用ArcGIS軟件中的克里金插值算法，需要將原始格式轉(zhuǎn)化為shapefile格式數(shù)據(jù)。測繪無人船直接測量的數(shù)據(jù)格式是CAD的DXF格式，使用ArcGIS軟件將其打開后，選中“*.dxf Point”后，右鍵點(diǎn)擊，在彈出窗口中選擇Data，然后是Export Data，就可以將其轉(zhuǎn)化為shapefile格式。

(2)其次是構(gòu)建克里金規(guī)則格網(wǎng)柵格數(shù)據(jù)。由于高頻點(diǎn)和低頻點(diǎn)都是點(diǎn)格式，無法直接利用，也無法直接生成面數(shù)據(jù)，這就需要使用克里金插值算法構(gòu)建規(guī)則格網(wǎng)柵格面。使用ArcGIS軟件中的Geostatistical Analyst(地統(tǒng)計(jì)分析工具)進(jìn)行計(jì)算，采用的是經(jīng)典的普通克里金插值算法。插值后的高頻面和低頻面見圖2所示，其中，a為高頻面結(jié)果圖，b為低頻面結(jié)果圖；亮度越暗，表示深度越深。

(3)使用高頻面和低頻面的插值，計(jì)算淤泥厚度圖。如果需要計(jì)算淤泥頂高和淤泥底高，則需要知道參數(shù)桿高和吃水深度。由于只需要獲得淤泥厚度，因此無需知道桿高和吃水深度，直接用高頻面減去低頻面，就能得到淤泥厚度。使用ArcGIS軟件中的Raster Calculator進(jìn)行運(yùn)算，得到的結(jié)果如圖3所示。使用10級分級圖顯示，亮度越高，則淤泥厚度約薄，亮度越暗，則淤泥厚度越厚。從結(jié)果圖可知，河道的兩側(cè)的淤泥厚度更厚，而河道中央的淤泥厚度則相對較薄，這是由于南京水西門的外秦淮河的水流較為湍急，中間的淤泥都被沖刷到河道的兩側(cè)，這就為今后該河段的清淤工作指明了方向，清理淤泥需要從河道兩側(cè)重點(diǎn)處理。

表1所示的是低頻面、高頻面以及淤泥厚度面的統(tǒng)計(jì)值。從中可以發(fā)現(xiàn)，低頻面可以近似指代淤泥底部，均值為3.58米；高頻面可以近似指代淤泥頂部，均值為7.76米。淤泥頂部的標(biāo)準(zhǔn)差為0.01，說明其淤泥分布較為均勻。

表1實(shí)驗(yàn)區(qū)淤泥面統(tǒng)計(jì)值

可見，使用測繪無人船進(jìn)行河底淤泥厚度的測量，不僅能夠減少人力物力的投資，降低測量人員的測繪風(fēng)險(xiǎn)，也能夠提高測繪的效率和精度，是未來的發(fā)展趨勢。從測量的結(jié)果可知，本實(shí)驗(yàn)所測的水西門地區(qū)外秦淮河河道的兩側(cè)的淤泥厚度更厚，而河道中央的淤泥厚度則相對較薄。