衛(wèi)星載荷磁場優(yōu)化計(jì)算方法研究

李 恪，姚崇斌，徐紅新，謝寶蓉

(上海航天電子技術(shù)研究所，上海 201109)

0 引言

地磁測量是衛(wèi)星弱磁探測技術(shù)的重要應(yīng)用之一，對整個(gè)地球的磁場環(huán)境描述有著重要的意義[1-4]。因衛(wèi)星與地球相距非常遠(yuǎn)，故衛(wèi)星所處空間的地磁感應(yīng)強(qiáng)度很小。但衛(wèi)星在軌工作時(shí)，整星內(nèi)部會形成一個(gè)基本穩(wěn)定的磁場。衛(wèi)星內(nèi)部磁場與地磁場相互作用形成磁場干擾力矩，嚴(yán)重影響姿態(tài)控制及地磁測量的精度。為減少星上磁場噪聲的影響，需要在地面上測量載荷的磁場，并做消磁處理。

現(xiàn)有的磁場測量方法大多采用數(shù)值計(jì)算方法，將被測物體劃分為許多小的磁場單元，對各個(gè)單元在不同方向上的磁場分量進(jìn)行積分，從而得到被測物體總的磁場，如有限元法、積分方程法等。有限元法要求在整個(gè)求解區(qū)域內(nèi)進(jìn)行離散剖分，計(jì)算數(shù)據(jù)的前處理較繁雜，且計(jì)算開域問題時(shí)需要截?cái)噙吔纾视邢拊ǖ氖褂梅秶芟蕖７e分方程法中系數(shù)矩陣元素的計(jì)算為三重體積分，計(jì)算較復(fù)雜，且當(dāng)源點(diǎn)和場點(diǎn)重合時(shí)，系數(shù)矩陣的計(jì)算帶有奇異積分。該方法的優(yōu)點(diǎn)是剖分簡單，待求解未知量少，但非對稱滿陣導(dǎo)致占用資源大，計(jì)算時(shí)間長。

地磁場是一個(gè)弱磁場，對于處于不同地磁場中的載荷，可認(rèn)為其磁化率基本不變。故針對數(shù)值計(jì)算方法的不足，根據(jù)實(shí)際應(yīng)用的需求，采用三維靜磁場逆計(jì)算方法的思路[5-7]，對傳統(tǒng)的數(shù)值計(jì)算方法加以優(yōu)化。

本算法利用某處地磁場下衛(wèi)星載荷的磁場測量數(shù)據(jù)，逆計(jì)算衛(wèi)星載荷的等效磁化率，根據(jù)等效磁化率計(jì)算衛(wèi)星載荷產(chǎn)生的感應(yīng)磁場，根據(jù)均勻磁化衛(wèi)星載荷的特點(diǎn)，將積分方程法中系數(shù)矩陣元素的體積分形式化為面積分形式，并消除積分計(jì)算的奇異性，建立基于單元表面積分的感應(yīng)磁場計(jì)算模型。最后，以實(shí)際載荷為實(shí)驗(yàn)對象，驗(yàn)證算法的有效性。相比傳統(tǒng)算法，在保證精度的前提下，極大地提高了工作效率。

1 基于表面積分的均勻磁化體磁場計(jì)算模型

基于體積分的積分方程法磁場計(jì)算公式為[2]

Hm(r)=

(1)

式中：Hm為鐵磁物體產(chǎn)生的磁場；M為鐵磁物體內(nèi)部磁化強(qiáng)度；r為場點(diǎn)的坐標(biāo)；r′為源點(diǎn)的坐標(biāo)；V為鐵磁物體的體積。由式(1)計(jì)算鐵磁物體的磁場時(shí)需要計(jì)算三重積分，計(jì)算較復(fù)雜，且當(dāng)場點(diǎn)和源點(diǎn)重合時(shí)，積分計(jì)算具有奇異性。

在不存在自由電流的區(qū)域中，從磁荷觀點(diǎn)出發(fā)，均勻磁化衛(wèi)星載荷產(chǎn)生的標(biāo)量磁位[8]可表示為

(2)

(3)

