基于邊界積分法的艦艇磁場(chǎng)模型參數(shù)確定

基于邊界積分法的艦艇磁場(chǎng)模型參數(shù)確定

0重點(diǎn)難以確定的磁場(chǎng)范圍難以延拓的歸位方法或原理消磁是提高潛水員耐磁性素質(zhì)的有效手段。評(píng)價(jià)消磁效果或者評(píng)估潛艇的磁性防護(hù)能力,一般都是針對(duì)潛艇水線以下的某種特定深度上某些特定點(diǎn)的磁感應(yīng)強(qiáng)度或潛艇上方的特定高度上的磁信號(hào)來(lái)進(jìn)行的。受實(shí)際測(cè)量條件(水深,風(fēng)浪和高度等)的限制,在實(shí)際測(cè)量中往往是很難直接獲得這類特定的磁場(chǎng)信息,所以需要從已知場(chǎng)點(diǎn)的磁場(chǎng)值來(lái)推算考核點(diǎn)上的磁場(chǎng)值。目前,用于艦艇磁場(chǎng)延拓的方法主要有有限元法、積分方程法和磁體模擬法。有限元和積分方程法都需要對(duì)艦船進(jìn)行剖分,其剖分的合理性和疏密程度嚴(yán)重影響著延拓結(jié)果。對(duì)于形體復(fù)雜且龐大的艦船而言,剖分是非常復(fù)雜困難的,并且得到的方程組也將非常龐大,運(yùn)算時(shí)間較長(zhǎng)。磁體模擬法常用的磁性模擬體主要有磁偶極子和磁橢球體,該方法易受人為經(jīng)驗(yàn)的影響,穩(wěn)定性和精度難以保證。人們?cè)?0世紀(jì)70年代將邊界積分方程用于對(duì)位場(chǎng)的研究,由于該方法只對(duì)求解區(qū)域邊界剖分,因而使問(wèn)題的維數(shù)降低,這就使代數(shù)方程組的維數(shù)大為減少,節(jié)約了計(jì)算成本。1磁體擴(kuò)展數(shù)學(xué)模型1.1各區(qū)域v內(nèi)磁場(chǎng)的各磁感應(yīng)強(qiáng)度計(jì)算如圖1所示,閉合面S將場(chǎng)域空間分成2部分V1和V,所有的磁源(包括鐵磁物質(zhì)和傳導(dǎo)電流)都被限制在V1內(nèi)。V可以看成是由2個(gè)閉合面限定的區(qū)域(從內(nèi)部限定的閉合面S和伸展到無(wú)限遠(yuǎn)的S∞),該區(qū)域內(nèi)沒(méi)有電流分布(J=0),標(biāo)量磁位φ滿足拉普拉斯方程??2φ=0.由于磁位φP在無(wú)窮遠(yuǎn)處定義為0,則直接邊界積分方程可以寫成:cPφP=∮S14πr?φ?ndS-14πr?φ?ndS?∮Sφ??nφ??n14πr14πrdS,(1)其中:cP對(duì)V內(nèi)的點(diǎn)為1,對(duì)光滑邊界S上的點(diǎn)為1/2,對(duì)區(qū)域V1內(nèi)的點(diǎn)為0;r是邊界面上任意點(diǎn)到場(chǎng)點(diǎn)的矢徑;S面上的外法線n從區(qū)域V指向V1.記φ*=?φ?nφ?=?φ?n,則從(1)式可以看出:V1內(nèi)的磁源對(duì)V中任意位置磁位φP的貢獻(xiàn)可以用邊界上的面積分來(lái)表示。若已知邊界面S上的等效源磁性參數(shù)(φ和φ*)后,可以由該式計(jì)算V中任何一點(diǎn)的標(biāo)量磁位φP.在得到該點(diǎn)的φP后,通過(guò)求取該點(diǎn)處的負(fù)梯度,并結(jié)合該點(diǎn)處的磁導(dǎo)率μ,即可得到該點(diǎn)處的磁感應(yīng)強(qiáng)度。BP=-μ??′φΡ=-μ4π′φP=?μ4π∮S??′1r?φ?n-φ??n1r?φ?n?φ??n1r1rdS,(2)其中,法線方向的偏微分是對(duì)邊界源點(diǎn)坐標(biāo)進(jìn)行的,而梯度??′是對(duì)場(chǎng)點(diǎn)坐標(biāo)進(jìn)行的。將邊界面S劃分為N個(gè)面積單元,則兩式可離散化為這就是邊界上的等效源磁性參數(shù)與無(wú)磁區(qū)域V中的磁感應(yīng)強(qiáng)度的正演關(guān)系。如果在區(qū)域V內(nèi)測(cè)量得到一組潛艇磁場(chǎng)感應(yīng)強(qiáng)度,就可以利用(4)式反演出邊界上等效的磁性參數(shù),再利用這組等效的磁性參數(shù)可以計(jì)算V區(qū)域內(nèi)任意位置的磁感應(yīng)強(qiáng)度。在反演的過(guò)程中,還需要考慮一些其他信息,如邊界面磁通連續(xù),φ和φ*的平衡,以及兩者在邊界上的連續(xù)性。1.2單元中心處的磁位平衡從測(cè)量值反演等效磁性參數(shù)時(shí),約束條件為:1)φ和φ*的平衡;2)磁通連續(xù)。如果將場(chǎng)點(diǎn)取在各個(gè)單元的中心,則cp=1/2(為了計(jì)算方便,一般取光滑的邊界),φ和φ*在邊界上的平衡,可以描述為φi為第i個(gè)單元中心處的磁位。由于所有磁源都被限制在閉合面S內(nèi),進(jìn)出S的磁力線條數(shù)相等。由高斯定理∮SBdS=0,經(jīng)推導(dǎo)并離散化后得到Ν∑j=1∫Sj∑j=1N∫Sj?φj?ndSj=0.(6)?φj?ndSj=0.(6)1.3型磁管磁模型的正演分析邊界形狀和邊界單元的劃分嚴(yán)重影響著正演和反演計(jì)算的速度和精度。出于單元?jiǎng)澐趾头e分計(jì)算簡(jiǎn)單考慮,一般都選取較規(guī)則的邊界(如長(zhǎng)方體表面,圖2).由于邊界S上沒(méi)有實(shí)際電荷和磁源,所以邊界上的φ和φ*連續(xù);同時(shí)邊界S存在拐點(diǎn)(如長(zhǎng)方體的棱),在這些拐點(diǎn)處φ*的方向未定。線性單元適合于φ和φ*在邊界上連續(xù)的場(chǎng)合,但是如果將邊界全部使用線性單元分割,解在拐點(diǎn)處的誤差就大。為了改善解的精度,本文將邊界元視為混合單元,即在每個(gè)單元內(nèi)φ是連續(xù)線性變化的,而φ*是定常的,此時(shí)(4)式和(5)式相應(yīng)的變?yōu)槭街?Q為單元節(jié)點(diǎn)數(shù);ωkiki和ωkjkj為加權(quán)系數(shù);φkiki和φkjkj為單元i和j的第k個(gè)節(jié)點(diǎn)處的磁位。ωkiki和Q由邊界面的具體形狀和分割方式?jīng)Q定,一般采用簡(jiǎn)單的四邊形元素(本文中為長(zhǎng)方形,見(jiàn)圖2),它由4個(gè)節(jié)點(diǎn)組成,此時(shí)Q=4,ωkiki可按下式計(jì)算:{ω1i=(1-ε1)(1-ε2)4,ω2i=(1+ε1)(1-ε2)4,ε3i=(1+ε1)(1+ε2)4,ε4i=(1-ε1)(1+ε2)4,-1≤ε1,ε2≤1.(9)?????????ω1i=(1?ε1)(1?ε2)4,ω2i=(1+ε1)(1?ε2)4,ε3i=(1+ε1)(1+ε2)4,ε4i=(1?ε1)(1+ε2)4,?1≤ε1,ε2≤1.(9)(7)式和(8)式中的積分項(xiàng)與場(chǎng)點(diǎn)Pi和邊界單元Sj有關(guān),積分方法與單元的具體形狀密切相關(guān),一般難以求得解析表達(dá)式,都采用數(shù)值積分方法。積分項(xiàng)的具體表達(dá)式可以參閱文獻(xiàn),這里不再贅述。對(duì)于(8)式中,應(yīng)當(dāng)注意的是:當(dāng)i=j時(shí),源點(diǎn)與節(jié)點(diǎn)在同一單元上,出現(xiàn)積分奇異,兩者的積分結(jié)果如下:(11)式只適用于圖2中的長(zhǎng)方形邊界單元,其中l(wèi)和d為長(zhǎng)方形的長(zhǎng)和寬。至此,作為延拓潛艇磁感應(yīng)強(qiáng)度所需的第1個(gè)步驟,我們完成了正演問(wèn)題的分析。延拓的具體過(guò)程可以描述為:在潛艇周圍空間測(cè)量得到m個(gè)不同位置的潛艇磁場(chǎng)三分量,則綜合(6)～(8)式形成方程組,寫成矩陣的形式為:Ax=b;求解方程組Ax=b反演得到邊界上一組等效的磁性參數(shù)后,再利用(7)式計(jì)算艦船在其他位置的磁場(chǎng)感應(yīng)強(qiáng)度。2模型求解和正則化參數(shù)的選取由上節(jié)模型得到的方程組都呈現(xiàn)嚴(yán)重的病態(tài)特性,磁感應(yīng)強(qiáng)度測(cè)量值的微小擾動(dòng)會(huì)導(dǎo)致解的很大變化。而在求解方程組時(shí),通常按線性最小二乘來(lái)求解,即求解法方程組ATAx=ATb.(12)若方程本身是病態(tài)的,則法方程會(huì)更加的病態(tài)。國(guó)內(nèi)外很多學(xué)者都對(duì)這類病態(tài)問(wèn)題進(jìn)行了研究,這些研究表明TSVD和Tikhonov正則化等方法可以較好的處理逆問(wèn)題中由系數(shù)矩陣病態(tài)特性而造成結(jié)果的不穩(wěn)定性。本文將采用Tikhonov正則化方法對(duì)方程組進(jìn)行求解。2.1Tikhonov正則化基本原理不妨假設(shè)實(shí)際磁場(chǎng)測(cè)量過(guò)程中測(cè)量誤差為ξ,理想情況下磁場(chǎng)測(cè)量值為c,則有b=c+ξ,Ax=c+ξ.對(duì)系數(shù)矩陣進(jìn)行奇異值分解A=USVT,則(12)式的解為x=k∑i=1∑i=1k(uΤic)σivi+(uTic)σivi+k∑i=1∑i=1k(uΤiξ)σivi,(13)其中ui和vi分別為U和V矩陣的第i列,σ1≥σ2≥…≥σk為A的奇異值,k=rank(A).分析(13)式可知,噪聲對(duì)真解的污染主要反映在右端第2項(xiàng),當(dāng)奇異值σi異常小時(shí),該項(xiàng)將被嚴(yán)重放大,甚至掩蓋真解(右端第1項(xiàng))。TSVD方法把容易造成不穩(wěn)定的較小奇異值直接截去,反映在物理實(shí)際中就是將該單元內(nèi)的磁化參數(shù)置為0,加劇了磁化參數(shù)的不連續(xù)性。而Tikhonov方法的處理思路:在奇異值上疊加一個(gè)極小的正實(shí)數(shù)α(正則化參數(shù)),抑制異常小的奇異值σi對(duì)噪聲的放大作用,達(dá)到正則濾波的目的。其相應(yīng)的解為xα=k∑i=1σi(uΤic)σ2i+αvi+k∑i=1σi(uΤiξ)σ2i+αvi.(14)從(14)式可以看出,若正則化參數(shù)α選取較大(過(guò)正則化),則xα與真解差別較大,模型的精確性降低;若正則化參數(shù)α選取較小(欠正則化),則無(wú)法抑制對(duì)噪聲的放大作用,模型的穩(wěn)定性下降。所以恰當(dāng)?shù)恼齽t化參數(shù)需要在模型的精確性和穩(wěn)定性之間做出一個(gè)合適的折中。2.2正則化參數(shù)的選取正則化參數(shù)的選取方法分為兩類:先驗(yàn)選取法和后驗(yàn)選取法。先驗(yàn)選取方法有Morozov偏差原理,該類方法需要知道與測(cè)量數(shù)據(jù)相關(guān)的先驗(yàn)信息(如測(cè)量數(shù)據(jù)的誤差水平),然而這些先驗(yàn)信息在潛艇磁場(chǎng)反演建模中非常少。后驗(yàn)選取方法主要有L曲線法,交叉校驗(yàn)(CV)和廣義交叉校驗(yàn)(GCV).這3種方法直接從測(cè)量數(shù)據(jù)出發(fā),通過(guò)預(yù)解方程組來(lái)確定正則化參數(shù),而不需要先驗(yàn)信息,適合本文中逆問(wèn)題的求解。1)L曲線法如果以‖xa‖為橫坐標(biāo),‖Axa-b‖為縱坐標(biāo)繪成的曲線為L(zhǎng)形狀,稱為L(zhǎng)曲線。L曲線反映了逆問(wèn)題模型的估計(jì)參數(shù)的精確性和穩(wěn)定性隨著正則化參數(shù)α的變化規(guī)律。一般選取L曲線上曲率最大的點(diǎn)(都為L(zhǎng)曲線的拐點(diǎn))所對(duì)應(yīng)的α作為正則化參數(shù),但是磁場(chǎng)延拓模型實(shí)際的L曲線在拐點(diǎn)附近集聚了大量的α,最優(yōu)的α值難以確定。所以本文借助極小化如下輔助函數(shù)來(lái)選取α.ρ(α)=lg(‖Axα-b‖)+lg(‖xα‖),(15)磁場(chǎng)延拓模型的L曲線和輔助函數(shù)ρ(α)曲線見(jiàn)圖3和圖4.2)CV和GCV定義如下的CV和GCV函數(shù):CVα=1nn∑i=1bi-?bi1-Ηii(α)2,(16)GCVα=1nn∑i=1bi-?bi1-tr(Η(α))/n2,(17)式中:Hα=A(ATA+αI)-1AT;■=Hαb;tr(H(α))為H(α)的跡;Hii(α)是矩陣H(α)的對(duì)角元。使函數(shù)CV或GCV取得最小值的α就是由交叉校驗(yàn)或廣義交叉校驗(yàn)所求得的正則化參數(shù)。潛艇磁場(chǎng)延拓模型的CV函數(shù)和GCV函數(shù)曲線見(jiàn)圖5和圖6.從圖4～6可以看出,CV函數(shù)確定的α最小,其次是GCV函數(shù),最大的是L曲線。大量分析發(fā)現(xiàn):CV函數(shù)確定的模型精度較高,但魯棒性較差;L曲線確定的模型魯棒性較強(qiáng),但精度較差;綜合比較之下,GCV函數(shù)選取的正則化參數(shù)比較恰當(dāng)。本文在船模實(shí)驗(yàn)中,采用GCV函數(shù)確定正則化參數(shù),取得了較為滿意的效果。3船模右病z1-12模以一個(gè)長(zhǎng)308cm、寬34cm的鐵質(zhì)結(jié)構(gòu)船模為實(shí)驗(yàn)對(duì)象,采用三分量磁通門傳感器陣列進(jìn)行測(cè)量。在4m×0.4m的測(cè)量范圍內(nèi),以等間距分別測(cè)量得到21×5點(diǎn)陣列上的磁感應(yīng)強(qiáng)度,如圖7所示(船模中心在測(cè)量平面上的投影與平面中心重合,測(cè)量平面高度分別為z1=57.8cm和z2=72cm).閉合邊界面取包含船模的長(zhǎng)方體表面,測(cè)量平面較其在水平面的投影稍大,整個(gè)邊界共劃分為280個(gè)單元。分別利用z1(z2)平面上的磁感應(yīng)強(qiáng)度三分量的測(cè)量數(shù)據(jù),建立磁場(chǎng)模型,并由模型計(jì)算得到了z2(z1)平面磁感應(yīng)強(qiáng)度的三分量值。定義相對(duì)殘差:‖Bm-Bc‖2/‖Bm,max‖2,推算誤差:Bci/Bm,max×100%,其中Bm和Bc分別為磁感應(yīng)強(qiáng)度的測(cè)量值向量和推算值向量,Bm,max為磁感應(yīng)強(qiáng)度測(cè)量值的最大值,Bci為第i點(diǎn)磁感應(yīng)強(qiáng)度的推算值。從表1可以看出:X分量的相對(duì)殘差最小,Z分量其次,Y分量最大,這是因?yàn)樵赮方向的磁場(chǎng)信息比X方向和Z方向少造成的。圖8和圖9分別給出了近場(chǎng)到遠(yuǎn)場(chǎng)及遠(yuǎn)場(chǎng)到近場(chǎng)的推算值與測(cè)量值的對(duì)比曲線(由于篇幅的原因,只畫出潛艇龍骨下的磁場(chǎng)比較圖,其他位置有同樣的性質(zhì))。從圖中可以看出,由近場(chǎng)推算遠(yuǎn)場(chǎng)時(shí),3個(gè)分量的推算誤差均可控制在6%以內(nèi),而從遠(yuǎn)場(chǎng)推算近場(chǎng)模型誤差稍大,但均可控制在9%以內(nèi)。圖10和圖11分別給出了經(jīng)正則化處理和未經(jīng)正則化處理2種情況下,從近場(chǎng)到遠(yuǎn)場(chǎng)和從遠(yuǎn)場(chǎng)到近場(chǎng),船模右舷z分量的推算值與測(cè)量值的對(duì)比曲線(其他位置和其他分量具有同樣的性質(zhì))。從圖中可以看出,采用GCV函數(shù)正則