基于DTOA/DOA和牛頓迭代法的震源定位方法研究*

楊俊峰，張丕狀

0 引言

隨著現(xiàn)代科學(xué)技術(shù)的發(fā)展，定位技術(shù)在航空、航天、交通、震源探測(cè)、海洋勘探等領(lǐng)域得到了廣泛地應(yīng)用。定位精度的要求越來越高，在硬件系統(tǒng)確定條件下，可以通過改進(jìn)算法提高系統(tǒng)的定位精度 (呂晶晶，姚金杰，2011;姚金杰，韓焱，2010)。本文提出了一種以牛頓迭代算法為基礎(chǔ)，利用無源時(shí)差定位 (DTOA)和測(cè)量波達(dá)方向(DOA)聯(lián)合定位得到微震源的估計(jì)位置的算法，本算法充分利用DTOA和DOA聯(lián)合定位算法的估計(jì)特性，能夠有效地解決迭代法的初始值問題，保證算法收斂并且提高迭代算法的收斂速度。通過設(shè)置建立的冗余函數(shù)表達(dá)式和結(jié)束迭代的精度范圍，作為判斷迭代是否結(jié)束的條件，這樣既能滿足一定的精度，又能提高計(jì)算速度，本文通過對(duì)計(jì)算機(jī)仿真、數(shù)據(jù)分析和計(jì)算時(shí)間來證明該算法的有效性。

1 目標(biāo)定位原理

為了充分利用各種觀測(cè)信息，以獲取更好的定位效果，我們將觀測(cè)數(shù)據(jù)中的角度信息與時(shí)間信息結(jié)合起來。本文運(yùn)用的DTOA和DOA聯(lián)合定位方法，主要利用主站增加對(duì)目標(biāo)俯仰角的測(cè)量，形成一個(gè)包含3個(gè)子系統(tǒng)的冗余定位系統(tǒng)，各子系統(tǒng)分別對(duì)目標(biāo)進(jìn)行定位。最后對(duì)定位結(jié)果進(jìn)行融合，將融合后的結(jié)果作為信源目標(biāo)的初始估計(jì)值，再通過牛頓迭代算法進(jìn)行精確定位。

1.1 DTOA定位原理

無源時(shí)差定位利用主站和輔站之間接收到輻射源信號(hào)的到達(dá)時(shí)間差來定位。獲得時(shí)差的信息之后，通過乘以波的傳播速度就可以得到信源到主站和輔站的距離差，在三維空間定位中，得到一組雙曲面方程組。在三維時(shí)差定位中至少需要4個(gè)測(cè)量站才能實(shí)現(xiàn)對(duì)目標(biāo)的定位 (陳玲，李少洪，2003)，如圖1所示。

圖1 三維時(shí)差定位示意圖Fig.1 Schematic diagram of 3D TDOA location

設(shè)震源位置坐標(biāo)為P(x，y，z)，基站坐標(biāo)分別為 Pi(xi，yi，zi)(i=0，1，2，3…n)。目標(biāo)到基站的距離為ri(i=0，1，2，3…n)。目標(biāo)到各分站與主站(基站P0為主站)的距離差為ri0。假設(shè)有4個(gè)基站，則根據(jù)時(shí)差定位原理有以下關(guān)系

式中，t0i表示信源信號(hào)到達(dá)各分站與主站之間的時(shí)間差，c表示信號(hào)傳播的速度。

1.2 DOA定位原理

測(cè)量波達(dá)方向技術(shù)在無源定位中占有極其重要的作用。隨著數(shù)字信號(hào)處理技術(shù)的發(fā)展，數(shù)字相位干涉儀的測(cè)向方法由于其具有速度快、成本低、精度高的特點(diǎn)在工程上獲得了廣泛的應(yīng)用。在原理上，相位干涉儀可以實(shí)現(xiàn)對(duì)單個(gè)脈沖測(cè)向(袁孝康，1999)。圖2顯示相位干涉儀 (單基線)測(cè)向原理圖。

