孟慶奎,周堅(jiān)鑫,舒晴,高維,徐光晶,王晨陽(yáng)

(1.自然資源部 航空地球物理與遙感地質(zhì)重點(diǎn)實(shí)驗(yàn)室，北京 100083; 2.中國(guó)自然資源航空物探遙感中心，北京 100083)

0 引言

實(shí)驗(yàn)室條件下，基于超導(dǎo)量子干涉器(superconducting quantum interference device，SQUID)的全張量磁力梯度測(cè)量系統(tǒng)的靈敏度可以達(dá)到幾個(gè)fT的量級(jí)，能實(shí)現(xiàn)對(duì)微弱磁性目標(biāo)體的定位與識(shí)別[1-2]，為目標(biāo)探測(cè)提供有效手段。21世紀(jì)以來(lái)，航空全張量磁力梯度測(cè)量系統(tǒng)研發(fā)及配套數(shù)據(jù)處理解釋理論是國(guó)內(nèi)外相關(guān)領(lǐng)域?qū)W者研究的熱點(diǎn)，德國(guó)、美國(guó)、澳大利亞等發(fā)達(dá)國(guó)家對(duì)全張量磁力梯度測(cè)量技術(shù)開展了大量探索研發(fā)工作。目前，國(guó)際上僅德國(guó)光學(xué)分子研究中心擁有實(shí)用化水平的航空全張量磁力梯度測(cè)量系統(tǒng)[3]，且該技術(shù)屬于西方國(guó)家限制出口的尖端技術(shù)，公開的文獻(xiàn)資料也非常有限。面對(duì)技術(shù)封鎖，我國(guó)自力更生，在“國(guó)家重大科研裝備研制項(xiàng)目”和“國(guó)家863計(jì)劃”支持下，中國(guó)科學(xué)院上海微系統(tǒng)與信息技術(shù)研究所及吉林大學(xué)研發(fā)團(tuán)隊(duì)分別研制出基于低溫和高溫SQUID的航空全張量磁力梯度測(cè)量系統(tǒng)，研發(fā)團(tuán)隊(duì)攻克了串?dāng)_、平臺(tái)震動(dòng)監(jiān)測(cè)、磁張量補(bǔ)償?shù)榷囗?xiàng)技術(shù)難關(guān)[4-9]，取得了長(zhǎng)足進(jìn)步，但離實(shí)用化還存在一定的差距。我國(guó)學(xué)者在全張量磁力梯度數(shù)據(jù)處理[10-11]、正反演[12-15]、綜合應(yīng)用[16-19]等方面也取得了系列成果。

航空全張量磁力梯度測(cè)量數(shù)據(jù)中包含非常復(fù)雜的運(yùn)動(dòng)噪聲，在頻譜圖中從低頻至高頻均有分布，并以白噪聲為主，如何有效壓制運(yùn)動(dòng)噪聲是一個(gè)較大的挑戰(zhàn)[3]。Stolz R總結(jié)了運(yùn)動(dòng)噪聲來(lái)源，包括儀器本身幾何位置變化、外部環(huán)境噪聲、飛機(jī)發(fā)動(dòng)機(jī)噪聲以及儀器本證噪聲[3]，航磁軟補(bǔ)償技術(shù)可以去除運(yùn)動(dòng)載體在地磁場(chǎng)中產(chǎn)生的恒定磁場(chǎng)、感應(yīng)磁場(chǎng)和渦流磁場(chǎng)等干擾場(chǎng)[20]，對(duì)于白噪聲卻是無(wú)能為力。傳統(tǒng)的數(shù)字濾波只能濾除指定頻段的噪聲，對(duì)于混疊在全張量磁力梯度有用信號(hào)中的噪聲不能有效分離?？柭鼮V波技術(shù)在航空重力測(cè)量中得到了應(yīng)用，實(shí)踐表明，該方法可以有效提取淹沒在載體高頻振動(dòng)噪聲中的有用信號(hào)[21-23]，將其推廣應(yīng)用到航空全張量磁力梯度數(shù)據(jù)處理中十分必要。張興東將卡爾曼濾波應(yīng)用到全張量磁力梯度數(shù)據(jù)處理中[9]，但是沒有明確狀態(tài)方程和觀測(cè)方程的搭建方式，而不合理的狀態(tài)方程會(huì)直接影響到卡爾曼濾波效果。筆者借鑒時(shí)間序列分析思想，基于最小差分方差建立狀態(tài)方程，并將卡爾曼濾波與數(shù)字濾波相結(jié)合，通過(guò)模型計(jì)算驗(yàn)證了方法的有效性，全區(qū)噪聲衰減因子優(yōu)于0.92，即能夠去除92%以上噪聲成分，全區(qū)均方誤差優(yōu)于10 pT/m，滿足航空全張量磁力梯度測(cè)量精度要求。卡爾曼濾波是一種快速、高效的遞推算法，利用實(shí)時(shí)測(cè)量信息對(duì)估計(jì)值進(jìn)行最優(yōu)估計(jì)，可直接應(yīng)用于航空全張量磁力梯度測(cè)量實(shí)時(shí)數(shù)據(jù)處理。

