基于點電荷模型的艦船靜電場反演算法研究

姜潤翔，林春生，龔沈光

(海軍工程大學 兵器工程系，湖北 武漢430033)

0 引言

為了全面了解不同海區(qū)、不同海深條件下的艦船靜電場分布，可在部分實測數(shù)據(jù)的基礎上，對場源模型進行反演，而后利用反演后的場源對其他條件下的電場進行推算［1-2］。上述問題的關鍵是場源模型的建立，考慮到艦船靜電場可視為若干個電偶極子的疊加［3-5］，若能利用某個測量深度的實測數(shù)據(jù)換算出每個電偶極子的信息(位置、強度和方向)，即可達到反演場源模型的目的。但是偶極子的信息往往是未知的，從而導致在利用電偶極子法反演源強度時，因所需要求解的未知量較多而引起的數(shù)據(jù)量大且計算復雜的問題［3-4］。為了避開反演源強度的問題，文獻［5］提出了直接由已知場點的電場值推算其他場點電場值的深度換算方法，為艦船靜電場的反演提供了一種新的思路，但該方法需要測量某一平面的靜電場信號，受測量條件的限制，難以在實際中應用，因此需要研究其他艦船靜電場的反演方法。

1 艦船靜電場

由于艦船通常是由不同金屬材料制成的，當其在海水(電解液)中航行時，不同材質的金屬部件在形成電解偶的過程中，會產生腐蝕電流，同時，陽極會逐漸被腐蝕而發(fā)生剝落。比如在“螺旋槳—船體”電解偶系統(tǒng)(如圖1所示，A 為陽極)中，陽極為鋼制船殼(受到腐蝕)，陰極為銅制螺旋槳。為了保護艦船殼體不受腐蝕，除采用涂層防腐外，現(xiàn)代艦船普遍采用了外加電流陰極保護(ICCP)系統(tǒng)產生保護電流進行防腐。腐蝕電流和防腐電流在艦船周圍產生的電場信號統(tǒng)稱為艦船靜電場信號［6］。

圖1 螺旋槳-船體電化學系統(tǒng)示意圖Fig.1 Electrochemical system formed by propeller and hull

對于一艘涂層相對完好的艦船，其電解偶電路主要有“螺旋槳—船殼”、“通海閥—船殼”、“聲納導流罩—船殼”等幾種形式，它們均是沿艦船縱向分布的，因此艦船靜電場可視為沿艦船縱向分布的偶極子或電流元疊加而成［7-9］，考慮到電流元可等效為電荷，進而靜電場可等效為一系列點電荷產生的電場疊加。

2 靜電場的點電荷建模方法

在對艦船靜電場進行建模時，假設水線以下船殼(閉曲面S)上的電荷密度為σ(S)，則S 在某一場點P 處的電位［10］可表示為

式中:γ 為海水電導率;K(S，P)為點P 與S 之間的距離函數(shù)。

將表面S 分為n 個小面源Si，i =1，2，…，n，每個小面源的電荷密度為σ(Si)，則電位U 可視為n 個小面源在場點P 處產生電位值的疊加［11］，即有

式中:Qi為面源Si上的電荷大小;ΔSi為等效面源的面積;K(ΔSi，P)為面源Si和場點P 之間的距離函數(shù)。

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

在0 階收斂半徑之外，忽略高次項，僅保持第一項，則(2)式可表示［12］為

式中:K(Si，P)為等效點電荷的坐標(xi，yi，zi)到場點P(xP，yP，zP)處的距離函數(shù)。

根據(jù)鏡像法理論可知，計算海水中單位點電荷在3 層(空氣—海水—海床)均勻介質條件下海水中任意場點的電位分布時，其距離函數(shù)為

式中:H 為海水深度;k =(γ -γ1)/(γ +γ1)為海底反射系數(shù)，γ1為海床電導率;m 為反射層數(shù)，實際計算中其上限值可取10 ～20［13］.

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

3 基于點電荷模型的靜電場反演方法

由第2 節(jié)分析可知，在利用點電荷疊加的方法對靜電場進行反演時，需要計算的未知參數(shù)主要包括:點電荷個數(shù)n、每個點電荷的坐標(xi，yi，zi)和Qi，未知參量較多，增加了場源反演模型計算的復雜性。若能利用先驗信息確定n 的大小和每個點電荷的坐標(xi，yi，zi)，則將減小待求解未知數(shù)的數(shù)量。

3.1 點電荷個數(shù)及坐標的確定

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

m=2 時某一肋骨截面等效點電荷的位置確定方法如下:建立如圖2所示的船體坐標系Oxyz，其中，y 軸為船體的橫向方向，z 軸為垂直方向，Od =B/2(B 為船寬)。為了使同一截面上等效點電荷的位置關于船體中軸線對稱，將等效點電荷e 和f 的坐標分別位于cd 和cg 的中垂線上是合適的(e 和f關于軸線對稱分布)，同時，為了保證所有測量點均在泰勒展開的收斂半徑之外，應保證由點電荷e 形成并經過c、d 和最近測量點s 的圓的半徑最小。利用上述方法，可計算出等效源點電荷e、f 的位置yi和zi.

