姜潤翔 胡英娣 龔沈光

（1.海軍工程大學兵器工程系，湖北 武漢430033；2.中國人民解放軍91872部隊，北京100071）

引 言

為了有效地解決船舶靜電場的深度換算問題，前提是建立正確的正演模型，在此基礎上，通過實測數(shù)據(jù)建立關于場源信息的超定方程組，而后利用最優(yōu)求解算法求解場源，最后利用反演后的場源對其他條件下的電場進行建模運算［1-2］.上述問題的難點在于正演模型的建立和場源最優(yōu)解的計算（超定方程組的求解）.在船舶靜電場的建模方面，目前主要有邊界元法［3］、有限元法［4］、偶極子源法［5］和電流元法［6］，其中，等效偶極子源和電流元法因相對求解的場源信息較少而被廣泛應用于船舶電場的深度換算領域［6-8］.由于船舶靜電場可近似為電流元場的疊加［9-10］，考慮到電流元可等效為點電源，提出了基于點電源的船舶靜電場反演方法.

1 基于點電源的靜電場建模方法

在對船舶靜電場進行建模時，可將靜電場視為由電荷密度為σ（S）的閉曲面S所產(chǎn)生的，則S在某一場點P處的電位可表示為［11］

式中：γ為海水電導率；K（S，P）為點P與S之間的距離函數(shù).

將表面S分為n個小面源Si，i＝1，2，…，n，每個小面源的電荷密度為σ（Si），則船舶電場真實源在某一場點P處的電位可近似表示為n個等效面源在該點產(chǎn)生電位的和，即

式中：Qi為面源Si上的電荷；ΔSi為等效面源的面積；K（Si，P）為面源Si和場點P之間的距離函數(shù).

由電多極矩的相關知識可知，如果場點P（坐標為（xp，yp，zp））距離等效面源Si區(qū)域較遠，可將K（Si，P）相對面源Si內(nèi)的某一點（xi，yi，zi）進行泰勒展開，得

由于在0階收斂半徑之外，可忽略高次項，僅保持第一項［12］，因此，式（2）可表示為

式中，K1（Si，P）為等效點源的坐標（xi，yi，zi）到場點P處的距離函數(shù).

根據(jù)式（4）可知，船舶的靜電場可等效為n個點源產(chǎn)生的電場和.在利用式（4）進行電場計算時，其核心參數(shù)K1（Si，P）在空氣-海水-海床三層均勻介質條件下為［13］

式中：H為海水深度；k＝（γ-γ1）／（γ＋γ1）為海底反射系數(shù)，γ1為海床電導率；m為反射層數(shù)，實際計算中其上限值可取10～20［14］.

值得注意的是，在利用式（4）計算場點P處的電位時，為了減小計算誤差，場點與等效點源之間的距離應大于函數(shù)K（Si，P）泰勒公式展開時的收斂半徑.

有了海水中任意一點的電位表達式（4），便可計算出該點的三分量電場表達式為

式中：

2 基于點電源的靜電場反演方法

由第1節(jié)分析可知，船舶的靜電場可等效為n個點源產(chǎn)生的電場和，因此，在利用點源疊加的方法對靜電場進行反演時，需要計算的未知參數(shù)主要包括：點源個數(shù)n、點源坐標（xi，yi，zi）和點源等效電荷Qi.若能利用先驗信息確定點源個數(shù)n和點源坐標（xi，yi，zi），則將簡化計算量.

2.1 點源位置計算

根據(jù)理論計算和實測數(shù)據(jù)分析可知，船舶靜電場電位關于中軸線對稱，因此，在反演時可將等效源分布在沿船體縱向方向的平行線上，平行線關于軸對稱分布，其數(shù)量m可以是任意的偶數(shù)（通常取m＝2），每條平行線與船體肋骨橫截面的交點即為等效點源的位置.