式中：dS為單位面元的投影面積；n為曲面S在r處的法向矢量。故均勻磁化衛(wèi)星載荷在空間產(chǎn)生的磁場強(qiáng)度可表示為

(4)

一般情況下，不規(guī)則衛(wèi)星載荷在外磁場中未被均勻磁化。若將衛(wèi)星載荷剖分成N個(gè)小單元，當(dāng)剖分單元體積足夠小時(shí)，則可認(rèn)為每個(gè)小單元是等效均勻磁化體。由磁場疊加原理可知，衛(wèi)星載荷在空間產(chǎn)生的磁場強(qiáng)度為所有等效均勻磁化體的疊加，即式(4)可表示為

(5)

考慮到每個(gè)剖分單元衛(wèi)星載荷的線性度和均勻性，設(shè)每個(gè)剖分單元滿足Mi=χiHi， 式(5)可進(jìn)一步表示為

(6)

式中:χ=[χ1,χ2,…,χN]為待優(yōu)化的參數(shù)組。一般情況下，該問題很難直接求解，但通過迭代過程，即多次計(jì)算不同的χ所對應(yīng)的目標(biāo)函數(shù)值，即可逼近問題的真實(shí)解。當(dāng)已知每個(gè)剖分單元的磁化率和磁場強(qiáng)度時(shí)，根據(jù)式(6)就可計(jì)算出衛(wèi)星載荷在空間任一點(diǎn)產(chǎn)生的感應(yīng)磁場。

假設(shè)衛(wèi)星載荷處于外磁場為He的空間，則空間任一點(diǎn)的磁場H=He+Hm。當(dāng)把場點(diǎn)配置在第j個(gè)剖分單元的中心時(shí)，可得到第j個(gè)剖分單元中心磁場強(qiáng)度的表達(dá)式為

(7)

寫成矩陣形式為AH=He，其中系數(shù)矩陣A由場點(diǎn)坐標(biāo)、剖分單元中心坐標(biāo)及單元磁化率決定，H為每個(gè)剖分單元中心的磁場強(qiáng)度。

已知衛(wèi)星載荷的磁性參數(shù)，根據(jù)式(7)建立每個(gè)單元中心處磁場強(qiáng)度的線性方程組，求解該方程組即可得到每個(gè)剖分單元的磁場強(qiáng)度，再由式(6)就可計(jì)算衛(wèi)星載荷在空間任一點(diǎn)產(chǎn)生的感應(yīng)磁場。由此可看出，系數(shù)矩陣元素的計(jì)算為面積分的形式，較傳統(tǒng)的體積分形式簡單；由于積分是在剖分單元表面進(jìn)行，而場點(diǎn)選取在剖分單元的中心，故當(dāng)源點(diǎn)和場點(diǎn)重合時(shí)，積分計(jì)算無奇異性。

2 基于微粒群算法的等效磁化率逆計(jì)算方法

2.1 微粒群算法的基本原理

微粒群算法中，所有粒子都有一個(gè)由被優(yōu)化函數(shù)決定的適應(yīng)值，每個(gè)粒子還有一個(gè)速度用于決定他們搜索的方向和距離，粒子們追隨當(dāng)前的最優(yōu)粒子在解空間中搜索。在求解過程中，微粒群初始化為一群隨機(jī)粒子(隨機(jī)解)，并通過迭代找到最優(yōu)解。在每一次迭代中，粒子通過跟蹤兩個(gè)極值來更新自身，第一個(gè)極值是粒子本身所找到的最優(yōu)解，即個(gè)體極值；另一個(gè)極值是整個(gè)種群目前找到的最優(yōu)解，即全局極值。找到這兩個(gè)最優(yōu)值后,粒子根據(jù)如下的公式來更新自己的速度和新的位置[9]

(8)

式中：vi為粒子的速度；xi為當(dāng)前粒子的位置；pbest為個(gè)體極值；gbest為全局極值；xrand為介于0與1之間的隨機(jī)數(shù)；c1、c2為學(xué)習(xí)因子；w為慣性權(quán)重；t為迭代的次數(shù)。每一維粒子的速度都會受限于最大速度vmax，如果粒子某一維更新后的速度超過vmax，那么這一維的速度就被限定為vmax。

