楊海東，白書(shū)華

(南昌理工學(xué)院，江西南昌 330044)

0 引 言

多普勒測(cè)速儀(Doppler Velocity Log，DVL)是一種測(cè)量船舶或水下平臺(tái)相對(duì)于水底航速的儀器，其利用向水底發(fā)射聲波和接收水底反射回波之間的頻率變化，計(jì)算平臺(tái)在橫向和縱向的速度[1-3]。多普勒測(cè)速儀進(jìn)行信號(hào)處理的關(guān)鍵是多普勒頻率的測(cè)量，其直接影響了測(cè)速精度。早期測(cè)頻方法主要是過(guò)零點(diǎn)檢測(cè)方法，包括雙計(jì)數(shù)器過(guò)零點(diǎn)檢測(cè)、數(shù)字式過(guò)零點(diǎn)檢測(cè)方法等[4]。當(dāng)前測(cè)頻主要采用復(fù)相關(guān)法，其主要依據(jù)是回波信號(hào)復(fù)相關(guān)的相位是多普勒頻率的一次函數(shù)[5]。根據(jù)采用的信號(hào)形式，多普勒測(cè)速儀可分為三類(lèi)：一是脈沖非相干技術(shù)，即采用CW脈沖信號(hào)；二是脈沖相干技術(shù)，采用CW脈沖對(duì)；三是寬帶編碼技術(shù)，采用 M 序列或其它偽隨機(jī)序列調(diào)制信號(hào)[6]。

不論采用何種信號(hào)形式，對(duì)回波的復(fù)相關(guān)處理后，計(jì)算其相位時(shí)會(huì)得到多普勒頻移的一次函數(shù)，對(duì)其斜率的估計(jì)誤差將直接影響多普勒測(cè)頻精度，從而影響測(cè)速精度。影響相位斜率估計(jì)精度的因素有三點(diǎn)：一是噪聲的引入；二是波束開(kāi)角造成區(qū)域多目標(biāo)回波的疊加；三是復(fù)相關(guān)處理時(shí)對(duì)回波的截取方式，這涉及到回波的底跟蹤處理。復(fù)相關(guān)相位斜率估計(jì)的傳統(tǒng)方法是使用最小二乘法，其本質(zhì)是使誤差的均方根達(dá)到最小。最小二乘法適用于誤差較小時(shí)，當(dāng)誤差較大時(shí)其估計(jì)結(jié)果與真值差異較大。隨機(jī)樣本一致(Random Sample Consensus，RANSAC)算法從一個(gè)樣本集中建立數(shù)學(xué)模型來(lái)估計(jì)有關(guān)參數(shù)，其優(yōu)點(diǎn)在于該方法可以適用于樣本集中含有大量異常數(shù)據(jù)的情況[7-8]。

1 多普勒頻移和回波表達(dá)式

1.1 多普勒頻移的兩種基本形式

聲波由聲源發(fā)出，由接收器接收，由此產(chǎn)生的多普勒頻移有兩種基本形式。一種是聲源相對(duì)接收器運(yùn)動(dòng)而接收器靜止，另一種是聲源靜止，接收器相對(duì)聲源運(yùn)動(dòng)。這兩種產(chǎn)生頻移的表達(dá)式分別為

式中：f0為發(fā)射聲波頻率；c為聲速；v1和f1分別為聲源靜止、接收器運(yùn)動(dòng)時(shí)的相對(duì)速度和接收回波的頻率；v2和f2分別為接收器靜止、聲源運(yùn)動(dòng)時(shí)的相對(duì)速度和接收回波頻率。相向運(yùn)動(dòng)時(shí)頻率增大，相背運(yùn)動(dòng)時(shí)頻率減小。

1.2 單目標(biāo)回波的一般表達(dá)式

實(shí)際應(yīng)用中聲源和接收器合置的情況更為普遍，此時(shí)聲源發(fā)射聲波，遇到目標(biāo)后反射，接收器接收回波并提取多普勒頻移。示意圖如圖1所示。

圖1 聲源與接收器合置的多普勒示意圖Fig.1 Doppler diagram for transceiver platform

聲源和目標(biāo)的速度分別為v1和v2，脈沖寬度為T(mén)，頻率為f0，初始相位為?0，脈沖幅度為A0。聲基陣發(fā)射信號(hào)為

由于多普勒頻移，可得到接收回波的信號(hào)形式為

其中：A為接收回波幅度；為接收時(shí)刻。與發(fā)射時(shí)刻t的關(guān)系為

