本發(fā)明涉及的是一種地磁位場數(shù)據(jù)處理方法，具體涉及到一種平面地磁異常場一步向上延拓平面模量梯度場的遞歸余弦變換法。

背景技術(shù)：

地磁延拓是一種重要的地磁位場數(shù)據(jù)處理方法，在地磁匹配定位、地磁觀測及其資料解釋等方面具有重要應(yīng)用。地磁基準(zhǔn)圖是實(shí)現(xiàn)地磁匹配定位的數(shù)據(jù)基礎(chǔ)。為進(jìn)行地磁匹配定位，一般需要獲得航行深度面上的地磁基準(zhǔn)圖。不同航行器的航行高度各不相同，通過直接測量方式得到各航行高度上的地磁數(shù)據(jù)是不現(xiàn)實(shí)的。大面積的地磁測量值多數(shù)為地磁異常場數(shù)據(jù)，地磁模量梯度測量雖已開始，但匹配定位需要大范圍導(dǎo)航區(qū)域的基準(zhǔn)圖，現(xiàn)有的地磁模量梯度測量數(shù)據(jù)難以滿足這一要求。根據(jù)已有的航測地磁數(shù)據(jù)或者海測地磁數(shù)據(jù)，通過延拓和導(dǎo)數(shù)的同時計算一步獲得地磁模量梯度場，是獲得航行面地磁模量梯度場的有效途徑之一。向上延拓還可融合不同高度的航磁數(shù)據(jù)，使之歸化到同一高度面。地磁觀測及其資料解釋方面，利用向上延拓可以進(jìn)行區(qū)域場與局部場、深源場與淺源場的分離。

地磁模量梯度場具有以下優(yōu)點(diǎn)：第一，地磁模量梯度與不同磁傳感器間的姿態(tài)無關(guān)，降低了多個傳感器測量軸配置精度的要求；第二，便于根據(jù)載體磁場分布特性，選擇載體磁場模量相等的地點(diǎn)合理布置磁傳感器，通過求取模量梯度消除載體磁場等的影響。第三，地磁模量梯度場有三個分量，可利用的信息源較地磁異常場多；第四，模量梯度場能有效地突出淺源異常。

cordelll和grauchv指出因有限截斷和離散采樣，非周期函數(shù)的離散傅里葉變換會導(dǎo)致頻譜的“零頻”過小而高頻過高；為此提出先將觀測數(shù)據(jù)拓展成周期函數(shù)，再進(jìn)行離散傅里葉變換。在拓展之前采用等效源方法對觀測數(shù)據(jù)進(jìn)行擴(kuò)邊(cordelll,grauchv.reconciliationofthediscreteandintegralfouriertransforms,geophysics,1982,47(2):337-243)。ricardy和blakelyrj對邊界效應(yīng)進(jìn)行了較深入分析，提出了一種經(jīng)驗方法提高頻譜的計算精度(amethodtominimizeedgeeffectsintwo-dimensionaldiscretefouriertransforms,geophysics,1988,53(8):1113-1117)。郭華等人利用理論公式推導(dǎo)證明實(shí)測梯度數(shù)據(jù)進(jìn)行延拓轉(zhuǎn)換處理的可行性，研究了地磁梯度延拓后的數(shù)據(jù)解釋規(guī)律，闡明了梯度數(shù)據(jù)向上延拓的理論意義(郭華，王平，朱春華，杜穎.向上延拓對航磁梯度數(shù)據(jù)的影響及其規(guī)律研究，地球物理學(xué)進(jìn)展，2015，30(3):1214-1223)。王小兵等人利用基于傅里葉變換的頻率域迭代法，測試了實(shí)際基準(zhǔn)圖的延拓精度和對誤差的方法效果(王小兵，張金生，王哲，王仕成.地磁匹配中位場頻率域延拓方法研究，探測與控制學(xué)報，2009,31(增刊)：29-33)。