2.2 計(jì)算方法與流程

由于外加磁場的影響，衛(wèi)星載荷的磁性狀態(tài)會根據(jù)外磁場的大小，沿著其磁化曲線工作于相應(yīng)的磁化點(diǎn)。電磁場逆問題是由果求因的問題[10]，在靜磁場計(jì)算中即為已知磁場測量數(shù)據(jù)來反演衛(wèi)星載荷的磁性參數(shù)。建立的逆問題數(shù)學(xué)模型可表示為

(9)

由式(9)可知，模型中有N個(gè)未知數(shù)，故將粒子的解空間設(shè)置成N維空間，每個(gè)粒子的坐標(biāo)為[χ1,χ2,…,χN]，粒子的適應(yīng)度函數(shù)取為式(9)中的J，每個(gè)粒子具有的速度v=[v1,v2,…,vN]，第t+1代粒子在第t代粒子的基礎(chǔ)上按式(8)更新自己的位置和速度，其優(yōu)化過程如下：

1) 設(shè)置相關(guān)初始參數(shù)，如粒子群規(guī)模m、慣性權(quán)重w、學(xué)習(xí)因子c1和c2、粒子最大飛行速度vmax、迭代結(jié)束條件；

2) 初始化粒子的位置及速度；

3) 評價(jià)各粒子的適應(yīng)度函數(shù)；

4) 對每個(gè)微粒，將其適應(yīng)度值與歷史最好位置pbest相比較，若當(dāng)前適應(yīng)度值更優(yōu)，則用當(dāng)前適應(yīng)度值更新pbest；

5) 將每個(gè)微粒的適應(yīng)度值與群體經(jīng)歷過的歷史最佳位置gbest相比較，若當(dāng)前群體中最好的適應(yīng)度值較好,則將其設(shè)置為新的gbest；

6) 根據(jù)式(8)更新各微粒的位置和速度；

7) 計(jì)算優(yōu)化目標(biāo)函數(shù)，若達(dá)到精度要求或迭代次數(shù)，則返回當(dāng)前最佳微粒的結(jié)果，算法結(jié)束，否則返回步驟3)。

3 實(shí)驗(yàn)方案與結(jié)果分析

3.1 實(shí)驗(yàn)方案

圖1 實(shí)驗(yàn)方案及現(xiàn)場示意圖Fig.1 Schematic diagrams of experimental scheme and experimental site

為考察不同場點(diǎn)、不同外磁場及不同磁化方向下計(jì)算的感應(yīng)磁場的準(zhǔn)確性，測量時(shí)，被測物體先沿南北方向滑動測量，再轉(zhuǎn)動測量架，使被測物體沿東西方向滑動測量。地磁平面的選取標(biāo)準(zhǔn)原則上是離被測物體越近越好，但是要考慮到其他設(shè)備的干擾，因此需要綜合考慮取測量偏差較小的距離。用上述相同的方法測量縱向載荷上端高度為h2、載荷下端高度為h3，以及橫向載荷下端高度為h4時(shí)的平面磁場數(shù)據(jù)，以用于與計(jì)算值比較。

3.2 實(shí)例計(jì)算過程說明

應(yīng)用本文算法計(jì)算衛(wèi)星載荷感應(yīng)磁場過程中，兼顧計(jì)算時(shí)間和計(jì)算精度，將載荷剖分為20×5個(gè)長方體單元，并假定每個(gè)剖分單元為等效均勻磁化體。應(yīng)用微粒群優(yōu)化算法反演計(jì)算磁化率的過程中，考慮到載荷材料的均勻性，認(rèn)為每個(gè)剖分單元在地磁場的作用下的磁化率相同，故式(8)可簡化為只含一個(gè)未知數(shù)的優(yōu)化問題，相應(yīng)的粒子的解空間為1維空間，即待優(yōu)化的參變量χ=[χ,χ,…,χ]。