某一肋骨截面等效點源的位置（m＝2）具體確定方法如下：建立如圖1所示的船體坐標系o-xyz，其中，y軸為船體的橫向方向，z軸為垂直方向，｜od｜＝B／2（B為船寬）.為了使同一截面上等效點源的位置關于船中軸線對稱，將等效點源e和f的坐標分別位于cd和cg的中垂線上是合適的，且e和f關于軸線對稱分布，為了保證所有測量點均在泰勒展開的收斂半徑之外，應保證由點源e形成并經(jīng)過c、d和最近測量點s的圓的半徑最小.利用以上方法，可方便地求出等效源點e、f的位置.需要說明的是，上述點源位置的確定方法并不是為了對點源進行準確的定位，而是一種滿足由式（4）進行等效計算條件下合理點源位置的計算方法.

圖1 兩條平行線等效源坐標示意圖

為了保證換算具有較高的精度，所有測量點應在K（S，P）泰勒公式展開時的收斂半徑之外，因此可將最近測量點s與等效點源e之間的距離R視為最大收斂半徑，在此條件下，可保證其他測量點均在收斂半徑的范圍之外.在上述基礎之上，根據(jù)圖2很容易確定點源的最小個數(shù)nmin為

式中：z0為e點的垂直坐標；T為船舶吃水深度；L為船舶水線長度為點源縱向方向的最大間隔；方括號表示選擇大于括號里數(shù)值的最小偶數(shù).

圖2 等效源間隔示意圖

實際計算過程中，為了提高精度，通常取n稍大于nmin.根據(jù)n即可方便地求出等效點源的位置信息：

若在實際計算過程中，發(fā)現(xiàn)等效點源的位置在船體橫截面積以外，將等效源分布在沿船體縱向方向的一條直線（在o-xz平面與中軸線平行）上是合適的，該直線與肋骨截面的交點即為等效點源的位置.與兩條平行線點源位置的確定方法類似，此時，只需保證點源形成并經(jīng)過d、g和最近測量點s的圓的半徑最小.為了確保某一深度平面的所有可能的測量點按照式（4）進行計算時，均具有較高的精度，可將過最近測量點s的圓的半徑最小轉變?yōu)辄c源與最近測量平面相切的問題（如圖3所示），得

式中，zh為距e最近的測量點s的深度.當h≤T時，取z0＝h，否則取z0＝T.

取R＝zh-z0，即有

圖3 一條平行線等效源坐標示意圖

根據(jù)式（10）、（11）可計算所需要的最小點源個數(shù)及每個點源的位置，由于所有點源在同一條直線上，由式（11）獲取的點源位置存在重合現(xiàn)象，但在計算時不影響，為了進一步減小誤差，在計算點源個數(shù)時，可取R1＝zh-z0

2.2 反演算法

在確定點源個數(shù)及其坐標后，可根據(jù)式（6）～（9）對測量點建立電場計算模型，假設源點個數(shù)為n，點源坐標分別為（xi，yi，zi），i＝1，2，…，n，測量點坐標及其電場值分別為（xj，yj，zj）和（Exj，Eyj，Ezj），j＝1，2，…，m.分別計算每個點源相對于每個測量點之間的距離函數(shù)K1x（Si，Pj）、K1y（Si，Pj）和K1z（Si，Pj），則可建立電場方程：

通過方程（14）的求解可得到Qi，在實際測量過程中有m＞n，方程（14）為冗余方程，此時可采用最小二乘法對方程進行求解，即式（14）轉化為最小泛函數(shù)

的求解.

為了避免超定方程的求解問題，對式（15）增加控制方程，得到

式中，α為控制參數(shù).

同時，式（16）的解應滿足電中性方程

利用最小值條件

得到n個方程.聯(lián)合式（17），得到n＋1個方程，其矩陣表達式為

定 義K1x（Si，Pj）＝ηij、K1y（Si，Pj）＝κij、K1z（Si，Pj）＝λij，則式（19）中的aii、aij、bi為

綜上所述，在獲得點源個數(shù)n、點源坐標（xj，yj，zj）和點源等效電荷Qi之后，便可對其偶極矩進行估算：

同時，可以利用式（6）對任意場點的電場信號進行計算.

3 典型算例

3.1 仿真船舶電場數(shù)據(jù)