圖2 單基線相位干涉儀側(cè)向原理圖Fig.2 Lateral principle graph of single-baseline phase interferometer

假設(shè)發(fā)射信號(hào)的波長(zhǎng)為λ，與天線視線軸的夾角為β，則到達(dá)相位干涉儀兩個(gè)測(cè)向天線1、2間的相位差為

式中，L為兩個(gè)測(cè)向天線之間的基線長(zhǎng)度，λ為入射信號(hào)的波長(zhǎng)，β為入射信號(hào)與天線視線軸的夾角。由式 (2)可以得到信號(hào)的波達(dá)角為

式中，φ是兩個(gè)陣元之間輸出的相位差，入射信號(hào)波長(zhǎng)的估計(jì)值λ由信號(hào)的頻率估計(jì)值計(jì)算得到，基線長(zhǎng)度L則可直接測(cè)量，在上述參數(shù)確定之后，就能根據(jù)式 (3)得到輻射源信號(hào)的入射角度β^。

2 算法

2.1 DTOA和DOA聯(lián)合定位算法求初始值

利用相位干涉儀測(cè)向方法，各輔站測(cè)得的目標(biāo)的俯仰角分別為

式中，αi(i=1，2，3)表示各輔站測(cè)得的俯仰角;z表示目標(biāo)輻射源的高度。

將式 (1)和式 (4)組成方程組 (其中，i=1，2)，稱為子系統(tǒng)Ⅰ，同理，將式 (1)和式(5)組成方程組 (其中，i=2，3)，稱為子系統(tǒng)Ⅱ，式 (1)和式 (6)組成方程組 (其中，i=1，3)，稱為子系統(tǒng)Ⅲ，這三個(gè)子系統(tǒng)均可完成對(duì)目標(biāo)的空間三維定位，我們選擇這三個(gè)子系統(tǒng)組成一個(gè)冗余定位系統(tǒng)，則該系統(tǒng)的定位方程組為:

其中，i=1，2，可改寫為

假設(shè)3個(gè)子系統(tǒng)所得到的目標(biāo)位置分別為X^i(i=1，2，3)，定位協(xié)方差矩陣分別為Pi(i=1，2，3)，各協(xié)方差矩陣之間互不相關(guān)。SWLS融合后的目標(biāo)位置為 (Petre，Li，2006;毛永毅，白菊蓉，2006)

2.2 牛頓迭代定位算法

牛頓法是一種使用導(dǎo)數(shù)的算法，每一步迭代方向都是沿著當(dāng)前點(diǎn)函數(shù)值下降的方向。因此，需設(shè)定合適的初值才能保證算法是否收斂及其收斂速度。根據(jù)牛頓迭代算法的條件，將所建立目標(biāo)定位模型轉(zhuǎn)換為 (柳輝，2007)

其雅克比矩陣為

當(dāng)雅克比矩陣為非奇異陣，則震源坐標(biāo)用牛頓迭代法表示為

2 計(jì)算機(jī)仿真及分析

在野外條件下，通過實(shí)驗(yàn)完成震源定位。首先測(cè)算震動(dòng)波在試驗(yàn)場(chǎng)地介質(zhì)中的傳播速度，再根據(jù)加速度傳感器接收到的震動(dòng)波數(shù)據(jù)，結(jié)合震源定位的原理及算法，得出震源的位置 (鄒翔宇，徐翊峰，2010)。

震動(dòng)波在該實(shí)驗(yàn)場(chǎng)地介質(zhì)中的傳播速度測(cè)算具體方法是:在距離待爆炸位置3 m、5 m、7 m、10 m、15 m的直線上分別依次布設(shè)5個(gè)加速度傳感器A、B、C、D、E，埋深均為0.5 m，在同一位置爆炸3次，通過數(shù)據(jù)采集儀采集震動(dòng)波數(shù)據(jù)，并將數(shù)據(jù)儲(chǔ)存在計(jì)算機(jī)里，利用MATLAB對(duì)信號(hào)進(jìn)行分析處理，計(jì)算出兩個(gè)加速度傳感器之間接收到震動(dòng)信號(hào)的時(shí)間差，根據(jù)兩個(gè)傳感器距離震源的距離差計(jì)算震動(dòng)波的傳播速度v=Δd/Δt。通過計(jì)算測(cè)得的震動(dòng)波傳播速度如表1所示。

