奚 暢 蔡志明 袁 駿

(海軍工程大學(xué)電子工程學(xué)院 武漢 430033)

0 引言

拖線陣聲吶對目標(biāo)搜索、定位、跟蹤、攻擊過程中，拖船常需要進(jìn)行戰(zhàn)術(shù)機動，這必然造成拖線陣陣形畸變。如果在陣列信號處理時仍假設(shè)拖線陣直線形態(tài)，將嚴(yán)重影響拖線陣對弱目標(biāo)的探測能力及方位估計精度，制約聲吶系統(tǒng)整體性能的發(fā)揮。因此，有效估計陣形是提高機動條件下拖線陣探測性能的關(guān)鍵。

國內(nèi)外學(xué)者提出了許多陣形估計方法，但應(yīng)用于工程實際時往往存在困難。聲學(xué)計算類陣形估計方法[1]無法回避的問題是陣列采樣數(shù)據(jù)的信噪比和快拍數(shù)難以滿足要求，并且不具備拖纜段陣形估計能力?；陉囍星度氲淖藨B(tài)傳感器的陣形估計方法[2]原理相對簡單，但估計效果受到傳感器數(shù)量和精度的限制，在硬件實施上代價較大?；诹黧w力學(xué)原理的陣形估計方法[3]不依賴外部聲學(xué)環(huán)境，不需要在硬件上對基陣進(jìn)行特殊設(shè)計，可以估計包括拖纜段和聲陣段的全陣流形。Ablow模型[4]是典型的流體力學(xué)類陣形估計方法，其缺陷在于估計頻率的選取、空間網(wǎng)格的劃分、較復(fù)雜的拖船機動均可能造成雅克比矩陣求逆困難，導(dǎo)致求解發(fā)散。

針對經(jīng)典流體力學(xué)類方法穩(wěn)定性和可靠性欠佳問題，本文基于流體力學(xué)原理，分析轉(zhuǎn)向機動過程中拖線陣上各點的運動特性，探究拖線陣上相鄰兩點沿陣切線方向的變化規(guī)律，建立一種工程應(yīng)用前景較好的陣形估計方法。

1 穩(wěn)態(tài)振蕩響應(yīng)特性計算

1.1 零浮力纜情況

Paidoussis[3]在慣性坐標(biāo)系下分析柔性細(xì)長圓柱體微元段的受力平衡情況，建立了流體中零浮力纜的運動方程(Paidoussis方程)為

式(1)中，m是單位長度的纜質(zhì)量，M是單位長度纜等體積的流體質(zhì)量，dc是纜直徑，L是纜長，U是纜軸向水流速度，Ct和Cn分別是纜的切向和法向阻力系數(shù)，C′t是纜尾部的形阻系數(shù)，當(dāng)尾部處于自由狀態(tài)時該系數(shù)為零。對式(1)進(jìn)行無因次處理可得

其中，

τ=tU/L是無因次時間變量；ξ=x/L是無因次纜上位置，ξ= 0 表示拖點，ξ= 1 表示纜尾；η(ξ,τ)是纜上位置和時間的函數(shù)，表示ξ位置在τ時刻的位移；?可為任意數(shù)。

拖船做簡諧運動且纜達(dá)到穩(wěn)態(tài)時，纜上各點均做相同頻率、不同振幅的簡諧運動。在式(2)中帶入拖船位移方程，得到零浮力纜的穩(wěn)態(tài)振蕩響應(yīng)公式[5]：

其中，

Jp是p階貝塞爾函數(shù)；λ是拖船簡諧運動的波長；ω是無因次頻率；|v(z)|表示歸一化振幅，tan?1[Im(v(z))/Re(v(z))]表示相對于拖船的相位差。

拖船回轉(zhuǎn)運動情況下，達(dá)到穩(wěn)態(tài)的纜上各點做相同圓頻率、不同半徑的回轉(zhuǎn)運動，且拖船回轉(zhuǎn)運動與簡諧運動的拖船位移方程相同。因此當(dāng)纜滿足拖纜坐標(biāo)系條件(螺旋狀纜首尾相位差小于π)時，穩(wěn)態(tài)振蕩響應(yīng)公式也可用于拖船回轉(zhuǎn)機動的情況，此時無因次頻率ω= πL/(2R)，R是回轉(zhuǎn)半徑。根據(jù)纜上各點歸一化振幅和相對于拖船的相位差，即可得到拖船回轉(zhuǎn)運動時纜的穩(wěn)態(tài)流形。