為了檢驗上述方法的有效性，利用英國研制的Beasy軟件（該軟件在船舶設計階段被用于對船舶電場信號進行評估［15］）對某船舶進行電場計算，其中，船舶水線長L＝117m，水線寬B＝15.6m，吃水深度T＝6.27m.船殼的材料為鋼，涂有防腐蝕涂層，螺旋槳材料為銅鎳鋁合金，假設該船舶船身涂層大致完好，但螺旋槳是裸露在海水中的，且船殼某處涂層有破損，螺旋槳和船殼破損區(qū)域即是被腐蝕區(qū)域.船身上有兩對位置對稱的輔助陽極，四個輔助陽極的坐標分別為（10.5，-4.9，0.8），（10.5，4.9，0.8），（55，-7，3.69），（55，7，3.69）.船身破損區(qū)域和陽極位置如圖4中標注所示.三維的船殼模型被離散為1 708個單元，船身的大部分單元為四邊形，在船首船尾處有少量單元為三角形，每個分割成的單元可以表示該部分的電壓、電流和電流密度等各種信息.

圖4 離散后的船殼模型

設海水電導率γ＝4S／m、海底電導率γ1＝0.01S／m，海水深度H＝500m.令兩對輔助陽極分別輸出-12A和-8A的保護電流，水面以下10 m和20m處的靜電場Ex、Ey、Ez和電位分別如圖5和圖6所示.

從圖5中選取200個點作為已知場點，其坐標如表1所示，電場值如圖7所示.

表1 場點坐標

3.2 反演及換算性能檢驗

根據(jù)2.1節(jié)中的點源坐標計算方法，可求出e點的位置y0＝-8.10m、z0＝8.343m，y0和z0均在船體橫截面積之外.因此，可取y0＝0m，利用式（12）計算得到z0＝1.958m.取R1＝zh-z0＝8.042m，根據(jù)式（10）可計算出最小點源個數(shù)為nmin＝30，本例中取n＝32，此時可求出等效點源的間隔為R1＝2L／n＝7.32m，從而等效點源的坐標值為

圖5 水深10m的電場及電位分布

圖6 水深20m的電場及電位分布布

圖7 不同正橫距下的船舶電場信號

根據(jù)式（19）、（20）對靜電場建立模型以求解Qi，值得注意的是，式（20）中的α選擇對算法的結果影響較大，計算過程中可將α依次選取為α1＝10，αs＝0.1×αs-1，s＝2，3，…，分別計算其等效電荷Qi，并計算sum（Q）＝Q1＋Q2＋…，Qn，若其值在一定范圍內(nèi)接近于0，并且穩(wěn)定，即可將α取為當前所設定的值.在此算例中，s＝10即α＝10-8時，sum（Q）接近于0，并且穩(wěn)定.求出Qi后，利用式（21）可計算出Mx＝-1140A·m、My＝0A·m和Mz＝-4.6·10-10A·m.

值得注意的是，模型反演出的點源位置和Qi應與實際船舶電荷分布相近，本例中等效電荷Qi的大小及其分布如圖8所示.從圖8可明顯發(fā)現(xiàn)：在10.97m、54.84m分別出現(xiàn)了電荷的極值，與輔助陽極的位置（10m、55m）相符，說明了模型的正確性；但在船殼破損區(qū)域處的電荷值不明顯，這是由于輔助陽極為點電極，因此在利用點電源進行反演時對其定位精度較高.若適當增加n的個數(shù)，對局部破損區(qū)域的定位精度有望得到提高.圖9是n＝44時反演出的等效電荷大小及其分布.從圖9可明顯發(fā)現(xiàn)，不但在10m和55m附近出現(xiàn)了電荷極值，在20m位置附近，同樣出現(xiàn)了電荷的極值，與破損區(qū)域的位置相近.

在求出所有的源點信息后，便可對任意場點的電場值進行計算，圖10和圖11分別為z＝10m時，y＝-10m和y＝10m的實際電場信號和換算結果.圖12為z＝20m，y＝-10m時的實際電場信號及換算結果.

圖8 n＝32時等效電荷大小及其分布

圖10 z＝10m，y＝-10m時的電場換算及實際結果

從圖10～12可明顯發(fā)現(xiàn)，利用z＝10m正橫距為y＝-6m和y＝6m的200個測量點的電場值，向同一深度（z＝10m）、不同正橫距（y＝-10 m、y＝10m）和同一正橫距（y＝-10m）、不同深度（z＝20m）的場點進行換算時，精度較高，換算結果可近似滿足實際的需要.但在輔助陽極位置附近，換算誤差明顯增大，這是由于輔助陽極位置附近的電場變化劇烈，利用式（5）對電場進行等效計算時誤差較大，為了提高計算精度，應適當增加泰勒公式展開的項數(shù).

