基于點(diǎn)電荷模型的艦船靜電場(chǎng)反演算法研究

姜潤(rùn)翔，林春生，龔沈光

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

0 引言

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

1 艦船靜電場(chǎng)

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

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

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

2 靜電場(chǎng)的點(diǎn)電荷建模方法

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

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

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

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

由電多極矩的相關(guān)知識(shí)可知，如果場(chǎng)點(diǎn)P 距離等效面源Si區(qū)域較遠(yuǎn)，可將K(ΔSi，P)相對(duì)面源Si內(nèi)的某一點(diǎn)(xi，yi，zi)進(jìn)行泰勒展開，得

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

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

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

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

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

3 基于點(diǎn)電荷模型的靜電場(chǎng)反演方法

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

3.1 點(diǎn)電荷個(gè)數(shù)及坐標(biāo)的確定

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

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

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

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

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

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

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

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

式中:zh為距e 最近的測(cè)量點(diǎn)s 的深度。當(dāng)h≤T時(shí)，取z0=h，否則取z0=T，此時(shí)可取即有

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

3.2 反演算法

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

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

的求解。

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

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

同時(shí)，(12)式的解還應(yīng)滿足電中性方程:

利用最小值條件:

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

式中:

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

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

4 典型算例

4.1 實(shí)測(cè)艦船靜電場(chǎng)數(shù)據(jù)

為了獲取實(shí)船周圍的電場(chǎng)信號(hào)，利用電場(chǎng)傳感器(Ag/AgCl 電極)對(duì)某小型運(yùn)輸船周圍的電位信號(hào)進(jìn)行了測(cè)量，其中，艦船水線長(zhǎng)L =73.5 m，寬B =8 m，吃水深度T=5.5 m，測(cè)量區(qū)域海深H=70 m，海水電導(dǎo)率γ=3.0 S/m，海底電導(dǎo)率γ1=0.003 S/m.建立測(cè)量坐標(biāo)系Oxyz(如圖4所示)，其中，船首的位置為原點(diǎn)O，x 方向?yàn)榇目v向方向，y 方向?yàn)榇淖笙戏较颍瑉 方向垂直向下。

圖4 測(cè)量坐標(biāo)系Fig.4 Measuring coordinate system

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

表1 測(cè)量點(diǎn)坐標(biāo)Tab.1 The coordinates of measuring points

為了檢驗(yàn)換算的效果，選取z =8.0 m 測(cè)量得到的102 點(diǎn)的電位信號(hào)為學(xué)習(xí)樣本(其電位如圖5所示)，將z=6.5 m 和z=22.5 m 的電位信號(hào)作為檢驗(yàn)樣本，從而達(dá)到檢測(cè)算法的目的。

圖5 z=8.0 m 時(shí)的實(shí)測(cè)艦船電位信號(hào)Fig.5 Measured electric potential signal for z=8.0 m

4.2 反演結(jié)果

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

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

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

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

圖7 z=6.5 m 時(shí)的換算與實(shí)測(cè)結(jié)果Fig.7 Measured and extrapolated electric potential signals for z=6.5 m

圖8 z=22.5 m 時(shí)的換算與實(shí)測(cè)結(jié)果Fig.8 Measured and extrapolated electric potential signals for z=22.5 m

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

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

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

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

5 結(jié)論

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

1)點(diǎn)電荷為標(biāo)量，而電偶極子為矢量(未知個(gè)數(shù)較多)，減少了反演時(shí)超定方程未知數(shù)的個(gè)數(shù)，場(chǎng)源的位置進(jìn)行了預(yù)先確定，進(jìn)一步減少了未知量的個(gè)數(shù)。

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

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

4)反演出的艦船電場(chǎng)源強(qiáng)度包含了船體表面電流大小的分布狀況，該信息可為艦船電場(chǎng)防護(hù)和隱身提供基礎(chǔ)數(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］林春生，龔沈光. 艦船物理場(chǎng)［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］陳聰，龔沈光，李定國(guó). 艦船靜態(tài)電場(chǎng)深度換算方法［J］. 哈爾濱工程大學(xué)學(xué)報(bào)，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］劉文寶，王向軍，稽斗. 基于電偶極子模型的艦船靜電場(chǎng)深度換算［J］. 空軍雷達(dá)學(xué)院學(xué)報(bào)，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］陳聰，李定國(guó)，龔沈光. 基于拉氏方程的艦船靜態(tài)電場(chǎng)深度換算［J］.電子學(xué)報(bào)，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.