式中，R為發(fā)射起始時(shí)刻聲源與目標(biāo)的距離。

對(duì)照發(fā)射聲波信號(hào)表達(dá)式(2)，可見(jiàn)回波頻率發(fā)生了改變，頻移量fd為

式中，fe為回波頻率。

1.3 多普勒測(cè)速儀回波表達(dá)式

針對(duì)多普勒測(cè)速儀，設(shè)其水平向前的速度為vx，指向前下方波束的俯角為(波束朝前下方發(fā)射，波束中心線(xiàn)與水平 向夾角為俯角)，則v1=vxcos 。而海底為靜止?fàn)顟B(tài)，前述的回波信號(hào)表達(dá)式中v2為0。

將v1=vxcos和v2=0代入式(3)、(4)、(5)可得

2 復(fù)相關(guān)測(cè)頻方法

2.1 基本流程

多普勒測(cè)頻有多種方法，當(dāng)前主流的方法是使用復(fù)相關(guān)法。該方法的主要依據(jù)是復(fù)相關(guān)相位的斜率為多普勒頻移。方法的主要步驟是：

① 正交解調(diào)

② 低通濾波

③ 做復(fù)相關(guān)處理

④ 求復(fù)相關(guān)相位的斜率

處理流程如圖3所示。

圖2 復(fù)相關(guān)法測(cè)頻流程圖Fig.2 Frequency measuring flow by complex correlation method

2.2 C W單脈沖的復(fù)相關(guān)形式

當(dāng)發(fā)射信號(hào)為CW單脈沖時(shí)，發(fā)射信號(hào)為

其中：T為發(fā)射信號(hào)脈沖寬度；f0為信號(hào)頻率；?0為初始相位。

接收回波經(jīng)過(guò)時(shí)延和多普勒頻移fd后，其信號(hào)形式可記為

正交解調(diào)(IQ)是指同相支路(I路)和正交支路(Q路)分別與回波信號(hào)相乘。同相支路是回波信號(hào)與相乘，為

兩路信號(hào)合并為復(fù)數(shù)信號(hào)，即

其復(fù)相關(guān)為

對(duì)cτ的取值范圍進(jìn)行分析，可得

由此可見(jiàn)，復(fù)相關(guān)的相位是多普勒頻移的一次函數(shù)。

2.3 CW脈沖對(duì)的復(fù)相關(guān)形式

設(shè)發(fā)射聲波是兩個(gè)先后的CW脈沖組成的脈沖對(duì)，發(fā)射信號(hào)為

其中：T為發(fā)射信號(hào)脈沖寬度；f0為信號(hào)頻率；?0為初始相位；Tτ為兩個(gè)脈沖的起始時(shí)刻之間的時(shí)間間隔。

接收回波經(jīng)過(guò)時(shí)延τ和多普勒頻移fd，并假設(shè)在短時(shí)內(nèi)回波中兩個(gè)脈沖的幅度仍然保持一致，且相位也沒(méi)有差異。其信號(hào)形式可記為

經(jīng)過(guò)正交解調(diào)、低通濾波后形成的復(fù)數(shù)信號(hào)形式為

其復(fù)相關(guān)為

當(dāng)Tτ> 2T時(shí)，復(fù)相關(guān)函數(shù)分為3段考慮：

① 當(dāng) ?T<τc

由上述分析可知，CW脈沖對(duì)形式回波的復(fù)相關(guān)函數(shù)包括三個(gè)分段項(xiàng)，各個(gè)分段項(xiàng)均符合“復(fù)相關(guān)函數(shù)的相位的斜率為多普勒頻移”這個(gè)結(jié)論。實(shí)際上中間一段為互相關(guān)，兩邊兩端為自相關(guān)。

采用CW脈沖對(duì)進(jìn)行測(cè)頻，主要是利用其中的互相關(guān)段，而采用CW單脈沖的復(fù)相關(guān)處理是自相關(guān)。與自相關(guān)相比，互相關(guān)的優(yōu)點(diǎn)主要有：

① 減小目標(biāo)區(qū)域中強(qiáng)點(diǎn)對(duì)測(cè)頻精度的影響；

② 減小某些彩色噪聲的影響。

3 利用RANSAC方法提高測(cè)頻精度