圖2 等效點電荷橫向位置示意圖Fig.2 Schematic diagram of position of equivalent point charge

需要說明的是，上述點電荷位置的確定方法并不是為了對點電荷進行準確的定位，而是在滿足由(4)式進行等效計算時的一種合理點源位置計算方法。

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

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

圖3 等效點電荷縱向間隔示意圖Fig.3 Schematic diagram of longitudinal separation of equivalent point charge

若在實際計算過程中，發(fā)現(xiàn)點電荷e 的坐標yi和zi在船體橫截面積以外，從近似的角度來看，將等效源的位置放置在船的中軸線上是合適的，為了保證圓的半徑最小，應保證圓與測量平面相切，則根據(jù)點電荷e 到點d 和測量平面的距離相等，有

式中:zh為距e 最近的測量點s 的深度。當h≤T時，取z0=h，否則取z0=T，此時可取即有

實際計算過程中，為了提高精度，通常取n 稍大于nmin. 根據(jù)n 即可方便地求出等效點源沿船體縱向方向的坐標:

3.2 反演算法

假定實際測量點的坐標及其電位值分別為(xj，yj，zj)、Uj，j=1，2，…，m，根據(jù)(5)式可計算每個點源(xi，yi，zi)，i=1，2，…，n 和每個測量點之間的距離函數(shù)Kij(Si，Pj)，結合(4)式，有

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

的求解。

對(11)式增加控制方程:

式中:α 為控制參數(shù)。

同時，(12)式的解還應滿足電中性方程:

利用最小值條件:

結合(13)式，可得到(n+1)×(n+1)階矩陣

式中:

對比(10)式和(15)式可知，原方程m ×n 的矩陣轉化為(n +1)×(n +1)的方陣，由于約束方程(13)式的引入，減小了Qi解的范圍，進一步簡化了運算。

由(15)式可計算出Qi，進而根據(jù)(5)式對任意場點建立正演模型，達到反演任意場點電位的目的。同時，在獲得點源個數(shù)和點源等效電荷之后，還可對艦船的等效偶極矩Mx、My和Mz進行估算，

4 典型算例

4.1 實測艦船靜電場數(shù)據(jù)

為了獲取實船周圍的電場信號，利用電場傳感器(Ag/AgCl 電極)對某小型運輸船周圍的電位信號進行了測量，其中，艦船水線長L =73.5 m，寬B =8 m，吃水深度T=5.5 m，測量區(qū)域海深H=70 m，海水電導率γ=3.0 S/m，海底電導率γ1=0.003 S/m.建立測量坐標系Oxyz(如圖4所示)，其中，船首的位置為原點O，x 方向為船的縱向方向，y 方向為船的左舷方向，z 方向垂直向下。

圖4 測量坐標系Fig.4 Measuring coordinate system

試驗中，為了獲取船周圍不同位置的電位信號，利用1 對Ag/AgCl 電極對3 個深度下(z =6.5 m，z=8.0 m和z=22.5 m)不同位置的電位信號進行測量，其中，參比電極放置在距離原點250 m 的位置(此處的電位可近似視為零電位)，將測量電極依次放置在事先設定的測量點處進行逐點測量。每種測量深度下，測量點數(shù)為102 個，分別布置在y =-4 m，y=0 m 和y=4 m 的水平線上，其測量點坐標如表1所示。

表1 測量點坐標Tab.1 The coordinates of measuring points

為了檢驗換算的效果，選取z =8.0 m 測量得到的102 點的電位信號為學習樣本(其電位如圖5所示)，將z=6.5 m 和z=22.5 m 的電位信號作為檢驗樣本，從而達到檢測算法的目的。

圖5 z=8.0 m 時的實測艦船電位信號Fig.5 Measured electric potential signal for z=8.0 m

4.2 反演結果

根據(jù)3.1 節(jié)中的點源坐標計算方法，可求出e 點的位置y0=2.6 m，z0=3.2 m，同時得到收斂半徑R=3.45 m，根據(jù)(6)式可計算出最小點源個數(shù)為nmin=30，本例中取n =36，此時可求出等效點源的間隔為Δx=2L/n=4.08 m，從而等效點源的坐標值為

根據(jù)(15)式和(16)式對靜電場建立模型以求解Qi，值得注意的是，(16)式中的α 選擇對算法的結果影響較大，可將α 依次選取為α1=10，αs=0.1 ×αs-1，s =2，3，…分別計算其等效電荷Qi，并計算sum(Q)=Q1+Q2+，…，Qn，若其值在一定范圍內接近于0，并且穩(wěn)定，即可將α 取為當前所設定的值。在此算例中，s 與sum(Q)的關系曲線如圖6所示，從圖6中可明顯發(fā)現(xiàn)在s=10 時，即α=10-8時，sum(Q)接近于0，并且穩(wěn)定，因此文中取α =10-8.待求出Qi后，利用(17)式可計算出Mx= -64.6 A·m，My= -4.53 A·m 和Mz= -2.59 ×10-10A·m.