目前的航空或海洋測量數(shù)據(jù)主要是地磁異常值，其向上延拓方法也多針對地磁異常場，如何由平面地磁異常場一步向上延拓獲得平面地磁場模量梯度場的問題還較少涉及。基于傅里葉變換的平面地磁異常向上延拓平面地磁模量梯度場存在兩種兩步法，一種是先由觀測平面上的地磁異常場轉(zhuǎn)換模量梯度場，再由模量梯度場向上延拓得到延拓平面上的模量梯度場；另一種先由觀測平面上的地磁異常場向上延拓得到延拓面上的地磁異常場，再由延拓面上的地磁異常場轉(zhuǎn)換得到模量梯度場。但這兩種兩步法都會使用兩步基于傅里葉變換的波譜變換與反變換，增加了這種算法的計算時間，不利于向上延拓與位場轉(zhuǎn)換的實(shí)時性。另外，傅里葉變換通常要求數(shù)據(jù)序列的空間周期性，而實(shí)際上數(shù)據(jù)序列兩個邊界值往往不相等。因此，使用傅里葉變換進(jìn)行位場延拓或轉(zhuǎn)換都會在邊界形成一個跳躍，增加了高頻分量，形成gibbs邊界效應(yīng)。余弦變換可看作實(shí)偶序列的傅里葉變換，變換的輸入和輸出都為實(shí)數(shù)。與傅里葉變換相比，余弦變換具有更高的能量壓縮性能，在一階馬爾科夫過程中依據(jù)最小均方誤差原則是最接近karhunen-loeve變換性能的，可減小gibbs邊界效應(yīng)。

技術(shù)實(shí)現(xiàn)要素：

本發(fā)明的目的在于提供一種能簡化位場延拓與轉(zhuǎn)換的步驟，提高一步位場延拓與轉(zhuǎn)換的精度的平面地磁異常場一步向上延拓平面模量梯度場的遞歸余弦變換法。

本發(fā)明的目的是這樣實(shí)現(xiàn)的：

步驟1、選取網(wǎng)格化參數(shù)，對平面地磁異常場測量數(shù)據(jù)進(jìn)行預(yù)處理，剔除超出合理范圍的粗大測量誤差；

步驟2、利用徑向基函數(shù)插值法，將觀測平面上的離散地磁異常場測量數(shù)據(jù)進(jìn)行網(wǎng)格化，獲得觀測平面上的網(wǎng)格化地磁異常測量數(shù)據(jù)δt(nx,ny,z0)，nx＝1,2,…,mx，ny＝1,2，…,my，觀測平面方程為z＝z0，向上延拓平面方程為z＝z0-δz，z軸垂直向下；

步驟3、先后對觀測平面上的地磁異常場進(jìn)行行向和列向的一維遞歸余弦變換，這兩個方向的計算次序可以交換，得到地磁異常場的二維余弦變換波譜δt^c(ux,uy,z0)；

步驟4、按式由觀測平面上的地磁異常場波譜δt^c(ux,uy,z0)計算延拓平面上的地磁模量梯度場的波譜和式中，r1wc(ux,uy,δz)為余弦變換下的由觀測平面地磁異常場計算延拓平面地磁模量梯度場的一步波數(shù)域轉(zhuǎn)換因子，其表達(dá)式為

步驟5、先后對進(jìn)行行向的移相π/2的一維遞歸反余弦變換和列向的一維遞歸反余弦變換，定義為行向正弦、列向余弦反變換，這兩個方向的計算次序可以交換，得到延拓平面上的地磁模量梯度場的水平x分量δtx(nx,ny,z0-δz)；

步驟6、先后對進(jìn)行行向的一維遞歸反余弦變換和列向的移相π/2的一維遞歸反正弦變換，定義為行向余弦、列向正弦反變換，這兩個方向的計算次序可以交換，得到延拓平面上的地磁模量梯度場δty(nx,ny,z0-δz)的水平y(tǒng)分量；