表1 震動(dòng)波傳播速度測(cè)算Tab.1 Measurement of vibratioon wave propagation velocity

通過加速度傳感器接收到的震動(dòng)波分析處理可大致得到震動(dòng)波的傳播規(guī)律，在距震源約7 m處的震動(dòng)信號(hào)衰減為彈性波。以此為依據(jù)在待測(cè)震源周圍布設(shè)加速度傳感器，傳感器基站的位置分別為 (6，3， －0.5)，(－4，6， －0.4)，(－6，－5， －0.42)， (7， －4， －0.28)， (0，5， －0.38)。數(shù)據(jù)采集儀的采樣率為10 kHz，采樣時(shí)間為6.553 8 s，采用手動(dòng)觸發(fā)方式。在同一地點(diǎn)爆炸3次，通過測(cè)算的震動(dòng)波傳播速度，傳感器基站間接收到震動(dòng)波的時(shí)差及震動(dòng)定位算法以解算出震源位置。加速度傳感器接收到的局部放大信號(hào)如圖3所示。

圖3 局部震動(dòng)信號(hào)Fig.3 Local vibration signal

表2 震源定位結(jié)果Tab.2 The source location results

利用時(shí)窗能量比法判斷初至波到時(shí)，對(duì)傳感器1接收到的震動(dòng)信號(hào)進(jìn)行預(yù)處理及時(shí)窗能量比，如圖4所示。

以同樣的方法判斷其它各傳感器接收到的初至波到時(shí)，結(jié)合各傳感器基站位置及定位算法解算出的震源位置如表2所示。對(duì)運(yùn)用本文算法所得定位結(jié)果進(jìn)行誤差分析如表3所示。

在定位過程中，由于傳感器基站的布設(shè)及埋深較淺，使得Z方向的誤差較大，加之震動(dòng)波在傳播過程中引入的誤差，傳感器基站站址測(cè)量誤差，初至波到時(shí)估計(jì)誤差等因數(shù)影響導(dǎo)致震源位置與實(shí)際位置之間有偏差。

從表3可知，本文算法定位結(jié)果的均方根誤差值遠(yuǎn)遠(yuǎn)小于DTOA和DOA聯(lián)合定位、LS-Newton法及Chan-Newton法定位結(jié)果的均方根誤差值，定位結(jié)果比較理想。

圖4 判斷初至波時(shí)刻圖Fig.4 Determine the time of the first arrival

表3 定位結(jié)果誤差分析Tab.3 Error analysis of positioning results

3 震源定位精度分析

通過對(duì)震源定位進(jìn)行精度分析，可以了解影響定位精度的因素以及這些因素對(duì)定位精度產(chǎn)生的影響，因此對(duì)震源定位進(jìn)行精度分析是很必要的。下文對(duì)在不同的布站形式、基站數(shù)目、初至波到達(dá)時(shí)間差估計(jì)誤差、及在不同站址測(cè)量誤差情況下的定位精度變化情況進(jìn)行仿真分析，仿真結(jié)果可供工程參考。

圖5 不同誤差情況下星型布站GDOP分布圖(a)a=0.03 m，b=0.5 ms;(b)a=0.01 m，b=50 usFig.5 GDOP distribution in star-model station position under different error conditions