用微粒群優(yōu)化算法進(jìn)行優(yōu)化求解時(shí)，參數(shù)設(shè)定為：微粒粒子數(shù)m=40，慣性權(quán)重w=0.9，學(xué)習(xí)因子c1=c2=2，最大速度vmax=10，迭代次數(shù)50。為考察微粒群算法反演磁化率的有效性， 逆計(jì)算15次。圖2為15次逆計(jì)算所得的等效磁化率和最終最優(yōu)粒子的適應(yīng)度曲線和一次逆計(jì)算時(shí)粒子適應(yīng)度的收斂曲線。由圖2可知，用微粒群優(yōu)化算法反演磁化率是有效的，且收斂速度較快。計(jì)算衛(wèi)星載荷感應(yīng)磁場時(shí)所用磁化率取15次逆計(jì)算得到的等效磁化率的平均值為140.8。圖3為衛(wèi)星載荷受地磁場水平分量縱向磁化和橫向磁化時(shí)，產(chǎn)生的感應(yīng)磁場三分量測量值和計(jì)算值對比曲線。由圖3可看出，測量值和計(jì)算值吻合較好，其最大相對誤差約為10.79%。

圖2 磁化率逆計(jì)算參數(shù)圖和收斂曲線Fig.2 Parameters and convergence curve of magnetic susceptibility inverse calculation

4 結(jié)束語

由均勻磁化體的磁化特點(diǎn)出發(fā)，建立了基于單元表面積分的衛(wèi)星載荷感應(yīng)磁場求解方法。與傳統(tǒng)基于體積分的磁場求解方法相比，該方法簡化了系數(shù)矩陣運(yùn)算，消除了源點(diǎn)和場點(diǎn)重合時(shí)的積分奇異性。以衛(wèi)星載荷的磁化率為未知量，建立了求解磁化率的數(shù)學(xué)模型，并運(yùn)用隨機(jī)類微粒群優(yōu)化算法對其求解。最后，以衛(wèi)星載荷為計(jì)算實(shí)例，驗(yàn)證了本文算法的可行性，計(jì)算值與測量值具有較好的一致性。結(jié)果表明，該方法完全可用于衛(wèi)星載荷感應(yīng)磁場的計(jì)算。

圖3 不同磁化方向、不同場點(diǎn)的感應(yīng)磁場Fig.3 Induction magnetic field at different magnetization directions and field points

[1] 吳衛(wèi)權(quán). 地磁場中衛(wèi)星不側(cè)置態(tài)下磁矩測試方法[J]. 上海航天, 2014, 31(1): 50-55.

[2] 徐文耀. 地磁場的三維巡測和綜合建模[J]. 地球物理學(xué)進(jìn)展, 2007, 22(4): 1035-1039.

[3] 傅綏燕， 宗秋剛, 濮祖蔭. 地磁活動與太陽活動對環(huán)電流離子成分與位置的影響[J]. 地球物理學(xué)報(bào), 2003, 46(6): 725-730.

[4] 林瑞淋, 張效信, 劉四清. 高緯磁層頂位形統(tǒng)計(jì)分析[J]. 地球物理學(xué)報(bào), 2003, 53(1): 1-9.

[5] PECHENKOV A N . A mathematical algorithm for solving the inverse problem of magnetostatic flaw detection[J]. Russian Journal of Nondestructive Testing, 2005, 41(11): 714-718.

[6] PECHENKOV A N. Numerical simulation of an inverse problem of magnetostatic flaw detection[J]. Russian Journal of Nondestructive Testing, 2005, 41(11): 719-723.

[7] 周克定. 工程電磁場專論[M]. 武漢: 華中工學(xué)院出版社, 1986: 28-29.

[8] 樊明武, 顏威利. 電磁場積分方程法[M]. 北京: 機(jī)械工業(yè)出版社, 1988: 65-68.

[9] 連麗婷, 楊明明, 孫開江. 基于改進(jìn)微粒群算法的磁靶電流參數(shù)優(yōu)化調(diào)整[J]. 艦船科學(xué)技術(shù), 2015, 37(3): 73-76.

[10] 黃卡瑪, 趙翔. 電磁場中的逆問題及應(yīng)用[M]. 北京: 科學(xué)出版社, 2005: 13-19.