步驟7、先后對進(jìn)行行向和列向的一維遞歸反余弦變換，這兩個方向的計算次序可以交換，得到延拓平面上的地磁模量梯度場δtz(nx,ny,z0-δz)的垂直z分量。

本發(fā)明的一種平面地磁異常場一步向上延拓平面模量梯度場的遞歸余弦變換法。先進(jìn)行平面地磁異常場測量數(shù)據(jù)預(yù)處理，剔除粗大誤差；利用徑向基函數(shù)插值法將觀測平面上不規(guī)則離散分布的地磁異常場測量數(shù)據(jù)進(jìn)行網(wǎng)格化，獲得網(wǎng)格化的地磁異常測量數(shù)據(jù)。再利用遞歸余弦變換將網(wǎng)格化的地磁異常測量數(shù)據(jù)轉(zhuǎn)換為余弦變換波譜；利用一步波數(shù)域轉(zhuǎn)換因子，由地磁異常余弦變換波譜計算地磁模量梯度場各分量的波譜。最后對水平x分量的波譜進(jìn)行行向移相π/2的遞歸反余弦變換，得到延拓平面上的地磁模量梯度場的水平x分量；對水平y(tǒng)分量的波譜進(jìn)行列向移相π/2的遞歸反余弦變換，得到延拓平面上的地磁模量梯度場的水平y(tǒng)分量；對垂直z分量的波譜進(jìn)行遞歸反余弦變換，得到延拓平面上的地磁模量梯度場的垂直z分量。其中對于并行計算而言，步驟5、6和7是可以進(jìn)行的。

本發(fā)明提出了一種基于遞歸余弦變換與反變換(含移相π/2的遞歸余弦反變換)的平面地磁異常場一步向上延拓平面地磁模量梯度場的方法，推導(dǎo)出余弦變換下的由觀測平面地磁異常場計算延拓平面地磁模量梯度場的一步波數(shù)域轉(zhuǎn)換因子。利用遞歸余弦變換與反變換(含移相π/2的遞歸余弦反變換)，由觀測平面上的地磁異常場一步計算延拓平面上的地磁模量梯度場。由于每一步位場轉(zhuǎn)換都有計算時間的消耗和計算誤差的存在，因此相比于兩步法，一步法的計算時間和計算誤差都得以減少，也簡化了計算流程。遞歸余弦變換與反變換可以適應(yīng)任意長度的數(shù)據(jù)序列，快速傅里葉變換與反變換通常要求數(shù)據(jù)序列的長度是基數(shù)的整數(shù)次冪，補(bǔ)零有時會增大不少的計算量。

本發(fā)明的方法能由平面地磁異常場直接地一步向上延拓，同時得到延拓平面上的地磁模量梯度場的三個分量；解決現(xiàn)有技術(shù)由觀測平面上的地磁異常場獲得向上延拓面上的地磁模量梯度場，需進(jìn)行向上延拓和地磁異常場轉(zhuǎn)換地磁模量梯度場的兩步轉(zhuǎn)換問題；兩步位場轉(zhuǎn)換通常需要兩步波譜變換和反變換，本方法僅需一步波譜變換和反變換，簡化了位場延拓與轉(zhuǎn)換的步驟，減少了計算時間；而且本方法采用遞歸余弦變換與反變換，在保證計算量小的同時，相比于傅里葉變換與反變換減小gibbs邊界效應(yīng)，提高了這種一步位場延拓與轉(zhuǎn)換的算法精度。