設(shè)震源范圍是x方向0～30 m，y方向是－5～+5 m，z方向－1 m，對(duì)4個(gè)傳感器基站為星型布局時(shí)進(jìn)行仿真分析。圖5a中存在a=0.03 m的傳感器基站站址測(cè)量誤差，b=0.5 ms的初至波到達(dá)時(shí)間差測(cè)量誤差。圖5b中存在a=0.01 m的傳感器基站站址測(cè)量誤差，b=50 us的初至波到達(dá)時(shí)間差測(cè)量誤差。由圖5b可以看出，減小傳感器基站站址和初至波達(dá)到時(shí)差測(cè)量誤差可有效地提高定位精度。

以4個(gè)傳感器基站為例，在相同的傳感器基站站址測(cè)量誤差a=0.01 m和初至波達(dá)到時(shí)差測(cè)量誤差b=50 us的情況下，對(duì)不同布設(shè)的定位精度進(jìn)行仿真，如圖6所示。從圖中4種布站下的GDOP圖中可知，星型布站的定位精度最高;平行四邊形布站的定位次之，其在X軸方向的定位精度要高于Y軸;菱形布站和倒三角布站的定位精度較差。

圖6 不同布站形式GDOP圖(a)星型布站;(b)平行四邊形布站;(c)菱形布站;(d)倒三角布站Fig.6 GDOP in different forms of station arrangement(a)star-model station position;(b)parallelogram-model station position;(c)dimod-model station position;(d)upside-down triangle-model station position

4 結(jié)語

本文結(jié)合DTOA和DOA聯(lián)合定位法估計(jì)性好和牛頓迭代法收斂速度快的優(yōu)點(diǎn)，提出了DTOA、DOA和牛頓迭代法相結(jié)合的地下震源定位算法。避免了DTOA和DOA聯(lián)合定位法受測(cè)量誤差影響較大和牛頓法對(duì)初始值選擇敏感的問題。經(jīng)過在野外條件下實(shí)驗(yàn)及仿真，實(shí)驗(yàn)結(jié)果表明本文算法有效地提高了震源定位的精度。在此基礎(chǔ)上對(duì)震源系統(tǒng)中存在不同站址測(cè)量誤差、不同初至波到達(dá)時(shí)差估計(jì)誤差、不同布站形式情況下的定位精度進(jìn)行了仿真分析，結(jié)果表明減小傳感器基站站址測(cè)量誤差和初至波達(dá)到時(shí)差測(cè)量誤差可有效地提高定位精度，星型布站的定位精度要優(yōu)于平行四邊形布站、菱形布站和倒三角布站的定位精度。

陳玲，李少洪.2003.無源側(cè)向測(cè)時(shí)差定位算法研究［J］．電子信息學(xué)報(bào)，(6):771－776.

柳輝.2007.解非線性方程組的牛頓迭代法及其應(yīng)用［J］．重慶工學(xué)院學(xué)報(bào)，21(8):95－98.

呂晶晶，姚金杰.2011.基于最小二乘法和牛頓迭代法的空中目標(biāo)定位［J］．微電子學(xué)于計(jì)算機(jī)，28(9):110－115.

呂明，郭士民.2007.基于數(shù)據(jù)融合的時(shí)差定位處理算法的應(yīng)用［J］．儀器儀學(xué)報(bào)，28(4):100－102.

毛永毅，白菊蓉.2006.空間四站時(shí)差定位中的模糊及無線研究［J］．電訊技術(shù)，(6):53－57.

姚金杰，韓焱.2010.基于粒子群和牛頓迭代法的目標(biāo)定位方法研究［J］．計(jì)算機(jī)應(yīng)用研究，27(5):1 700 －1 713.

袁孝康.1999.相位干涉儀側(cè)向定位研究［J］．上海航天，(3):1－7.

鄒翔宇，徐翊峰.2010.基于無線傳感器網(wǎng)絡(luò)的礦震震源定位［J］．煤礦安全，(10):53－57.

Carter G C.1987.Coherence and Time Delay Estimation［J］．Proc IEEE，75(2):236－255.

Petre S，Li J.2006.Source Localization from Range-Diffrernce Measurements［J］．IEEE Signal Processing Magzine，63(11):63 － 65.