1.2 拖線陣情況

式(3)只適用于單一物理屬性的纜，拖線陣的拖纜段和聲陣段的物理屬性不同，計算拖線陣穩(wěn)態(tài)振蕩響應(yīng)特性時需要分別計算拖纜段和聲陣段的振蕩特性，再將二者拼接起來。

計算拖纜段穩(wěn)態(tài)特性時，應(yīng)采用文獻(xiàn)[5]方法，首先利用式(4)對切向及法向阻力系數(shù)進(jìn)行調(diào)整，其次通過式(5)所示的單位長度纜切向阻力公式計算直行狀態(tài)下拖纜段切向阻力以及聲陣段切向阻力，進(jìn)而由式(6)計算得到聲陣段形阻系數(shù)Ct′。將調(diào)整后的阻力系數(shù)Cn?、Ct?以及形阻系數(shù)Ct′帶入式(3)即可得到拖纜段振蕩響應(yīng)特性。

式(4)中，?和θ分別是纜與流體之間的垂直夾角和水平夾角。纜與流體的水平夾角是時變的，且與無因次頻率、纜上位置相關(guān)，難以從理論角度分析。海試數(shù)據(jù)表明，對于不同無因次頻率的簡諧運動陣上各點平均水平夾角在1.1?～18.9?區(qū)間內(nèi)[5]，本文在計算時假設(shè)θ= 10?。垂直夾角可由放纜長度和嵌入纜中的深度傳感器數(shù)值計算得到。

計算聲陣段穩(wěn)態(tài)特性時，需要將拖纜段等效為一定長度的聲陣，將聲陣段和拖纜等效段作為一體進(jìn)行計算，得到整體的響應(yīng)特性再截取聲陣段部分即可，等效原則是對于某一波長的振蕩，拖纜等效段尾端與拖纜段尾端的振幅相等，即保證聲陣段首端的運動軌跡一致。

2 轉(zhuǎn)向機動中的拖線陣運動規(guī)律

實際偵測過程中，常用的戰(zhàn)術(shù)機動模式是在直行的基礎(chǔ)上調(diào)整操舵角度，令拖船以一固定轉(zhuǎn)彎半徑進(jìn)行轉(zhuǎn)向，調(diào)整到指定航向后將舵角歸零繼續(xù)直航，完整航跡由圓弧和圓弧兩端的切線組成，本文將此稱為轉(zhuǎn)向機動。

對于轉(zhuǎn)彎角度較大的轉(zhuǎn)向機動，拖船轉(zhuǎn)向后拖線陣上各點依次受到影響，脫離直行穩(wěn)態(tài)并逐漸進(jìn)入回轉(zhuǎn)穩(wěn)態(tài)；拖船結(jié)束轉(zhuǎn)向開始直行后，陣上各點依次脫離回轉(zhuǎn)穩(wěn)態(tài)并逐漸變?yōu)橹毙蟹€(wěn)態(tài)。因此可將拖船轉(zhuǎn)向機動過程中拖線陣上各點運動狀態(tài)變化情況歸納為“直行穩(wěn)態(tài)-過渡態(tài)-回轉(zhuǎn)穩(wěn)態(tài)-過渡態(tài)-直行穩(wěn)態(tài)”。陣上各點在回轉(zhuǎn)穩(wěn)態(tài)的運動特性可利用第1 節(jié)方法計算得到。以陣尾端為例，各運動階段示意圖如圖1所示。

圖1 拖線陣運動狀態(tài)示意圖Fig.1 Motion state of towed linear array

圖1中尾端在A點脫離直行穩(wěn)態(tài)，在B點進(jìn)入回轉(zhuǎn)穩(wěn)態(tài)，在C點達(dá)到圓弧頂點，在D點脫離回轉(zhuǎn)穩(wěn)態(tài)，在E點進(jìn)入直行穩(wěn)態(tài)。用TAB表示從A點運動到B點所需時間，由陣形仿真結(jié)果可知，TAC與TCE基本相等，TBC與TCD基本相等，因此可假設(shè)脫離/進(jìn)入直行穩(wěn)態(tài)和回轉(zhuǎn)穩(wěn)態(tài)的時間關(guān)于圓弧頂點對稱，仿真算例將在第4節(jié)給出。