本發(fā)明與現(xiàn)有技術(shù)比較具有以下優(yōu)點(diǎn)：提出的一種平面地磁異常場一步向上延拓垂向模量梯度場的遞歸余弦變換法，具有計算量小及算法精度較高等優(yōu)點(diǎn)，解決現(xiàn)有技術(shù)由觀測平面上的地磁異常場獲得向上延拓平面上的地磁模量梯度場的三分量，需進(jìn)行向上延拓和地磁異常場轉(zhuǎn)換模量梯度場的兩步運(yùn)算問題；同時本發(fā)明僅需一步波譜變換和反變換就可同時得到延拓平面上的地磁模量梯度場的三分量，簡化了位場處理的步驟，減少了位場處理的計算量；而且在進(jìn)行波譜變換和反變換過程中采用遞歸余弦變換與反變換，在保證算法計算量小的同時，減小gibbs邊界效應(yīng)，提高了由平面地磁異常場一步向上延拓平面地磁模量梯度場的轉(zhuǎn)換精度。遞歸余弦變換與反變換可以適應(yīng)任意長度的數(shù)據(jù)序列，快速傅里葉變換與反變換通常要求數(shù)據(jù)序列的長度是基數(shù)的整數(shù)次冪，補(bǔ)零有時會增大不少的計算量。

附圖說明

圖1是平面地磁異常場一步向上延拓平面地磁模量梯度場的余弦變換法的方框圖。

圖2是平面地磁異常場一步向上延拓平面地磁模量梯度場的余弦變換法的流程圖。

圖3(a)-圖3(i)是球體模型數(shù)據(jù)對應(yīng)的傅里葉變換與余弦變換一步向上延拓模量梯度圖。其中：圖3(a)δtx(z0-δz)；圖3(b)圖3(c)圖3(d)δty(z0-δz)；圖3(e)圖3(f)圖3(g)δtz(z0-δz)；圖3(h)圖3(i)

圖4(a)-圖4(i)是球體和長方體混合模型數(shù)據(jù)對應(yīng)的傅里葉變換與余弦變換一步向上延拓模量梯度圖。其中：圖4(a)δtx(z0-δz)；圖4(b)圖4(c)圖4(d)δty(z0-δz)；圖4(e)圖4(f)圖4(g)δtz(z0-δz)；圖4(h)圖4(i)

具體實(shí)施方式

下面結(jié)合附圖舉例對本發(fā)明進(jìn)行詳細(xì)描述：

步驟1、平面地磁異常場測量數(shù)據(jù)的預(yù)處理。設(shè)定測量誤差及地磁異常值的合理范圍，對觀測平面上的地磁異常場測量數(shù)據(jù)超出這一合理范圍的粗大誤差數(shù)據(jù)予以剔除。

步驟2、平面地磁異常場測量數(shù)據(jù)的網(wǎng)格化處理。選取網(wǎng)格化參數(shù)，利用徑向基函數(shù)插值法，將觀測平面上離散的地磁異常場測量數(shù)據(jù)進(jìn)行網(wǎng)格化，獲得網(wǎng)格化的地磁異常測量數(shù)據(jù)。

徑向基函數(shù)插值法是一種精確度高的插值法，適用于對大量點(diǎn)數(shù)據(jù)進(jìn)行插值計算，具有較高的預(yù)測精度，能較好地反映數(shù)據(jù)的變化；它是利用基函數(shù)確定已知數(shù)據(jù)點(diǎn)到內(nèi)插網(wǎng)格節(jié)點(diǎn)的最佳權(quán)重。待插值點(diǎn)的徑向基函數(shù)平面插值表達(dá)式為

式中，(x,y)為插值點(diǎn)的坐標(biāo)，f0(x,y)＝c1xx+c1yy+c0，(xj,yj)為采樣點(diǎn)的坐標(biāo)，λj為對應(yīng)于每一個采樣點(diǎn)的權(quán)值，為徑向基函數(shù)，||·||為歐幾里德范數(shù)。選取薄板樣條函數(shù)作為徑向基函數(shù)，令d·j＝||x-xj,y-yj||，則有

由式(1)得到，采樣點(diǎn)的已知值所給定的方程約束條件為

設(shè)f(x,y)具有二次連續(xù)導(dǎo)數(shù)，則能量函數(shù)表示為：

最小化能量函數(shù)，得到所需滿足的正交條件為

為了求解式(1)中的系數(shù)c0、c1x、c1y及權(quán)值λj，聯(lián)立式(3)和式(5)得到關(guān)于系數(shù)c0、c1x、c1y及權(quán)值λj的線性方程組為

