利用拖線陣運動特性的陣形估計方法

奚 暢 蔡志明 袁 駿

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

0 引言

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

國內外學者提出了許多陣形估計方法，但應用于工程實際時往往存在困難。聲學計算類陣形估計方法[1]無法回避的問題是陣列采樣數據的信噪比和快拍數難以滿足要求，并且不具備拖纜段陣形估計能力。基于陣中嵌入的姿態傳感器的陣形估計方法[2]原理相對簡單，但估計效果受到傳感器數量和精度的限制，在硬件實施上代價較大。基于流體力學原理的陣形估計方法[3]不依賴外部聲學環境，不需要在硬件上對基陣進行特殊設計，可以估計包括拖纜段和聲陣段的全陣流形。Ablow模型[4]是典型的流體力學類陣形估計方法，其缺陷在于估計頻率的選取、空間網格的劃分、較復雜的拖船機動均可能造成雅克比矩陣求逆困難，導致求解發散。

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

1 穩態振蕩響應特性計算

1.1 零浮力纜情況

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

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

其中，

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

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

其中，

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

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

1.2 拖線陣情況

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

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

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

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

2 轉向機動中的拖線陣運動規律

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

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

圖1 拖線陣運動狀態示意圖Fig.1 Motion state of towed linear array

圖1中尾端在A點脫離直行穩態，在B點進入回轉穩態，在C點達到圓弧頂點，在D點脫離回轉穩態，在E點進入直行穩態。用TAB表示從A點運動到B點所需時間，由陣形仿真結果可知，TAC與TCE基本相等，TBC與TCD基本相等，因此可假設脫離/進入直行穩態和回轉穩態的時間關于圓弧頂點對稱，仿真算例將在第4節給出。

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

3 拖線陣上各點響應時間計算

拖船轉向機動的過程中，陣上各點受到拖船轉向的影響時脫離直行穩態，受到拖船直行的影響時脫離回轉穩態，若能計算陣上各點對于拖船轉向和直行的響應時間，即可確定圖1中A點和D點。

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

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

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

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

3.1 近似拖直情況

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

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

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

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

3.2 陣形畸變情況

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

零浮力纜法向阻力表達式[7]為

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

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

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

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

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

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

4 實時陣形估計方法

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

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

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

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

回轉穩態時的拖線陣首尾相位差可用穩態振蕩響應模型計算得到，設∠EOH=χ。尾端在回轉穩態階段轉過的角度為

進而可得尾端在過渡段的時長

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

因此，根據當前拖船運動狀態和陣上各點對于拖船機動的響應時間，可判斷陣上各點所處運動階段。再利用沿陣方向差變化速度計算相鄰兩點沿陣方向差。最后假設陣上相鄰兩點之間纜呈直線，利用沿陣方向差和相鄰兩點間距遞推計算陣上各點位置坐標，實現陣形的實時估計。

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

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

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

5 實驗結果與分析

5.1 仿真驗證

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

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

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

表2 拖線陣參數Table2 Parameters of towed linear array

假設某點航跡的曲率小于5×10?5時，此點處于直行穩態，大于此點航跡曲率最大值的90%時處于回轉穩態。依次計算陣上每個點脫離/進入直行穩態與達到圓弧頂點時間差的比值(圖1中TAC/TCE)，并計算此點進入/脫離回轉穩態與達到圓弧頂點時間差的比值(圖1中TBC/TCD)，結果如圖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 附近，可以近似認為脫離/進入直行穩態和回轉穩態的時刻關于圓弧頂點對稱。

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

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

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

圖7 拖船轉向后拖線陣航跡曲率變化情況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所示拖線陣參數，假設聲陣段上均勻嵌入41 個陣元，拖船航速6 kn，令拖船以一定轉彎半徑轉過120?，以Ablow 模型計算結果為真實陣形，利用本文方法得到估計陣形。對于不同轉彎半徑的機動情況，相比于未校正陣形(直線陣)，計算陣形校正后陣列輸出信噪比的增大和方位估計誤差的減小，結果如圖10所示。

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

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

5.2 海試數據運用

海試數據來源于2017年8月在我國東中國海海域實施的一次綜合性水聲試驗，試驗中拖船航速大約5.1 kn，全陣長1200 m，在130 s 時間內完成一次約30?的轉向機動。不做陣形校正以及用本文方法校正陣形后的方位歷程圖如圖11、圖12所示。

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

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

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

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

圖13 功率譜圖Fig.13 Power spectrum

6 結論

本文分析了轉向機動過程中拖線陣上各點的運動特性，探究拖線陣上相鄰兩點沿陣切線方向的變化規律，根據當前拖線陣上各點的運動狀態計算相鄰兩點沿陣方向差，實現由局部到整體的陣形估計。計算機仿真和海上實驗數據驗證表明算法基本有效，與傳統的流體力學類陣形估計方法相比具有更高的穩定性和更好的工程應用前景。