基于上述假設(shè)，若能求得B、D兩點中任一點和A、E兩點中任一點時刻，即可得到各運動階段的分界時刻。

3 拖線陣上各點響應(yīng)時間計算

拖船轉(zhuǎn)向機動的過程中，陣上各點受到拖船轉(zhuǎn)向的影響時脫離直行穩(wěn)態(tài)，受到拖船直行的影響時脫離回轉(zhuǎn)穩(wěn)態(tài)，若能計算陣上各點對于拖船轉(zhuǎn)向和直行的響應(yīng)時間，即可確定圖1中A點和D點。

對于式(2)所示無因次Paidoussis 方程， 考慮纜尾處于自由狀態(tài)即C′t= 0 的情況，ξ ∈[0,1?dc/(2CtL))時a+bξ <0，β2?(a+bξ)>0，此時式(2)屬于雙曲型二階偏微分方程。常規(guī)拖線陣滿足dc/(2CtL)?1，因此可以近似認(rèn)為Paidoussis方程是雙曲型的。

對于雙曲型偏微分方程，如果初始時刻在某區(qū)間存在擾動，則經(jīng)過一段時間后擾動影響的區(qū)域由方程的兩條特征線決定[6]，初始擾動區(qū)間發(fā)生變化時將特征線沿x軸平移即可，示意圖如圖2所示。

圖2 擾動影響區(qū)域示意圖Fig.2 Affected area of the disturbance

拖線陣運動過程中，拖船擾動的影響會由拖點逐漸傳至纜尾，以拖船擾動時刻為初始時刻，將特征線沿x軸平移至與原點相交，則可用特征線表示陣上各點受到拖船擾動影響的時間。

3.1 近似拖直情況

計算拖線陣近似拖直情況下陣上各點對于拖船機動的響應(yīng)時間，需要分別計算聲陣段對于聲陣段首端機動的響應(yīng)時間以及拖纜段對于拖船機動的響應(yīng)時間。

聲陣段是零浮力纜，可直接利用Paidoussis 方程。式(2)的特征方程為

解此特征方程可得到描述特征線的微分方程如式(8)所示，由特征線即可得到聲陣段上各點對聲陣段首端機動的響應(yīng)時間。

拖纜段密度通常大于水密度且拖纜段尾端并非自由狀態(tài)，利用1.2 節(jié)方法將調(diào)整后的阻力系數(shù)以及形阻系數(shù)帶入式(2)，進(jìn)而求解特征線，可得到拖纜段對于拖點機動的響應(yīng)時間。

3.2 陣形畸變情況

拖船停止轉(zhuǎn)向開始直行時，拖線陣處于畸變狀態(tài)，陣與流體在水平面內(nèi)的夾角較大。而Paidoussis方程中表示微元段法向阻力的項(式(1)中最后一項)經(jīng)過線性化近似，只適用于陣流夾角較小的情況，Rispin 通過實驗證明陣流夾角大于3?時線性化會帶來較大誤差[5]。因此，陣形畸變情況下拖線陣上各點對于拖船直行的響應(yīng)時間，需要另行分析。

零浮力纜法向阻力表達(dá)式[7]為

其中，CDp是法向壓差系數(shù)，Cf是法向摩擦系數(shù)，θ表示水平方向陣流夾角。Lopes 等[8]將微元段纜受力平衡方程從歐拉坐標(biāo)系變換到拉格朗日坐標(biāo)系，展開sinθ和θ發(fā)現(xiàn)二者小于2 階的項相等，因此可將式(9)中的sinθ替換為θ。

以拖船停止轉(zhuǎn)向開始直行時為初始時刻，用θ(x,t)表示陣上x位置在t時刻的陣流夾角。已知處于自由狀態(tài)的纜尾端陣流夾角為0，拖線陣處于回轉(zhuǎn)穩(wěn)態(tài)時陣流夾角從拖點到纜尾逐漸減小，假設(shè)陣流夾角沿拖線陣線性變化，則初始時刻陣上各點的陣流夾角為