線性方程組的系數(shù)矩陣是一個對稱正定矩陣。求解得到系數(shù)c0、c1x、c1y及權(quán)值λj，將其代入式(1)得到觀測面上的網(wǎng)格化地磁異常場數(shù)據(jù)，nx＝1,2,…,mx，ny＝1,2,…,my，觀測面方程為z＝z0，z軸垂直向下。

步驟3、利用遞歸余弦變換由δt(nx,ny,z0)計算平面地磁異常場的余弦變換波譜定義

步1)按式(8)，由δt(ux,uy,z0)計算

式中，ux＝nx，θux＝(ux-1)π/mx，與由遞歸公式(9)確定。

式中，mx＝1,2,…,kx，kx＝[(mx+1)/2]，符號[·]表示取整。

步2)按式(11)，由計算δt^c(ux,uy,z0)，

式中，uy＝ny，θuy＝(uy-1)π/my，與由遞歸公式(12)確定。

式中，my＝1,2,…,ky，ky＝[(my+1)/2]，符號[·]表示取整。

步驟4、分別按式(14)、(15)和(16)，由觀測面上的地磁異常場的余弦變換波譜δt^c(ux,uy,z0)計算延拓平面上的地磁模量梯度場三分量的波譜和

步驟5、對進(jìn)行行向的移相π/2的一維遞歸反余弦變換和列向的一維遞歸反余弦變換，得到延拓平面上的地磁模量梯度場的x分量δtx(nx,ny,z0-δz)。

步1)按式(17)，由計算

式中，θnx＝(nx-1/2)π/mx，由遞歸公式(18)確定。

步2)按式(19)，由計算δtx(nx,ny,z0-δz)，

式中，θny＝(ny-1)π/my，與由遞歸公式(20)確定。

步驟6、對進(jìn)行行向的一維遞歸反余弦變換和列向的移相π/2的一維遞歸反正弦變換，得到延拓平面上的地磁模量梯度場的y分量δty(nx,ny,z0-δz)；

步1)按式(21)，由計算

式中，θnx＝(nx-1/2)π/mx，和由遞歸公式(22)確定。

步2)按式(23)，由計算δty(nx,ny,z0-δz)，

式中，θny＝(ny-1)π/my，由遞歸公式(24)確定。

步驟7、對進(jìn)行行向和列向的一維遞歸反余弦變換，得到延拓平面上的地磁模量梯度場的z分量δtz(nx,ny,z0-δz)。

步1)按式(25)，由計算

式中，與由遞歸公式(26)確定。

步2)按式(27)，由計算δtz(nx,ny,z0-δz)，

式中，與由遞歸公式(28)確定。

利用遞歸余弦變換與反變換由平面地磁異常場一步向上延拓，得到平面地磁模量梯度場三分量，平面地磁異常場一步向上延拓平面地磁模量梯度場的余弦變換法的方框圖如圖1所示。平面地磁異常場一步向上延拓平面地磁模量梯度場的余弦變換法的算法流程圖如圖2所示。

為直接反映平面地磁異常場一步向上延拓平面地磁模量梯度場的三分量的位場轉(zhuǎn)換效果，定義一個無量綱的相對誤差指標(biāo)和

傅里葉變換與反變換下平面地磁模量梯度場j分量(j＝x,y,z)的位場轉(zhuǎn)換相對誤差為

式中，δtj(z0-δz)為平面地磁模量梯度場j分量的真實(shí)值，為使用傅里葉變換與反變換得到的平面地磁模量梯度場j分量的延拓值，符號表示求平均值，符號||表示求絕對值。

余弦變換與反變換下平面地磁模量梯度場j分量的位場轉(zhuǎn)換相對誤差為

式中，為使用余弦變換與反變換得到的平面地磁模量梯度場j分量的延拓值。