復(fù)相關(guān)法測(cè)頻的核心在于求復(fù)相關(guān)相位的斜率來(lái)估計(jì)多普勒頻率。受多目標(biāo)回波疊加以及噪聲影響，實(shí)際回波的復(fù)相關(guān)相位并不是理想的線(xiàn)段，而包含了多種誤差。傳統(tǒng)上使用最小二乘法在誤差樣本中估計(jì)線(xiàn)段，這適用于誤差較小的情況。當(dāng)誤差較大時(shí)，尤其當(dāng)樣本中存在錯(cuò)誤值或異常點(diǎn)時(shí)，使用最小二乘法的結(jié)果會(huì)與真值產(chǎn)生較大差別。

RANSAC方法[7-8]利用樣本的內(nèi)在約束關(guān)系剔除錯(cuò)誤樣本。其目的是從一個(gè)數(shù)據(jù)集S中估計(jì)出一個(gè)較為精確的模型，其基本思路如下：

(1) 隨機(jī)從數(shù)據(jù)集S中選取s個(gè)點(diǎn)組成一個(gè)樣本，使用該樣本計(jì)算樣本模型；

(2) 利用上述樣本模型確定距離閾值dth內(nèi)的數(shù)據(jù)點(diǎn)集Si，即為該樣本模型的一個(gè)支撐。Si是S的一個(gè)子集。Si中的數(shù)據(jù)點(diǎn)稱(chēng)為內(nèi)點(diǎn)，Si補(bǔ)集中的數(shù)據(jù)點(diǎn)稱(chēng)為外點(diǎn)；

(3) 如果Si中內(nèi)點(diǎn)數(shù)量大于數(shù)量閾值Nth，則使用Si的所有點(diǎn)重新估計(jì)模型并結(jié)束；

(4) 如果Si中內(nèi)點(diǎn)數(shù)量小于閾值Nth，則選擇一個(gè)新的樣本子集，重復(fù)以上過(guò)程。

經(jīng)過(guò)N次試驗(yàn)選擇最大的支撐集Si，使用Si的所有點(diǎn)重估模型。

在多普勒測(cè)頻中，這個(gè)數(shù)據(jù)集即是回波信號(hào)的復(fù)相關(guān)相位序列，或是該序列中的一段數(shù)據(jù)，樣本模型即是待求參數(shù)多普勒頻移。由于理想情況下多普勒頻移是復(fù)相關(guān)相位的一次函數(shù)，因此該方法即是從樣本集中估計(jì)直線(xiàn)，其基本步驟如下：

(1) 隨機(jī)選擇兩點(diǎn)，確定一條直線(xiàn)；

(2) 將到該直線(xiàn)距離小于閾值的點(diǎn)確定為內(nèi)點(diǎn)，其余點(diǎn)設(shè)為外點(diǎn)；

(3) 如果內(nèi)點(diǎn)數(shù)量小于閾值，重新選擇樣本；

(4) 內(nèi)點(diǎn)數(shù)量大于閾值后取樣結(jié)束，使用內(nèi)點(diǎn)集利用最小二乘法重新估計(jì)直線(xiàn)。

估計(jì)得到的直線(xiàn)斜率即為多普勒頻移。

4 仿真結(jié)果

為了驗(yàn)證上述方法的有效性，進(jìn)行區(qū)域多目標(biāo)回波仿真，并進(jìn)行測(cè)速處理?；囅蛐鼻跋路桨l(fā)射波束，這是典型的Janus X型配置中的一個(gè)波束。示意圖如圖3所示。采用CW脈沖對(duì)的信號(hào)形式進(jìn)行仿真，具體參數(shù)如表1所示。發(fā)射波形、回波波形如圖4、5所示。

圖3 發(fā)射波束幾何示意圖Fig.3 Geometry of transmitting beam

表1 仿真參數(shù)表Table 1 Simulation parameter table

圖4 發(fā)射信號(hào)波形Fig.4 Transmitted signal waveform

利用復(fù)相關(guān)法求得的相位如圖6所示。從圖6中可以看出，復(fù)相關(guān)中包含了三個(gè)近似直線(xiàn)的段落(圖6中A、B、C三段)。其中A段和C段分別是脈沖對(duì)中前一個(gè)脈沖和后一個(gè)脈沖的回波自相關(guān)，B段為兩個(gè)脈沖回波的互相關(guān)，實(shí)際中用B段進(jìn)行處理。

圖6 復(fù)相關(guān)的相位圖Fig.6 Phase diagram of complex correlation