1 卡爾曼濾波方程建立

1.1 卡爾曼濾波基本理論

濾波估計(jì)經(jīng)歷了最小二乘、溫納濾波和卡爾曼濾波的發(fā)展而不斷完善。1960年，卡爾曼提出離散系統(tǒng)卡爾曼濾波，次年與布西合作，將這種濾波方法推廣到連續(xù)時(shí)間系統(tǒng)中，從而形成了卡爾曼濾波設(shè)計(jì)理論[24]?？柭鼮V波是一種時(shí)間域?yàn)V波方法，采用狀態(tài)空間方法描述系統(tǒng)，算法采用遞推形式，數(shù)據(jù)存儲(chǔ)量小，不僅可以處理平穩(wěn)隨機(jī)過(guò)程，也可以處理多維和非平穩(wěn)隨機(jī)過(guò)程。

卡爾曼濾波采用如下空間模型描述動(dòng)態(tài)系統(tǒng)：

X(k)=ΦX(k-1)+ΓW(k-1)，

(1)

Y(k)=HX(k)+V(k)，

(2)

式中：k為離散時(shí)間，系統(tǒng)在時(shí)刻k的狀態(tài)為X(k)；Y(k)為對(duì)應(yīng)狀態(tài)的觀測(cè)信號(hào)；W(k)為輸入的白噪聲；V(k)為觀測(cè)噪聲；Φ為狀態(tài)轉(zhuǎn)移矩陣；Γ為噪聲驅(qū)動(dòng)矩陣；H為觀測(cè)矩陣。式(1)為狀態(tài)方程，式(2)為觀測(cè)方程。

(3)

P(k,k-1)=ΦP(k-1)ΦT+ΓQΓT，

(4)

K(k)=P(k,k-1)HT[HP(k,k-1)HT+R]-1，

(5)

(6)

P(k)=[I-K(k)H]P(k,k-1)，

(7)

1.2 全張量磁力梯度數(shù)據(jù)卡爾曼濾波方程建立

全張量磁力梯度儀在每個(gè)采樣點(diǎn)測(cè)量6個(gè)梯度值，定義卡爾曼濾波器的狀態(tài)序列為：

X=(Bxx,Bxy,Bxz,Byy,Byz,Bzz)T，

(8)

各個(gè)測(cè)量時(shí)刻的測(cè)量序列為：

Z=(Bxx,Bxy,Bxz,Byy,Byz,Bzz)T，

(9)

因?yàn)闇y(cè)量序列與狀態(tài)序列變量相同，所以可以容易的定義觀測(cè)矩陣為：

H=diag(1,1,1,1,1,1)，

(10)

式中：diag表示對(duì)角矩陣。狀態(tài)轉(zhuǎn)移矩陣Φ的定義是一個(gè)難點(diǎn)，因全張量磁力梯度儀測(cè)量的數(shù)據(jù)是隨時(shí)間變化的變量，故可借鑒時(shí)間序列分析思想，計(jì)算包含白噪聲的梯度測(cè)量值的各階次差分方差，方差最小的階次差分的殘差可視為白噪聲，根據(jù)該階次差分方程可構(gòu)建狀態(tài)方程。經(jīng)大量模型計(jì)算，一階差分方差最小，故可定義狀態(tài)轉(zhuǎn)移矩陣為：

Φ=diag(1,1,1,1,1,1)。

(11)

本文中卡爾曼濾波進(jìn)行迭代的本質(zhì)是利用兩個(gè)正態(tài)分布的融合仍是正態(tài)分布的特性而進(jìn)行的，所指的兩個(gè)正態(tài)分布分別為系統(tǒng)輸入噪聲和測(cè)量噪聲的正態(tài)分布，系統(tǒng)輸入噪聲是指我們建立的系統(tǒng)狀態(tài)方程本身存在的系統(tǒng)誤差，測(cè)量噪聲是指測(cè)量?jī)x器精度、外界溫度變化等因素產(chǎn)生的測(cè)量誤差。融合后的正態(tài)分布所選擇的權(quán)值是指卡爾曼增益矩陣，其計(jì)算公式可參見式(5)。

2 模型算例分析

2.1 球體模型正演

自然界中存在諸多尺度有限的磁性地質(zhì)體，當(dāng)其中心埋深遠(yuǎn)遠(yuǎn)大于其直徑時(shí)，地質(zhì)體所引起的磁異常特征與球體磁場(chǎng)特征基本吻合，所以研究球體的磁場(chǎng)特征具有一定的實(shí)際意義。根據(jù)重磁位泊松公式，由引力位計(jì)算磁場(chǎng)三分量，進(jìn)而對(duì)磁場(chǎng)三分量在3個(gè)不同方向(X、Y、Z)求導(dǎo)便可得到球體的全張量梯度。