拖船停止轉(zhuǎn)向開始直行后，假設(shè)陣上各點陣流夾角隨時間線性減小，若初始時刻拖點處陣流夾角為θ0，利用全增量公式結(jié)合式(10)可得

由于?θ/?t ?1，且θ=?y/?x+U?1?y/?t[9]，式(12)可近似表示為

將式(11)和式(13)帶入式(9)，進(jìn)而用其替換式(1)中最后一項，可得處于回轉(zhuǎn)穩(wěn)態(tài)的纜在拖直過程中的運動方程，求其特征方程和特征線，即可得到聲陣段對于聲陣段首端直行的響應(yīng)時間。

計算畸變的拖纜段對于拖船直行的響應(yīng)時間時，拖纜段尾端非自由狀態(tài)導(dǎo)致式(10)不成立，假設(shè)回轉(zhuǎn)穩(wěn)態(tài)的拖纜段尾端陣流夾角為θ1，應(yīng)將初始時刻拖纜段上各點陣流夾角表示為式(14)，其中θ0和θ1可利用拖線陣回轉(zhuǎn)穩(wěn)態(tài)特性計算得到：

4 實時陣形估計方法

首先探究陣上各點過渡態(tài)時長的計算方法，以陣尾端為例，拖船轉(zhuǎn)向機動過程中拖船和陣尾端航跡示意圖如圖3所示。拖船運動到B點時開始轉(zhuǎn)向，在E點結(jié)束轉(zhuǎn)向開始直行，尾端在A點脫離直行穩(wěn)態(tài)，C點進(jìn)入回轉(zhuǎn)穩(wěn)態(tài)，D點脫離回轉(zhuǎn)穩(wěn)態(tài)，F(xiàn)點進(jìn)入直行穩(wěn)態(tài)，拖船在E點結(jié)束轉(zhuǎn)向時尾端位于H點，O點為圓心。拖船經(jīng)過某點的時間用T加上標(biāo)表示(如TB)，尾端經(jīng)過某點的時間用T加下標(biāo)表示(如TA)。

圖3 拖船和陣尾端航跡示意圖Fig.3 Track of tug and array end

設(shè)拖船轉(zhuǎn)向角度∠BOE=φ，角速度為μ。以拖船開始轉(zhuǎn)向的時刻為起始時刻，則TB= 0，TE=TH=φ/μ。

通過第3 節(jié)方法計算得到尾端對于拖船轉(zhuǎn)向和直行的響應(yīng)時間分別為tf和tl，則TA=tf，TD=TH+tl=φ/μ+tl。由于尾端在圓弧CD段角速度與拖船在圓弧BE段角速度相等，因此∠HOD=tlμ。

回轉(zhuǎn)穩(wěn)態(tài)時的拖線陣首尾相位差可用穩(wěn)態(tài)振蕩響應(yīng)模型計算得到，設(shè)∠EOH=χ。尾端在回轉(zhuǎn)穩(wěn)態(tài)階段轉(zhuǎn)過的角度為

進(jìn)而可得尾端在過渡段的時長

拖線陣處于直行穩(wěn)態(tài)時陣上各點沿陣切線方向一致，處于回轉(zhuǎn)穩(wěn)態(tài)時陣上相鄰兩點沿陣方向的差值恒定，可利用穩(wěn)態(tài)振蕩響應(yīng)模型計算得到。假設(shè)陣上相鄰兩點的沿陣方向差在過渡段線性變化，利用陣上各點在過渡段的時間以及回轉(zhuǎn)穩(wěn)態(tài)時的方向差，即可得到過渡段陣上相鄰兩點沿陣方向差的變化速率。根據(jù)之前假設(shè)，處于過渡態(tài)的相鄰兩點沿陣方向差增大和減小的速度相同。

因此，根據(jù)當(dāng)前拖船運動狀態(tài)和陣上各點對于拖船機動的響應(yīng)時間，可判斷陣上各點所處運動階段。再利用沿陣方向差變化速度計算相鄰兩點沿陣方向差。最后假設(shè)陣上相鄰兩點之間纜呈直線，利用沿陣方向差和相鄰兩點間距遞推計算陣上各點位置坐標(biāo)，實現(xiàn)陣形的實時估計。

