基于Hartley變換的地磁場延拓技術

武立華,劉志海，孟 霆，黃 玉

(哈爾濱工程大學 理學院,黑龍江 哈爾濱 150001)

1 磁異常數據的Hartley變換

Hartley變換是一種類似于傅里葉變換的實數域積分變換方法， Hartley變換在圖像處理、模式識別等許多領域內已被廣泛利用[1-5]，但很少用于位場延拓處理中. 本文將Hartley變換應用到磁異常延拓積分迭代法中來提高運算效率.

設磁場場源位于平面z=0之下，z軸正向豎直向下，z<0空間中的磁異常分布f(x,y,z)是調和函數，z=0平面上的磁異常分布為已知觀測數據，求解z<0空間中的磁異常分布函數f(x,y,z). 根據狄利克雷邊值問題進行求解，得到磁異常延拓積分表達式為

(1)

令

將(1)式轉化為卷積形式：

f(x,y,z)=f(ξ,η,0)*φ(ξ,η).

(2)

由文獻[6]可得φ(ξ,η)的Hartley變換結果為

(3)

根據Hartley變換的卷積性質可以得到對應的Hartley變換形式：

H(u,v,z)=H(u,v,0)·φH(u,v),

(4)

式中，H(u,v,z)表示所求解f(x,y,z)的Hartley變換形式，H(u,v,0)表示已知磁異常分布數據f(x,y,0)的Hartley變換形式. 對(4)式進行Hartley逆變換可以得到所求平面的磁異常分布：

φH(u,v)cas (ux)cas (vy)dudv.

(5)

2 模型數據仿真

使用球形磁體模型計算得到z=195 m平面上的磁異常理論值,如圖1所示.

分別通過傳統的波數域迭代法和Hartley變換迭代法將z=0平面上的磁異常數據向下延拓至z=195 m平面. 圖2即為2種方法所得到的延拓結果的磁異常分布等值線圖.

(a)利用Hartley變換的迭代法

(b)傳統快速傅里葉變換迭代法圖2 利用2種迭代法向下延拓195 m后的磁異常分布

從圖2可以看出，經過2種迭代法處理后所得數據分布幾乎完全相同，這是由于在數據處理方面，Hartley變換和傅里葉變換的作用是一致的. 為了體現這種相似性，也為了更好地描述向下延拓的穩定性，在此引入3種誤差指標：絕對平均誤差、平均相對誤差和均方根誤差. 對圖2中的2種延拓結果分別進行誤差計算，2種迭代法的誤差相同：絕對平均誤差為0.31 nT，平均相對誤差為0.45%，均方根誤差為2.5 nT，但Hartley變換迭代法用時2.035 4 s，傅里葉變換法用時3.254 6 s，因此Hartley變換迭代法效率更高.

3 實驗檢驗

文中所用真實數據下載自美國國家地球物理數據中心(NGDC)，其原始數據為一定范圍平面內的磁異常數據. 設其原始平面為z=0平面，其磁異常數據等值線圖如圖3所示. 設定網格間距為100 m，計算網格為512×512. 為了進行向下延拓算法測試，首先利用向上延拓算法將原始數據平面向上延拓1 000 m，得到z=-1 000 m平面的磁異常分布數據如圖4所示.

圖3 z=0觀測平面理論磁異常分布

圖4 向上延拓得到的z=-1 000 m平面磁異常分布

接下來進行s=1，n=100的迭代運算，分別使用Hartley變換迭代法和傳統迭代法將z=-1 000 m平面的磁異常數據向下延拓1 000 m，得到如圖5所示的z=0平面上的延拓數據，將其和觀測平面理論數據(即原始平面數據)進行對比并計算誤差，同時對運算時間進行統計，2種迭代法的誤差相同：絕對平均誤差為7.2 nT，平均相對誤差為0.25%，均方根誤差為17.1 nT，Hartley變換迭代法用時4.286 3 s，傅里葉變換法用時10.255 6 s.

由圖5可知，通過2種方法進行向下延拓后的結果與原始數據基本一致. 此外3種誤差的值都很小，在實際需求可接受的范圍內，但相比理論實驗誤差有所增大，這是由于真實數據中存在噪聲而產生的，但2種方法對噪聲還是有一定的抗干擾能力的. 對比2種方法的運算時間可以發現，基于Hartley變換的迭代法效率要高出很多，因此優于傳統迭代法.

(a)Hartley變換迭代法延拓結果

(b)傳統迭代法延拓結果圖5 利用2種迭代法向下延拓1 000 m后得到的磁異常分布

4 結束語

通過理論模型和真實數據的仿真測試對Hartley變換算法進行檢驗，并利用絕對平均誤差、平均相對誤差和均方根誤差3種誤差統計指標對算法的誤差進行了統計，從而對Hartley變換的迭代法和傳統快速傅里葉變換迭代法進行了比較，證明利用Hartley變換的磁異常向下延拓迭代法的可行性和高效性.