得到復(fù)相關(guān)相位后，提取其中的互相關(guān)段落進(jìn)行斜率估計(jì)，即得多普勒頻移，利用測(cè)得多普勒頻移計(jì)算得到基陣速度。為了驗(yàn)證本文方法的有效性，進(jìn)行200次仿真，統(tǒng)計(jì)最小二乘法和本文方法估計(jì)的均值及均值的標(biāo)準(zhǔn)差，并分為兩種情形。

情形一是嚴(yán)格提取互相關(guān)段落，仿真結(jié)果如圖7所示。

情形二是模糊提取互相關(guān)段落，使得出現(xiàn)部分野值。仿真結(jié)果如圖8所示。

對(duì)上述結(jié)果進(jìn)行定量統(tǒng)計(jì)，結(jié)果如表2所示。

圖7 嚴(yán)格提取互相關(guān)段落的估計(jì)速度Fig.7 Estimated velocity using the strictly extracted cross correlation section

圖8 模糊提取互相關(guān)段落的估計(jì)速度Fig.8 Estimated velocity using the fuzzy extracted cross correlation section

表2 最小二乘法和RANSAC估計(jì)速度均值和標(biāo)準(zhǔn)差Table 2 Mean and standard deviation of the velocity estimated by Least Square Method and RANSAC

從表2中可以看出，當(dāng)嚴(yán)格提取互相關(guān)段落時(shí)，最小二乘法和 RANSAC表現(xiàn)相當(dāng)；當(dāng)模糊提取互相關(guān)段落時(shí)，從標(biāo)準(zhǔn)差可以看出，RANSAC的穩(wěn)定性遠(yuǎn)勝于最小二乘法。

5 結(jié) 論

多普勒測(cè)速儀信號(hào)處理的關(guān)鍵是測(cè)頻，當(dāng)前主流的測(cè)頻方法是復(fù)相關(guān)法，該方法利用復(fù)相關(guān)相位是多普勒頻移的一次函數(shù)來(lái)估計(jì)多普勒頻移，進(jìn)而計(jì)算得到平臺(tái)速度。傳統(tǒng)的最小二乘法估計(jì)頻移適用于誤差較小且較為集中的情形，當(dāng)出現(xiàn)野值或誤差較大，尤其是復(fù)相關(guān)相位中的互相關(guān)段落提取較為模糊時(shí)，最小二乘法估計(jì)結(jié)果的穩(wěn)定性較差。為此提出了利用RANSAC來(lái)估計(jì)多普勒頻移的方法，該方法在復(fù)相關(guān)相位中存在較大誤差、野值時(shí)也有很好的穩(wěn)定性。通過(guò)仿真驗(yàn)證了本文提出方法的有效性。

[1]Krishnamurthy P, Khorrami F. A self-aligning underwater navigation system based on fusion of multiple sensors including DVL and IMU[C]//2013 9th Asian Control Conference (ASCC), 2013.

[2]Jerker Taudien, Sven G Bilén. Quantifying error sources that affect long-term accuracy of Doppler velocity logs[C]//OCEANS 2016 MTS/IEEE Monterey, 2016.

[3]HU Zhen, WANG Lei, YE Cong, et al. Research of acoustic survey devices applied on deepsea vehicles[C]//2016 IEEE/OES China Ocean Acoustics (COA), 2016: 1-5.

[4]田坦. 水下定位與導(dǎo)航技術(shù)[M]. 北京: 國(guó)防工業(yè)出版社, 2007.

TIAN Tan. Underwater positioning and navigation technologies[M]. Beijing: National Defense Industry Press, 2007.

[5]杜利利. 多普勒計(jì)程儀測(cè)速技術(shù)研究與實(shí)現(xiàn)[D]. 北京: 中國(guó)科學(xué)院大學(xué), 2011.

DU Lili. Research and realization of velocity-measuring for doppler velocity log[D]. Beijing: University of Chinese Academy of Sciences, 2011.

[6]馬海濤. 多普勒計(jì)程儀測(cè)流技術(shù)研究[D]. 北京: 中國(guó)科學(xué)院大學(xué),2011.

MA Haitao. Research of velocity estimation by DVL[D]. Beijing:University of Chinese Academy of Sciences, 2011.

[7]Peter C Niedfeldt, Kyle Ingersoll, Randal W Beard. Comparison and analysis of recursive-ransac for multiple target tracking[J].IEEE Transactions on Aerospace and Electronic Systems,2017(99): 1-1.

[8]Zahra Hossein-nejad, Mehdi Nasri. Image registration based on SIFT features and adaptive RANSAC transform[C]//2016 International Conference on Communication and Signal Processing (ICCSP), 2016: 1087-1091.