至此，可將陣形估計算法流程歸結(jié)如下(流程圖如圖4所示)。第一步，根據(jù)拖船轉(zhuǎn)彎半徑和拖線陣物理屬性，利用式(3)計算拖線陣的穩(wěn)態(tài)回轉(zhuǎn)陣形。第二步，利用式(8)計算陣上各點對于拖點轉(zhuǎn)向和直行的響應(yīng)時間。第三步，利用式(16)計算過渡段時長，結(jié)合穩(wěn)態(tài)回轉(zhuǎn)陣形，得到陣上相鄰兩點方向差在過渡段的變化速率。第四步，根據(jù)當(dāng)前拖船運動狀態(tài)，計算陣上相鄰兩點沿陣方向差，得到陣列位置坐標(biāo)。由于避免了傳統(tǒng)流體力學(xué)類陣形估計方法中雅克比矩陣求逆等步驟，因而所提方法具有更好的穩(wěn)定性。

圖4 陣形估計算法流程圖Fig.4 Algorithm flowchart of array shape estimation

由式(16)可知，過渡段時間長度與拖船轉(zhuǎn)向角度無關(guān)，因此陣上各點對拖船轉(zhuǎn)向/直行的響應(yīng)時間和沿陣方向差變化速度僅由拖船航速和轉(zhuǎn)彎半徑?jīng)Q定。對于轉(zhuǎn)彎角度較小的情況，靠近陣尾端的點尚未進(jìn)入回轉(zhuǎn)穩(wěn)態(tài)就受到拖船直行的影響向直行穩(wěn)態(tài)過渡，可認(rèn)為這些點受到拖船直行的影響后，沿陣方向差停止增大，并以相同的變化速度逐漸減小。

5 實驗結(jié)果與分析

5.1 仿真驗證

假設(shè)拖線陣中的拖纜段和聲陣段均為光滑圓柱體，由Ansys R16.0 軟件計算其阻力系數(shù)如表1所示。

表1 光滑圓柱體阻力系數(shù)Table1 Drag coefficients of smooth cylinder

Ablow 模型經(jīng)過海試驗證，具有一定的可靠性，可以將其計算結(jié)果作為真實值。采用表1、表2所示拖線陣參數(shù)，令拖船以一定轉(zhuǎn)彎半徑轉(zhuǎn)過150?，用Ablow 模型計算轉(zhuǎn)向機動過程中的拖線陣陣形。將陣上某點的估計位置按時間順序連接可得到此點航跡，進(jìn)而可計算航跡的曲率。繪圖顯示拖點、全陣1/4位置、全陣1/2位置、全陣3/4位置、纜尾的航跡曲率，如圖5所示。

表2 拖線陣參數(shù)Table2 Parameters of towed linear array

假設(shè)某點航跡的曲率小于5×10?5時，此點處于直行穩(wěn)態(tài)，大于此點航跡曲率最大值的90%時處于回轉(zhuǎn)穩(wěn)態(tài)。依次計算陣上每個點脫離/進(jìn)入直行穩(wěn)態(tài)與達(dá)到圓弧頂點時間差的比值(圖1中TAC/TCE)，并計算此點進(jìn)入/脫離回轉(zhuǎn)穩(wěn)態(tài)與達(dá)到圓弧頂點時間差的比值(圖1中TBC/TCD)，結(jié)果如圖6所示。

圖5 陣上若干點航跡的曲率Fig.5 Curvature of the track of several points on array

圖6 各階段起止時刻沿圓弧頂點比例Fig.6 Ratio of the start and end time of each stage to the apex

由圖6可知，時間差比值基本上在1 附近，可以近似認(rèn)為脫離/進(jìn)入直行穩(wěn)態(tài)和回轉(zhuǎn)穩(wěn)態(tài)的時刻關(guān)于圓弧頂點對稱。

繪圖顯示拖船轉(zhuǎn)向后以及拖船直行后陣上各點航跡曲率變化情況，并用第2 節(jié)所述方法分別計算拖線陣上各點對于拖船機動的響應(yīng)時間理論值，結(jié)果如圖7、圖8所示，圖中白色虛線為響應(yīng)時間理論值。