圖11 z＝10m，y＝10m時的電場換算及實際結果

圖12 z＝20m，y＝-10m時的電場換算及實際結果

4 結 論

為了對船舶靜電場進行反演，提出了基于點電源的船舶靜電場深度換算方法，該方法首先利用先驗信息對點源的位置進行了確定，避免了以往算法中源位置未知引起的未知參數(shù)過多而導致的局部最小解問題，通過仿真數(shù)據(jù)檢驗發(fā)現(xiàn)，該算法換算效果較好，為船舶靜電場反演及換算方法提供了一種新的思路.

［1］林春生，龔沈光.艦船物理場［M］.北京：兵器工業(yè)出版社，2007：237-242.

［2］HUBBARD J C，BROOKS S H，TORRRANCE B C.Practical measures for reduction and management of the electro-magnetic signatures of in-service surface ships and submarines［C］／／Underwater Defense Technology Conference，1996：64-65.

［3］ZAMANI N G.Boundary element simulation of the cathodic protection system in a prototype ship［J］.Applied Mathematics and Computation，1988，26：119-123.

［4］DOIG P，F(xiàn)LEWITT P E J.A finite difference numerical analysis of galvanic corrosion［J］.Electrochem，1979，126（12）：2057-2063.

［5］LU Xincheng，CHEN Hunqing，MA Shoujun，et al.The electromagnetic field generated by time-harmonic electric current element in shallow sea［C］／／Undersea Defense Technology Conference.Hamburg，2006.

［6］陳 聰，龔沈光，李定國.艦船靜態(tài)電場深度換算方法［J］.哈爾濱工程大學學報，2009，30（6）：719-722.CHEN Cong，GONG Shenguang，LI Dingguo.The method of extrapolation of the static electric field of ships［J］.Journal of Harbin Engineering University，2009，30（6）：719-722.（in Chinese）

［7］劉文寶，王向軍，稽 斗.基于電偶極子模型的艦船靜電場深度換算［J］.空軍雷達學院學報，2010，24（6）：436-439.LIU Wenbao，WANG Xiangjun，JI Dou.Conversion of static electric field depth of ships based on electric dipole model［J］.Journal of Air Force Radar Academy，2010，24（6）：436-439.（in Chinese）

［8］陳 聰，李定國，龔沈光.基于拉氏方程的艦船靜態(tài)電場深度換算［J］.電子學報，2010，38（9）：2025-2029.CHEN Cong，LI Dingguo，GONG Shenguang.Research on the extrapolation of static electric field of ships based on laplace equation［J］.Acta Electronica Sinica，2010，38（9）：2025-2029.（in Chinese）

［9］DYMARKOWSKI K，UCZCIWEK J.The extremely low frequency electromagnetic signature of the electric field of the ship［C］／／Underwater Defense Technology Conference.London，2001：1-6.

［10］BOSTICK F，SMITH H，BOEHL J.The Detection of ULF-ELF Emissions from Moving Ships［R］.New York：State Academic Educational Institutions，1977：13-24.

［11］SUSANU T，LINGVAY I，STOIAN F.The finite differences method to determine the potential function values in the three-dimensional space around a seagoing ship［C／OL］.1991［2013-07-15］.http：／／www.fict.ro／～ep／publications＿files／elecship＿98＿slsp.pdf.

［12］VANDERLINDE J.Classical Electromagnetic Theory［M］.2nd ed.London：Springer，2005：33-35.

［13］JONES D L，BURKE C P.The DC filed components of horizontal and vertical electric dipole sources immersed in three-layered stratified media［J］.Annales Geophysicae，1997，15：503-510.

［14］BANNISTER P.New Formulas for HED，HMD，VED，and VMD Subsurface-to-Subsurface Propagation［R］.USA：Naval Underwater Systems Center，1984：1-38.

［15］ADEY R，BAYNHAM J.Predicting corrosion related electrical and magnetic fields using BEM［C］／／Underwater Defense Technology.Tracey Westwood，2000：386-391.