模型參數(shù)：測(cè)區(qū)范圍2 000 m×2 000 m，測(cè)量間距5 m，球心位置(1 000 m，1 000 m，550 m)，球體半徑500 m，總磁化強(qiáng)度傾角65°，偏角11°，總磁化強(qiáng)度0.5 A/m。正演模擬結(jié)果見圖1所示。

2.2 卡爾曼濾波效果

在上述正演結(jié)果中加入高斯白噪聲V(k)，視為觀測(cè)噪聲，V(k)均值為零，標(biāo)準(zhǔn)差為0.032 nT/m。加入噪聲的球體模型全張量磁力梯度視為測(cè)量值，如圖2所示，可見測(cè)量結(jié)果全區(qū)域都受到高斯白噪聲影響，且較嚴(yán)重。遵循以上卡爾曼濾波狀態(tài)方程和測(cè)量方程，建立遞推公式，對(duì)測(cè)量結(jié)果進(jìn)行卡爾曼濾波。在具體操作中，為模擬航空全張量磁力梯度實(shí)時(shí)測(cè)量與實(shí)時(shí)處理，從第1條測(cè)線(即剖面Y=0)的第1個(gè)測(cè)點(diǎn)開始，隨著數(shù)據(jù)的采集，時(shí)間的推移，逐點(diǎn)進(jìn)行卡爾曼濾波處理，直到全部數(shù)據(jù)采集完畢，同時(shí)完成濾波處理，體現(xiàn)了濾波技術(shù)的實(shí)時(shí)性。在不改變異常形態(tài)的前提下，為得到更加光滑的濾波效果，卡爾曼濾波后進(jìn)行小尺度中值濾波，最終的濾波結(jié)果見圖3～6所示。

為了定量評(píng)價(jià)本文所提出方法的濾波效果，進(jìn)一步評(píng)價(jià)卡爾曼濾波的降噪性能，引入Pajot等[25]定義的噪聲衰減因子β來(lái)估計(jì)噪聲衰減量：

β=[var(Bzs-B)-var(Bkf-B)]/var(Bzs-B)，

(12)

式中:Bzs表示含噪聲測(cè)量值，B表示理論值，Bkf表示卡爾曼濾波后的值。計(jì)算6個(gè)梯度張量的全區(qū)噪聲衰減因子，計(jì)算結(jié)果見圖6a所示，可見β值優(yōu)于0.92，即能夠去除92%以上噪聲成分。同時(shí)，計(jì)算6個(gè)梯度張量的全區(qū)均方誤差，結(jié)果見圖6b所示，均優(yōu)于10 pT/m。另外，計(jì)算了典型剖面(Y=500 m)上6個(gè)梯度張量的噪聲衰減因子(0.921 2～0.966 8)和均方誤差(6.2～8.5 pT/m)，具體結(jié)果見表1所示。以上定量評(píng)價(jià)結(jié)果表明，本文所提方法的濾波效果能夠滿足當(dāng)前航空全張量磁力梯度測(cè)量精度要求。

圖1 球體模型全張量磁力梯度Fig.1 Full tensor magnetic gradient of sphere model

圖2 加噪球體模型全張量磁力梯度Fig.2 Full tensor magnetic gradient of noisy sphere model

圖3 濾波后全張量磁力梯度立體圖Fig.3 Filtered full tensor magnetic gradient of sphere model

圖4 Bxx、Byy和Bzz濾波前后對(duì)比(剖面Y=500 m)Fig.4 Comparison of Bxx ，Byyand Bzz with before and after filtering(profile Y=500 m)

表1 典型剖面(Y=500 m)卡爾曼濾波效果定量統(tǒng)計(jì)

圖5 Bxy和Bxz和Byz濾波前后對(duì)比(剖面Y=500 m)Fig.5 Comparison of Bxy ，Bxz and Byzwith before and after filtering (profile Y=500 m)

圖6 全區(qū)噪聲衰減因子及全區(qū)均方誤差Fig.6 Full area noise-reduction factor and mean square error

3 結(jié)論

1) 本文通過(guò)理論模型驗(yàn)證了卡爾曼濾波方法在航空全張量梯度測(cè)量數(shù)據(jù)中的有效性，結(jié)果顯示全區(qū)噪聲衰減因子優(yōu)于0.92，全區(qū)均方誤差優(yōu)于10 pT/m。

2) 由于因卡爾曼濾波對(duì)數(shù)據(jù)的使用率低及噪聲隨機(jī)特征影響，單純應(yīng)用卡爾曼濾波不能得到平滑的預(yù)測(cè)結(jié)果，需要串聯(lián)其他適合的平滑濾波方法，達(dá)到更好的濾波效果。

3) 卡爾曼濾波可以利用實(shí)時(shí)測(cè)量信息對(duì)估計(jì)值進(jìn)行最優(yōu)估計(jì)，具備快速、高效的遞推特性，建議應(yīng)用于航空全張量磁力梯度實(shí)時(shí)測(cè)量數(shù)據(jù)處理，也可為其他航空地球物理實(shí)時(shí)測(cè)量數(shù)據(jù)處理提供參考。