由圖7、圖8分析可知，陣上各點受到拖船轉(zhuǎn)向影響后航跡曲率逐漸增大，受到拖船直行影響后航跡曲率逐漸減小，特征線方法計算得到的拖線陣響應(yīng)時間與Ablow 模型得到的航跡曲率開始變化的時間基本一致，表明第2節(jié)提出的方法較為有效。

拖線陣聲吶通過波束形成得到目標(biāo)相對于聲陣段中心的舷角，而聲吶最終輸出的是目標(biāo)相對于本船的方位角。需利用本船航向及估計的聲陣段首尾陣元連線方向，將目標(biāo)相對于聲陣段的舷角轉(zhuǎn)化為相對于本船的舷角，再轉(zhuǎn)化為相對于本船的方位角。因此，估計陣形與真實陣形各自首尾陣元連線的法向之差，可以表示陣形估計誤差造成的目標(biāo)方位估計誤差，示意圖如圖9所示。

圖7 拖船轉(zhuǎn)向后拖線陣航跡曲率變化情況Fig.7 Curvature change of the array after the tug turn

圖8 拖船直行后拖線陣航跡曲率變化情況Fig.8 Curvature change of the array after the tug goes straight

圖9 方位估計誤差示意圖Fig.9 Schematic diagram of DOA error

采用表1、表2所示拖線陣參數(shù)，假設(shè)聲陣段上均勻嵌入41 個陣元，拖船航速6 kn，令拖船以一定轉(zhuǎn)彎半徑轉(zhuǎn)過120?，以Ablow 模型計算結(jié)果為真實陣形，利用本文方法得到估計陣形。對于不同轉(zhuǎn)彎半徑的機動情況，相比于未校正陣形(直線陣)，計算陣形校正后陣列輸出信噪比的增大和方位估計誤差的減小，結(jié)果如圖10所示。

圖10 信噪比和方位估計誤差的變化情況Fig.10 Changes of SNR and DOA error

由圖10 可知，對于上述假設(shè)，當(dāng)轉(zhuǎn)彎半徑在300 m～600 m 區(qū)間內(nèi)時，利用本文方法校正陣形可使陣列輸出信噪比增大約1.57 dB，方位估計誤差減小約35.6?。

5.2 海試數(shù)據(jù)運用

海試數(shù)據(jù)來源于2017年8月在我國東中國海海域?qū)嵤┑囊淮尉C合性水聲試驗，試驗中拖船航速大約5.1 kn，全陣長1200 m，在130 s 時間內(nèi)完成一次約30?的轉(zhuǎn)向機動。不做陣形校正以及用本文方法校正陣形后的方位歷程圖如圖11、圖12所示。

圖11 不做陣形估計的方位歷程圖Fig.11 Bearing time record without array shape estimation

圖12 用本文方法校正陣形的方位歷程圖Fig.12 Bearing time record using proposed method

圖11 中不做陣形校正的左右舷方位歷程圖互為鏡像，從中選取兩個目標(biāo)進(jìn)行分析，目標(biāo)1 的舷角從88?變化到75?，目標(biāo)2 的舷角從146?變化到125?。由圖12 分析可知，用本文方法進(jìn)行陣形估計后，兩個目標(biāo)在右舷歷程圖上的軌線變得聚集且清晰，可以判斷目標(biāo)均位于右舷。以第360 s 為例，分析陣形校正前后目標(biāo)信噪比的變化，此時刻右舷功率譜圖如圖13所示。

由圖13 可知，通過陣形校正，可使目標(biāo)1 和目標(biāo)2 的輸出信噪比分別提升約1.82 dB 和1.46 dB，說明估計陣形與真實陣形匹配較好，間接證明了本文所提算法的有效性和可行性。

圖13 功率譜圖Fig.13 Power spectrum

6 結(jié)論

本文分析了轉(zhuǎn)向機動過程中拖線陣上各點的運動特性，探究拖線陣上相鄰兩點沿陣切線方向的變化規(guī)律，根據(jù)當(dāng)前拖線陣上各點的運動狀態(tài)計算相鄰兩點沿陣方向差，實現(xiàn)由局部到整體的陣形估計。計算機仿真和海上實驗數(shù)據(jù)驗證表明算法基本有效，與傳統(tǒng)的流體力學(xué)類陣形估計方法相比具有更高的穩(wěn)定性和更好的工程應(yīng)用前景。