在求出所有的點電荷的信息后，便可對任意場點的電位值進行計算，圖7和圖8分別為測量深度z=6.5 m 和z =22.5 m 的換算結果和實測結果，圖中橫坐標為測量點坐標，縱坐標為測量點處的電位值，虛線為實測結果，實線為換算結果。從圖7和圖8中可明顯發(fā)現(xiàn)，實測結果和換算結果相似度較高。

圖6 s 與sum(Q)的曲線圖Fig.6 s vs. sum(Q)

圖7 z=6.5 m 時的換算與實測結果Fig.7 Measured and extrapolated electric potential signals for z=6.5 m

圖8 z=22.5 m 時的換算與實測結果Fig.8 Measured and extrapolated electric potential signals for z=22.5 m

為了定量的衡量反演結果的精度，將反演得到的電位Uj與實際電位Vj之間的誤差定義為

從而得到不同測量線上的換算誤差值σ 如表2所示。

從表2中可清晰發(fā)現(xiàn)，水深z =22.5 m 的反演精度明顯高于z=6.5 m 的換算精度，這是由于算法本身的特點所決定的，因為z =22.5 m 上的測量點均在z=8.0 m 上測量點所建等效點源的泰勒展開收斂半徑范圍之外，而z=6.5 m 上的測量點在等效點源的收斂半徑范圍之內，因此由(3)式轉化到(5)式時誤差較大。由上述分析可知，在利用本文方法反演艦船靜電場時，為了獲得較高的精度，應盡量以近距離的電場值反演遠距離的電場值。

表2 不同測量線上的換算誤差值Tab.2 Extrapolated errors of different measuring lines

5 結論

針對艦船靜電場的反演問題，提出了一種基于點電荷法的反演算法，其本質在于艦船靜電場可等效為一系列點電荷產生的電場疊加。它與以往基于電偶極子的反演算法相比，其主要特點為:

1)點電荷為標量，而電偶極子為矢量(未知個數(shù)較多)，減少了反演時超定方程未知數(shù)的個數(shù)，場源的位置進行了預先確定，進一步減少了未知量的個數(shù)。

2)由于點電荷之和滿足電中性的條件，增加了約束方程，減小了超定方程的求空間，避免陷入局部最小解。

3)算法反演精度較高，實測數(shù)據(jù)檢驗結果表明:該算法由水深8 m 的測量數(shù)據(jù)向水深22.5 m 進行反演時，換算精度較高，最大誤差值僅為3.7%.

4)反演出的艦船電場源強度包含了船體表面電流大小的分布狀況，該信息可為艦船電場防護和隱身提供基礎數(shù)據(jù)。

References)

［1］Hubbard J C，Brooks S H，Torrrance B C. Practical measures for reduction and management of the electro-magnetic signatures of inservice surface ships and submarines［C］∥Underwater Defence Technology Conference. London，UK:Royal Navay，1996:64 -65.

［2］林春生，龔沈光. 艦船物理場［M］. 北京:兵器工業(yè)出版社，2007:237 -242.LIN Chun-sheng，GONG Shen-guang. Physical fields of ships［M］. Beijing:Publishing House of Ordnance Industry，2007:237-242. (in Chinese)

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

［4］劉文寶，王向軍，稽斗. 基于電偶極子模型的艦船靜電場深度換算［J］. 空軍雷達學院學報，2010，24(6):436 -439.LIU Wen-bao，WANG Xiang-jun，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)

［5］陳聰，李定國，龔沈光. 基于拉氏方程的艦船靜態(tài)電場深度換算［J］.電子學報，2010，38(9):2025 -2029.CHEN Cong，LI Ding-guo，GONG Shen-guang. 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)

［6］Daya Z A，Hutt D L，Richards T C. Maritime electromagnetism and DRDC management research［R］. Canada:Defence R＆D Canada-Atlantic，2005:1 -80.

［7］Dymarkowski K，Uczciwek J. The extremely low frequency electromagnetic signature of the electric field of the ship［C］∥Underwater Defence Technology Conference. London，UK:SAAB Defence and Security，2001:1 -6.

［8］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.

［9］Mcgrath J N，Tighe-Ford D J，Hodgkiss L. Scale modeling of a ship＇s impressed-current cathodic protection system［J］. Corrosion Prevention ＆ Control，1985，4:36 -39.

［10］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(4):503 -510.

［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］. 2th ed.London，UK:Springer，2005:33 -35.

［13］King R W P.The electromagnetic field of a horizontal electric dipole in the presence of a three-layered region:supplement［J］.Journal of Applied Physics，1993，74(8):4845 -4848.

［14］Ditchfield R W，Mcgrath J N，Tighe-Ford D J.Theoretical validation of the physical scale modeling of the electrical potential characteristics of marine impressed current cathodic protection［J］.Journal of Applied Electrochemistry，1995，25(1):54 -60.

［15］Degiorgi V G，Brebbia C A，Adey R A. Simulation of electrochemical processes II［C］∥Adey R，Baynham J M W. Predicting corrosion related signature:2nd International Conference on the Simulation Electochemical Processes. Southampton:WIT Press，2007:368 -371.