楊俊峰，張丕狀

0 引言

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

1 目標定位原理

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

1.1 DTOA定位原理

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

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

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

式中，t0i表示信源信號到達各分站與主站之間的時間差，c表示信號傳播的速度。

1.2 DOA定位原理

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

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

假設發(fā)射信號的波長為λ，與天線視線軸的夾角為β，則到達相位干涉儀兩個測向天線1、2間的相位差為

式中，L為兩個測向天線之間的基線長度，λ為入射信號的波長，β為入射信號與天線視線軸的夾角。由式 (2)可以得到信號的波達角為

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

2 算法

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

利用相位干涉儀測向方法，各輔站測得的目標的俯仰角分別為

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

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

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

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

2.2 牛頓迭代定位算法

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

其雅克比矩陣為

當雅克比矩陣為非奇異陣，則震源坐標用牛頓迭代法表示為

2 計算機仿真及分析

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

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

表1 震動波傳播速度測算Tab.1 Measurement of vibratioon wave propagation velocity

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

圖3 局部震動信號Fig.3 Local vibration signal

表2 震源定位結果Tab.2 The source location results

利用時窗能量比法判斷初至波到時，對傳感器1接收到的震動信號進行預處理及時窗能量比，如圖4所示。

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

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

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

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

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

3 震源定位精度分析

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

圖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

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

以4個傳感器基站為例，在相同的傳感器基站站址測量誤差a=0.01 m和初至波達到時差測量誤差b=50 us的情況下，對不同布設的定位精度進行仿真，如圖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 結語

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

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

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

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

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

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

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

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

鄒翔宇，徐翊峰.2010.基于無線傳感器網(wǎng)絡的礦震震源定位［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.