實(shí)驗一：觀測平面為z＝z0＝0，觀測區(qū)域的x和y坐標(biāo)范圍都為[-5110m,5110m]，網(wǎng)格數(shù)為256×256，向上延拓的高度為δz＝-300m。六個球體磁源的位置分別為[0,0,500m]、[2000m,800m,400m]、[1800m,-1600m,500m]、[2000m,-1800m,600m]、[-1000m,2000m,500m]和[1500m,0,1000m]，磁化強(qiáng)度分別為1.6×10⁵a/m、1.4×10⁵a/m、1.3×10⁵a/m、1.5×10⁵a/m、1.2×10⁵a/m和1.8×10⁵a/m，球半徑分別為10m、12m、8m、15m、20m和18m，磁傾角分別45°、60°、75°、45°、60°和25°，磁偏角分別為45°、35°、15°、30°、45°和45°。地磁場的磁傾角和磁偏角分別為60°和30°。沒有磁測噪聲時，延拓平面上的地磁模量梯度三分量的模型值、基于傅里葉變換與反變換的延拓值和基于余弦變換與反變換的延拓值如圖3所示，計算得到兩種變換下的相對誤差分別是和這在理論上表明基于余弦變換與反變換的延拓地磁模量梯度精度高于基于傅里葉變換與反變換的精度。

觀測面的磁測噪聲為零均值、方差為10nt的高斯白噪聲時進(jìn)行50次蒙特卡洛仿真實(shí)驗得到延拓的平均相對誤差和球體模型數(shù)據(jù)的一步向上延拓地磁模量梯度的仿真實(shí)驗表明，基于余弦變換與反變換的延拓地磁模量梯度精度高于基于傅里葉變換與反變換的精度。

實(shí)驗二：在上述仿真實(shí)驗區(qū)域添加兩個長方體磁源。第一個長方體的中心位置為[0,400m,800m]，長、寬和高分別為30m、20m和15m，磁化強(qiáng)度為0.02a/m，磁傾角和磁偏角分別為45°和55°。第二個長方體的中心位置為[-1500m,1000m,600m]，長、寬和高分別為25m、18m和12m，磁化強(qiáng)度為0.012a/m，磁傾角和磁偏角分別為60°和40°。沒有磁測噪聲時，延拓平面上的地磁模量梯度三分量的模型值、基于傅里葉變換與反變換的延拓值和基于余弦變換與反變換的延拓值如圖4所示，計算得到兩種變換下的相對誤差和這在理論上表明基于余弦變換與反變換的延拓地磁模量梯度精度高于基于傅里葉變換與反變換的精度。

觀測面的磁測噪聲為零均值、方差為10nt的高斯白噪聲時進(jìn)行50次蒙特卡洛仿真實(shí)驗得到延拓的平均相對誤差和球體和長方體模型數(shù)據(jù)的一步向上延拓地磁模量梯度的仿真實(shí)驗同樣表明，基于余弦變換與反變換的延拓地磁模量梯度精度高于基于傅里葉變換與反變換的精度。

本發(fā)明有益效果說明如下：

由于步驟5、步驟6和步驟7對應(yīng)的三次波譜反變換在使用并行計算時可以同時進(jìn)行，因此可以認(rèn)為平面地磁異常場轉(zhuǎn)換向上延拓平面地磁模量梯度場的一步法共需要一次波譜變換和一次波譜反變換。而通常的兩種兩步法，需要兩次波譜變換和兩次波譜反變換。分析比較表明了使用同種波譜變換條件下，一步法的計算量約為兩步法計算量的一半。球體模型數(shù)據(jù)和球體與長方體混合模型數(shù)據(jù)仿真比較了基于傅里葉變換和余弦變換的一步向上延拓法的計算結(jié)果，結(jié)果表明余弦變換法減小gibbs邊界效應(yīng)，提高了算法精度。遞歸余弦變換與反變換對數(shù)據(jù)序列的長度沒有要求，可應(yīng)用于任意長度的數(shù)據(jù)序列。