“多維”多尺度重磁位場數(shù)據(jù)融合方法

劉福香， 王萬銀， 紀(jì)曉琳， 熊盛青

1 長安大學(xué)重磁方法技術(shù)研究所， 長安大學(xué)地質(zhì)工程與測繪學(xué)院， 長安大學(xué)西部礦產(chǎn)資源與地質(zhì)工程教育部重點實驗室， 西安 710054 2 中國自然資源航空物探遙感中心, 北京 100083 3 自然資源部航空地球物理與遙感地質(zhì)重點實驗室, 北京 100083

0 引言

重、磁勘探以其經(jīng)濟快速、綠色環(huán)保，覆蓋范圍廣、橫向分辨率高等特點在解決地質(zhì)問題中得到了廣泛應(yīng)用.為了更好地解決某一地質(zhì)問題，往往在不同時期采集了不同尺度(比例尺、精度)、不同“維度”(衛(wèi)星、航空、地面、海洋及井中等)的重、磁位場數(shù)據(jù).然而，由于存在觀測精度和觀測位置及數(shù)據(jù)基準(zhǔn)等差異，同一研究區(qū)采集的位場數(shù)據(jù)存在著精度差異、高度差異及數(shù)據(jù)基準(zhǔn)差異等問題.如何將這些不同精度、不同比例尺、不同“維度”和不同基準(zhǔn)的重、磁位場數(shù)據(jù)融合成為統(tǒng)一基準(zhǔn)、同一“維度”(后簡寫為維度)的位場數(shù)據(jù)就成為重、磁位場數(shù)據(jù)融合的核心問題之一.根據(jù)位場數(shù)據(jù)融合中數(shù)據(jù)觀測比例尺和觀測維度的差異，本文將位場數(shù)據(jù)融合分為單尺度和多尺度、單維和多維融合問題.單尺度指相同網(wǎng)格距、相同比例尺和相同精度的位場數(shù)據(jù)，多尺度指不同網(wǎng)格距、不同比例尺或不同精度的位場數(shù)據(jù)；單維指同一個連續(xù)觀測面的位場數(shù)據(jù)，多維指不同高度連續(xù)觀測面或間斷觀測面的位場數(shù)據(jù).

單尺度重、磁位場數(shù)據(jù)具有相同的觀測位置和觀測精度，不同來源的單尺度位場數(shù)僅存在數(shù)據(jù)基準(zhǔn)的差異，其融合主要考慮數(shù)據(jù)基準(zhǔn)的統(tǒng)一；而多尺度數(shù)據(jù)由于存在觀測精度或觀測比例尺的差異，需要在單尺度數(shù)據(jù)基準(zhǔn)統(tǒng)一的基礎(chǔ)上考慮精度校正.

單尺度數(shù)據(jù)融合方法主要有加權(quán)平均和主成分分析法等加權(quán)類方法(薛重生等，1997；張麗莉等，2003)，其在沒有基準(zhǔn)差異的單尺度數(shù)據(jù)融合中計算效果良好，若數(shù)據(jù)基準(zhǔn)存在較大差異，不考慮背景差校正的加權(quán)融合則不滿足位場數(shù)據(jù)融合要求.多尺度數(shù)據(jù)融合方法主要分為兩大類，第一類是插值方法，不同精度的數(shù)據(jù)同時進(jìn)行插值計算，主要方法有Kriging插值、加權(quán)最小曲率插值、Shepard插值、樣條插值等(Beylkin and Cramer，2002；劉翔等，2016；劉善偉等，2018a；王虎彪等，2018；柯寶貴等，2018；支澳威和陳華根，2019)，此類方法與單尺度數(shù)據(jù)中的加權(quán)計算方法思想類似，能夠考慮不同精度數(shù)據(jù)的差異，但直接插值計算結(jié)果會受到多尺度數(shù)據(jù)基準(zhǔn)差影響，不同精度數(shù)據(jù)邊界處形成梯級帶；第二類是由不同分辨率的圖像融合方法延伸而來，包括基于小波變換(Ranchin and Wald，1993；Deighan and Watts，1997；郭達(dá)志和方濤，1997；Kuroishi and Keller，2005；汪海洪，2005；李雙成等，2006；Douma and De Hoop，2007；Herrmann et al.，2008；吳健生和劉苗，2008；王海江，2012；楊揚，2013；畢奔騰等，2016；譚沛森，2018；劉善偉等，2018b；陳勐韜等，2020)和Curvelet變換(Candes and Donoho，2005；Nencinia et al.，2007；張恒磊和劉天佑，2011)的融合方法，這類方法將數(shù)據(jù)分解為不同頻段的信息，通過對多來源數(shù)據(jù)頻段的選擇和重新組合完成數(shù)據(jù)背景差的校正和不同精度數(shù)據(jù)的融合，其在不同分辨率的圖像融合中應(yīng)用效果良好，但由于不同精度的位場數(shù)據(jù)都包含全段頻率的信息，進(jìn)行頻率分解再組合的方法會損失部分高精度數(shù)據(jù)細(xì)節(jié)信息，不符合位場數(shù)據(jù)優(yōu)選高精度數(shù)據(jù)的原則.

單維數(shù)據(jù)不存在觀測位置差異，其融合按照上述分類直接考慮單尺度或多尺度融合方法，而多維數(shù)據(jù)則需在單尺度融合和多尺度融合的基礎(chǔ)上再考慮觀測空間位置的差異.在實際的數(shù)據(jù)處理工作中，有學(xué)者不考慮觀測位置的差異，將多維數(shù)據(jù)當(dāng)作單維數(shù)據(jù)進(jìn)行處理，這不符合位場數(shù)據(jù)隨空間高度位置變化的規(guī)律.

多維數(shù)據(jù)融合可以通過位場延拓方法實現(xiàn)維度統(tǒng)一后，再進(jìn)行上述單尺度和多尺度數(shù)據(jù)的融合計算(吳招才等，2018；余春平，2019)；而目前針對多維數(shù)據(jù)的融合方法主要分為統(tǒng)計法和解析法兩大類.統(tǒng)計法主要有最小二乘配置法，利用多維數(shù)據(jù)建立協(xié)方差函數(shù)后進(jìn)行系統(tǒng)偏差估計和精度估計，并在此基礎(chǔ)上逐漸發(fā)展了一些迭代算法、正則化方法等(Tscherning and Rapp，1974；Forsberg，1987；楊元喜，1996；Tscherning et al.，1998；Strykowski and Forsberg，1998；Tziavos et al.，1998；Olesen et al.，2002；Kern et al.，2003；Huang et al.，2006；翟振和和孫中苗，2010；劉曉剛等，2012；黃謨濤等，2013b；潘宗鵬等，2013；孫文等，2016；王灼華等，2017；趙啟龍，2017)，該方法能夠?qū)⒉煌^測面上的多種類型數(shù)據(jù)融合起來，但需以足夠分辨率的觀測數(shù)據(jù)為基礎(chǔ)，且經(jīng)驗協(xié)方差函數(shù)的選取復(fù)雜、大型矩陣的求解過程使得該方法適用性和穩(wěn)定性受到限制.許多學(xué)者致力于解析法的研究(Lehmann，1993；Sideris，1996；Li and Sideris，1997；寧津生等，1998；Kern et al.，2003；郝燕玲等，2007；Eicker，2008；Wittwer，2010；黃謨濤等，2013a，2015；王燚等，2015；吳懌昊等，2016；劉金釗等，2020)，該方法在迭代中利用不同觀測面上解析計算的位場數(shù)據(jù)殘差修正觀測數(shù)據(jù)或模型系數(shù)，引入分步平差、推估內(nèi)插、迭代下延、點質(zhì)量等技術(shù)，在大區(qū)域的重、磁力場計算中應(yīng)用廣泛，能夠省去統(tǒng)計法中求解大型矩陣逆的問題，在迭代解析計算中完成多維多尺度位場數(shù)據(jù)的融合，但迭代中修正模型數(shù)據(jù)的方法無法全部保留高精度數(shù)據(jù)，且仍存在迭代穩(wěn)定性的問題.

針對上述研究現(xiàn)狀，本文提出了一套適用于多維多尺度重、磁位場數(shù)據(jù)融合的實用性方法；對數(shù)據(jù)維度、基準(zhǔn)、尺度的差異，分別采用位場延拓、回歸分析、精度加權(quán)平均等方法來解決高度不同、基準(zhǔn)不同、精度不同或維度不同、尺度不同的重、磁位場數(shù)據(jù)融合問題，以期得到統(tǒng)一維度、基準(zhǔn)及尺度的重、磁位場數(shù)據(jù)，為進(jìn)一步處理、反演和解釋奠定數(shù)據(jù)基礎(chǔ).

1 多維多尺度重、磁位場數(shù)據(jù)融合方法

為了解決多維多尺度位場數(shù)據(jù)融合問題，需研究尺度統(tǒng)一方法、基準(zhǔn)差統(tǒng)一方法(單尺度、多尺度位場數(shù)據(jù)融合方法)和維度統(tǒng)一方法(多維單尺度位場數(shù)據(jù)融合方法).

單尺度數(shù)據(jù)融合方法，采用相關(guān)分析和回歸分析校正基準(zhǔn)差解決數(shù)據(jù)基準(zhǔn)差異；多尺度數(shù)據(jù)融合方法，采用不同精度數(shù)據(jù)加權(quán)解決觀測尺度和觀測精度的差異；多維多尺度位場數(shù)據(jù)融合方法，在以上兩類數(shù)據(jù)融合的基礎(chǔ)上采用位場延拓方法解決觀測維度的差異(圖1).

圖1 多維多尺度重、磁位場數(shù)據(jù)融合計算流程示意圖Fig.1 The schematic diagram of multi-dimension and multi-scale gravity and magnetic potential data fusion calculation process

1.1 維度統(tǒng)一方法

本文采用位場延拓實現(xiàn)多維多尺度重、磁位場數(shù)據(jù)融合中的維度統(tǒng)一.本次所使用的位場延拓方法以空間域積分-迭代法(徐世浙，2006)為基礎(chǔ)算法，將不穩(wěn)定的向下延拓問題轉(zhuǎn)換為迭代的穩(wěn)定向上延拓問題，每次迭代中向上延拓和偏差疊加更新過程是影響計算精度和效率的關(guān)鍵，這兩步驟也影響維度統(tǒng)一過程的計算精度和效率，因此成為進(jìn)行多維數(shù)據(jù)融合的基礎(chǔ).基于此，本文采用滑動窗口加曲面切片技術(shù)(劉芬和王萬銀，2019)提高計算效率，采用向上延拓濾波技術(shù)提高計算穩(wěn)定性.

(1)位場向上延拓

空間域由原始點(ξ,η,ζ)延拓至計算點(x,y,z)的基本公式(曾華霖，2005)為：

U(x,y,z)=

(1)

以上解析積分在程序?qū)崿F(xiàn)中使用的是數(shù)值積分的矩形積分公式，則(1)式變?yōu)椋?/p>

(2)

其中，Δξ和Δη分別為x和y方向的點距，H(x,y,z,ξi,ηj,ζij)為核函數(shù)，其值如(3)式，

H(x,y,z,ξi,ηj,ζij)=

本文采用滑動窗口和曲面切片的向上延拓計算技術(shù)(劉芬，2019)提高計算效率.根據(jù)位場延拓計算式(式(1))看出，距離越遠(yuǎn)的場源對計算點產(chǎn)生的異常影響越小.因此根據(jù)延拓高度在計算點一定窗口范圍內(nèi)計算，一般選擇20倍的延拓高度以上，且小于剖面長度的一半(劉芬等，2019)，可在保證精度的前提下減少計算量.

對于延拓高度相同的點(圖2中P點和Q點)，在相同的窗口大小內(nèi)，其核函數(shù)矩陣是相等的；平面到平面延拓(圖2由ΓA到ΓS)，可以計算得到窗口內(nèi)核函數(shù)矩陣，然后直接代入式(2)完成每個計算點的延拓.

圖2 向上延拓中滑動窗口和曲面切片示意圖Fig.2 The schematic diagram of sliding window and surface slice in upward continuation

若延拓面為曲面(圖2中ΓC)，用一組平面切片Slice將其切割后，采用平面的滑動窗口核函數(shù)計算方法，首先計算得到每個切片高度上的固定窗口內(nèi)的核函數(shù)矩陣，落在切片上的計算點(Q點和M點)直接代入，分布在切片之間的計算點(N點)的位場值UN則采用反距離加權(quán)插值得到，其加權(quán)公式(劉芬，2019)為：

(4)

其中，UM、UQ分別為M點和Q點的位場值，hMN、hQN為N點分別與M點、Q點在垂直方向的距離.

由此，通過在空間域向上延拓計算中加入滑動窗口和曲面切片技術(shù)可在保證精度的前提下大大減少計算量，提高多維重、磁位場數(shù)據(jù)融合計算的效率.

設(shè)計理論模型檢驗向上延拓計算方法技術(shù)的正確性和有效性.場源由五個直立六面體組成，空間位置和物性參數(shù)(表1)及其水平展布位置(圖3a)已知；模型基礎(chǔ)觀測曲面高度起伏為-1000～-2500 m(觀測曲面一，z坐標(biāo)向下為正，圖3b)，觀測平面為z=0 m，x方向坐標(biāo)范圍0～26000 m，y方向坐標(biāo)范圍0～26000 m，兩個方向網(wǎng)格間距都為100 m，在觀測曲面一和z=0 m平面上正演計算得到重力異常和磁力異常(圖3c、d、e、f).

本次測試z=0 m平面數(shù)據(jù)向上延拓至觀測曲面一(-1000～-2500 m)的計算精度和效率.分別進(jìn)行無窗口的全區(qū)域延拓計算、加滑動窗口的延拓計算以及滑動窗口結(jié)合曲面切片的延拓計算，分別統(tǒng)計了延拓計算至觀測曲面一的重力異常(圖4a、b、c)和磁力異常(圖4d、e、f)與理論重、磁力異常(圖3c、d)的均方根誤差、相對均方根誤差以及計算時間(表2，配置Intel(R) Core(i5) CPU 2.80 GHz，內(nèi)存12.0 G的計算機).由表2可以看出延拓計算的重、磁力異常與理論值的偏差都較小，滑動窗口和曲面切片的計算技術(shù)能夠保證延拓計算的精度；無窗口的全區(qū)域的延拓計算耗時最長，加滑動窗口的延拓計算耗時次之，滑動窗口結(jié)合曲面切片的延拓計算耗時最短，驗證了本次采用的滑動窗口結(jié)合曲面切片計算技術(shù)能夠提高計算效率.這為迭代法向下延拓的多維重、磁位場數(shù)據(jù)融合中維度統(tǒng)一處理奠定了計算基礎(chǔ).

表1 理論模型體空間位置坐標(biāo)及物性參數(shù)Table 1 The space position coordinates and physical property parameters of theoretical models

表2 理論模型向上延拓計算誤差統(tǒng)計表Table 2 The error table of upward continuation calculation of theoretical models

圖3 理論模型平面位置、觀測面一起伏及重、磁力異常等值線圖(a) 理論模型空間位置平面圖； (b) 觀測曲面一高程等值線圖； (c) 模型體在觀測曲面一上引起的重力異常等值線圖； (d) 模型體在觀測曲面一上引起的磁力異常等值線圖； (e) 模型體在z=0 m平面上引起的重力異常等值線圖； (f) 模型體在z=0 m平面上引起的磁力異常等值線圖.圖中黑色框線表示五個模型體的邊界位置.Fig.3 The plane position of theoretical models, observation surface One and its gravity and magnetic anomaly contour maps(a) The plane position of theoretical models; (b) The contour map of observation surface One elevation; (c) The contour map of gravity anomaly caused by the models on the observation surface One; (d) The contour map of magnetic anomaly caused by the models on the observation surface One; (e) The contour map of gravity anomaly caused by the models on the z=0 m plane; (f) The contour map of magnetic anomaly caused by the models on the z=0 m plane. The black solid lines represent the edge location of the five models.

圖5 積分-迭代法向下延拓示意圖Fig.5 The schematic diagram of integral-iterative method for downward continuation

(2)位場向下延拓

本文采用積分-迭代法(徐世浙，2006)，將向下延拓轉(zhuǎn)換為迭代法的向上延拓.如圖5所示，將觀測面ΓB上的位場值UB(x,y,z)向下延拓至ΓA主要進(jìn)行以下幾個步驟：

(5)

(6)

ΔU(n)′B(x′,y′,z′)=

(7)

其中濾波器H就是向上延拓的核函數(shù)(式(3))；相應(yīng)地，迭代中的校正更新公式(式(6))變?yōu)椋?/p>

圖6 觀測曲面二地形起伏Fig.6 The contour map of observation surface Two elevation

采用理論模型測試向下延拓方法技術(shù)的準(zhǔn)確性和精度.觀測曲面二與觀測曲面一坐標(biāo)范圍一致，z坐標(biāo)在觀測曲面一的基礎(chǔ)上抬升3000 m(z坐標(biāo)從-3000 m到-4500 m變化，圖6)，在其正演的重、磁力異常(圖7a、b)基礎(chǔ)上加入各自數(shù)據(jù)幅值5%的隨機噪聲，即重力異常數(shù)據(jù)噪聲均值為0 mGal、幅值為-0.049～0.055 mGal，磁力異常數(shù)據(jù)噪聲均值為0 nT、幅值為-27.488～27.888 nT，利用加噪后的數(shù)據(jù)(圖7c、d)向下延拓至不同的平面(z=0 m、z=-1000 m、z=-2000 m、z=-3000 m)，得到的向下延拓結(jié)果(圖8)，計算的均方根誤差和相對均方根誤差見表3.由圖8可以看出，向下延拓計算結(jié)果符合位場值隨空間變化的規(guī)律，隨著觀測面高度與模型體距離的增加，重、磁力異常的幅值減小.由表3可以看出，原始數(shù)據(jù)存在噪聲時，進(jìn)行低通濾波處理后向下延拓計算效果良好，各個延拓結(jié)果的計算誤差??；且隨著延拓距離的增加，向下延拓計算精度逐漸降低，這與位場隨空間變化規(guī)律相符合.

圖7 理論模型觀測曲面二重、磁力異常等值線圖(a) 重力異常； (b) 磁力異常； (c) 加噪聲后的重力異常； (d) 加噪聲后的磁力異常.Fig.7 The contour maps of gravity and magnetic anomaly of theoretical models on the observation surface Two(a) The gravity anomaly; (b) The magnetic anomaly; (c) The gravity anomaly that added noise; (d) The magnetic anomaly that added noise.

圖8 觀測曲面二向下延拓計算重、磁力異常等值線圖(a) z=0 m重力異常； (b) z=-1000 m重力異常； (c) z=-2000 m重力異常； (d) z=-3000 m重力異常； (e) z=0 m磁力異常； (f) z=-1000 m磁力異常； (g) z=-2000 m磁力異常； (h) z=-3000 m磁力異常.Fig.8 The contour map of downward continuation results from observation surface Two(a) The gravity anomaly on the z=0 m plane; (b) The gravity anomaly on the z=-1000 m plane; (c) The gravity anomaly on the z=-2000 m plane; (d) The gravity anomaly on the z=-3000 m plane; (e) The magnetic anomaly on the z=0 m plane; (f) The magnetic anomaly on the z=-1000 m plane; (g) The magnetic anomaly on the z=-2000 m plane; (h) The magnetic anomaly on the z=-3000 m plane.

表3 理論模型重、磁力異常向下延拓計算誤差統(tǒng)計表Table 3 The error table of gravity or magnetic anomaly downward continuation calculation of theoretical models

同時，為檢驗上述位場延拓計算方法在多維重、磁位場數(shù)據(jù)融合中的應(yīng)用效果，本文將向下延拓的計算結(jié)果(圖8)分別與原始理論重、磁力異常進(jìn)行融合(圖9a、c)，同時對比不進(jìn)行位場延拓直接“融合”的結(jié)果(即將多維數(shù)據(jù)當(dāng)作單維數(shù)據(jù)直接處理，圖9b、d)，并統(tǒng)計融合計算誤差(表4).結(jié)果顯示，在多維數(shù)據(jù)融合計算中，不考慮數(shù)據(jù)之間的維度差異直接進(jìn)行“融合”是不正確，而本文采用的位場延拓方法可以有效統(tǒng)一不同數(shù)據(jù)的維度差異，計算誤差小，是進(jìn)行多維多尺度位場數(shù)據(jù)融合的維度統(tǒng)一的有效方法.

表4 理論模型多維重、磁力異常融合計算誤差統(tǒng)計表Table 4 The error table of multi-dimension gravity or magnetic anomaly fusion of theoretical models

圖9 多維數(shù)據(jù)融合計算重、磁力異常立體圖(a) 多維重力異常進(jìn)行維度統(tǒng)一后的融合結(jié)果； (b) 多維重力異常當(dāng)作單維直接進(jìn)行處理的計算結(jié)果； (c) 多維磁力異常進(jìn)行維度統(tǒng)一后的融合結(jié)果； (d) 多維磁力異常當(dāng)作單維直接進(jìn)行處理的計算結(jié)果.Fig.9 The 3D surface map of multi-dimension gravity or magnetic anomalies fusion results(a) The fusion results of multi-dimension gravity anomalies after dimension unification; (b) The calculation results of multi-dimension gravity anomalies regarded as single-dimension data; (c) The fusion results of multi-dimension magnetic anomalies after dimension unification; (d) The calculation results of multi-dimension magnetic anomalies regarded as single-dimension data.

1.2 基準(zhǔn)差統(tǒng)一方法

本文采用回歸分析校正基準(zhǔn)的方法技術(shù)解決單尺度數(shù)據(jù)和多尺度數(shù)據(jù)的基準(zhǔn)差異問題.不同組數(shù)據(jù)間進(jìn)行回歸分析的前提是其存在顯著性相關(guān)，故需對重疊區(qū)的多來源數(shù)據(jù)進(jìn)行相關(guān)顯著性檢驗.

在單尺度位場數(shù)據(jù)融合中，由于不存在精度或比例尺的差異，在理論上選取任一組數(shù)據(jù)作為基準(zhǔn)數(shù)據(jù)都可，其余數(shù)據(jù)都向此基準(zhǔn)數(shù)據(jù)校正；而對于存在精度或比例尺差異的多尺度數(shù)據(jù)，須選擇精度最高或比例尺最大的一組數(shù)據(jù)作為校正的基準(zhǔn).

統(tǒng)計學(xué)領(lǐng)域中，相關(guān)分析是用來考查兩組地位平等且無因果關(guān)系的變量相互變化的關(guān)聯(lián)關(guān)系；回歸分析是研究一個因變量與一個或多個自變量之間數(shù)量關(guān)系的統(tǒng)計方法.應(yīng)用至位場數(shù)據(jù)融合中，相關(guān)分析中為兩組重、磁位場數(shù)據(jù)，例如不同時期采集的同一地區(qū)的維度相同、尺度相同但存在數(shù)據(jù)基準(zhǔn)差異的重力異常，記為gA、gB，采用Pearson相關(guān)計算公式，即(9)式，通過計算相關(guān)系數(shù)r判定兩組數(shù)據(jù)之間的線性相關(guān)程度(Pearson,1895)：

(9)

不同來源的重、磁位場數(shù)據(jù)存在的數(shù)據(jù)基準(zhǔn)差主要源自觀測中基點不同，因此在建立回歸模型過程中考慮構(gòu)建較簡單的回歸模型，如一元線性回歸或一元二次回歸即可擬合不同數(shù)據(jù)基準(zhǔn)的差異.求解方法是在建立合適的回歸模型條件下，以最小二乘法求解模型統(tǒng)計數(shù)，獲得回歸方程，從而依據(jù)其中一組作為基準(zhǔn)的數(shù)據(jù)對其他位場數(shù)據(jù)進(jìn)行校正.

例如，對于N元線性回歸，簡單回歸模型表示為(Galton，1886)：

gA=β0+β1gB+…+βNgB+ε,

(10)

模型中，gA是gB的函數(shù)部分加誤差項，函數(shù)部分反映了數(shù)據(jù)gB與gA之間的變化對應(yīng)關(guān)系，誤差項ε為隨機變量，反映了不同時期采集的位場數(shù)據(jù)gB和gA之間的線性關(guān)系之外的隨機干擾，β0、β1、…、βN為回歸模型的參數(shù)，由最小二乘估計得到.

(11)

(12)

根據(jù)這一基本原理，進(jìn)行多組數(shù)據(jù)之間的一元線性回歸分析的處理是選擇一個回歸的數(shù)據(jù)體gA(x,y)，依次用其他數(shù)據(jù)gB(x,y),gC(x,y),…,gM-1(x,y)與其進(jìn)行回歸校正，并建立簡單線性回歸模型形式為：

(13)

即只用數(shù)據(jù)本身進(jìn)行校正，得到每個數(shù)據(jù)體與選取的基準(zhǔn)數(shù)據(jù)的回歸系數(shù)后，進(jìn)行相應(yīng)的計算即完成了數(shù)據(jù)基準(zhǔn)的回歸分析校正.

圖10 觀測曲面三地形起伏Fig.10 The contour map of observation surface Three elevation

選取單尺度理論模型數(shù)據(jù)測試該回歸分析方法在消除基準(zhǔn)差中的正確性和有效性.理論模型場源幾何參數(shù)與物性參數(shù)見表1，為模擬不同時期采集的存在數(shù)據(jù)基準(zhǔn)差異的不同范圍的單尺度重力、磁力異常數(shù)據(jù)，設(shè)計兩個不同觀測曲面(觀測曲面一和觀測曲面三)，觀測曲面三(圖10)是在觀測曲面一(圖3b)的基礎(chǔ)上存在空白區(qū)，兩觀測曲面的重復(fù)區(qū)占總范圍的37%，設(shè)計觀測面三原始觀測異常(數(shù)據(jù)體2)和觀測面一原始觀測異常(數(shù)據(jù)體1)之間存在線性的基準(zhǔn)差，其關(guān)系式如式(14)所示：

(14)

對兩組重力異常(圖11a、b)和磁力異常(圖11c、d)進(jìn)行回歸分析的基準(zhǔn)校正，分別得到校正后重力異常和磁力異常(圖12a、b)和計算誤差(圖12c、d)，可以看到其校正后重、磁力異常形態(tài)和變化趨勢符合位場數(shù)據(jù)的特點，且形態(tài)和幅值均與理論重力異常(圖3c)和磁力異常(圖3d)吻合很好，計算誤差小；重力異常校正計算的均方根誤差小于2.008×10-4mGal，相對均方根誤差為0.0079%，磁力異常校正計算的均方根誤差小于0.135 nT，相對均方根誤差為0.0087%，計算精度高，說明本次采用的回歸分析能夠?qū)纬叨任粓鰯?shù)據(jù)基準(zhǔn)差進(jìn)行很好的校正.

多尺度理論模型場源幾何參數(shù)與物性參數(shù)見表1.為模擬多尺度數(shù)據(jù)特點，在觀測曲面一的基礎(chǔ)上，改變數(shù)據(jù)間隔為50 m并設(shè)置空白區(qū)形成觀測曲面四(圖13a)，改變數(shù)據(jù)間隔為200 m形成觀測曲面五(圖13e)，數(shù)據(jù)基準(zhǔn)差關(guān)系式如式(14)，同時加入不同大小的隨機噪聲模擬兩者觀測精度的差異，重力異常數(shù)據(jù)噪聲均值為0 mGal、幅值為-0.028～0.027 mGal(圖13b)和-0.048～0.045 mGal(圖13c)，磁力異常數(shù)據(jù)噪聲均值為0 nT、幅值為-5.826～5.944 nT(圖13f)和-11.861～12.982 nT(圖13g)；兩組數(shù)據(jù)覆蓋范圍不同，存在數(shù)據(jù)重復(fù)區(qū)(數(shù)據(jù)重復(fù)區(qū)占總區(qū)域37%).

對重力異常、磁力異常分別進(jìn)行回歸分析校正得到原觀測面上背景差校正后的重、磁力異常(圖13d、h)，可以看到其與理論重力異常(圖3c)和磁力異常(圖3d)的形態(tài)和幅值均吻合很好，基準(zhǔn)校正的效果受觀測比例尺和隨機誤差的影響小，重力異常計算的均方根誤差為0.013 mGal、相對均方根誤差為0.433%，磁力異常計算的均方根誤差為3.336 nT、相對均方根誤差為0.175%，計算精度較高.同時，進(jìn)行了不同空白范圍的測試，即保持觀測面五的范圍不變，觀測面四范圍占總區(qū)域比例為9.46%、25.00%、37.86%、47.92%、71.59%，并進(jìn)行誤差統(tǒng)計(表5)，可以看出兩組數(shù)據(jù)重疊區(qū)不同時，多尺度位場數(shù)據(jù)基準(zhǔn)校正計算結(jié)果穩(wěn)定性好，重疊區(qū)覆蓋范圍越廣，校正計算精度越高.

1.3 尺度統(tǒng)一方法

尺度統(tǒng)一需考慮數(shù)據(jù)觀測精度或比例尺的不同，在重疊區(qū)需根據(jù)精度或比例尺加權(quán)計算；為滿足后期數(shù)據(jù)處理和解釋的數(shù)據(jù)需求以及高精度數(shù)據(jù)的完全保留，需將多尺度數(shù)據(jù)統(tǒng)一至最高精度的位場數(shù)據(jù)尺度，本文采用最小曲率方法進(jìn)行網(wǎng)格化(王萬銀和邱之云，2011).

重、磁位場數(shù)據(jù)是地下密度、磁性不均勻體的綜合反映，高精度數(shù)據(jù)包含區(qū)域場信息和更豐富的細(xì)節(jié)信息，為在數(shù)據(jù)融合能保留高精度數(shù)據(jù)，在數(shù)據(jù)重疊區(qū)精度加權(quán)的原則為：高精度數(shù)據(jù)權(quán)重為1，低精度權(quán)重為0；若精度相同，則權(quán)重相同(即取平均值).需要注意的是，在進(jìn)行位場延拓時，延拓的距離越大，計算的誤差相應(yīng)也越大，因此不同高度數(shù)據(jù)經(jīng)位場延拓?fù)Q算至同一個平面或曲面后可信度也有所差別，所以，多維單尺度數(shù)據(jù)重疊區(qū)也存在精度差異，需要加權(quán)處理.

圖11 單尺度原始重力、磁力數(shù)據(jù)等值線圖(a) 重力異常一； (b) 重力異常二； (c) 磁力異常一； (d) 磁力異常二.白色區(qū)域為數(shù)據(jù)空白區(qū).Fig.11 The contour map of original single-scale gravity and magnetic data(a) The first gravity anomaly; (b) The second gravity anomaly; (c) The first magnetic anomaly; (d) The second magnetic anomaly. The white areas indicate that there has no data in the areas.

圖12 單尺度重力、磁力異?；貧w分析計算結(jié)果及誤差等值線圖(a) 重力異?；貧w分析結(jié)果； (b) 磁力異?；貧w分析結(jié)果； (c) 重力異?；貧w分析計算誤差； (d) 磁力異常回歸分析計算誤差.Fig.12 The contour map of regression results of single-scale gravity and magnetic anomalies(a) The regression results of gravity anomalies; (b) The regression results of magnetic anomalies; (c) The error of regression results of gravity anomalies; (d) The error of regression results of magnetic anomalies.

圖13 多尺度觀測面、原始數(shù)據(jù)及回歸分析結(jié)果等值線圖(a) 觀測面四； (b) 觀測面四重力異常； (c) 觀測面五重力異常； (d) 重力異常校正結(jié)果； (e) 觀測面五； (f) 觀測面四磁力異常； (g) 觀測面五磁力異常； (h) 磁力異常校正結(jié)果.白色區(qū)域為數(shù)據(jù)空白區(qū).Fig.13 The contour map of multi-scale observation surface, raw data and regression result(a) The observation surface Four; (b) The gravity anomaly on observation surface Four; (c) The gravity anomaly on observation surface Five; (d) The regression results of gravity anomalies; (e) The observation surface Five; (f) The magnetic anomaly on observation surface Four; (g) The magnetic anomaly on observation surface Five; (h) The regression results of magnetic anomalies. The white areas indicate that there has no data in the areas.

表5 理論模型多尺度重、磁力異?；貧w分析校正誤差統(tǒng)計表Table 5 The error table of multi-scale gravity or magnetic anomaly regression of theoretical models

根據(jù)多尺度位場數(shù)據(jù)融合的方法原理，選取已經(jīng)完成維度統(tǒng)一和基準(zhǔn)差統(tǒng)一的單尺度數(shù)據(jù)(圖12)和多尺度數(shù)據(jù)(圖13d、h)，進(jìn)行單尺度和多尺度融合(圖14).由圖14可以看出，位場數(shù)據(jù)融合計算結(jié)果能夠保留高精度、大比例尺數(shù)據(jù)的細(xì)節(jié)信息，同時對研究區(qū)內(nèi)空白區(qū)進(jìn)行了有效補充，不同數(shù)據(jù)接邊位置光滑，沒有梯級帶和突跳.同時，為測試重疊區(qū)范圍對多尺度數(shù)據(jù)融合結(jié)果的影響，設(shè)計不同重疊區(qū)模型測試并統(tǒng)計了融合計算結(jié)果的誤差(表6)，重疊區(qū)覆蓋范圍不同時，融合計算誤差均很小，相對均方根誤差均小于1%，計算穩(wěn)定；且重疊區(qū)范圍越廣，融合精度越高.因此，本文提出的多尺度位場數(shù)據(jù)融合方法在數(shù)據(jù)重疊區(qū)范圍不同時均可得到較好的融合計算結(jié)果.

最終設(shè)計多維多尺度理論模型測試其位場數(shù)據(jù)融合計算的正確性和有效性.場源參數(shù)如表1，設(shè)計3個不同高度、不同尺度、數(shù)據(jù)基準(zhǔn)有差異的不同覆蓋范圍的觀測面，觀測面一(圖3b)，觀測面六(在觀測面一的基礎(chǔ)上z坐標(biāo)增大1000 m，z坐標(biāo)范圍為0～-1500 m，設(shè)置空白區(qū))，觀測面七(在觀測面一的基礎(chǔ)上z坐標(biāo)減小2000 m，z坐標(biāo)范圍為-3000～-4500 m；網(wǎng)格間距調(diào)整為200 m并設(shè)置空白區(qū))，以其理論重力異常(圖15b)和磁力異常(圖15c)進(jìn)行融合測試.

采用本文的多維多尺度重、磁位場數(shù)據(jù)融合方法.首先將多維數(shù)據(jù)以維度統(tǒng)一方法計算至同一個觀測面上變?yōu)閱尉S數(shù)據(jù)，接著分別統(tǒng)一單維多尺度數(shù)據(jù)的基準(zhǔn)差異和尺度差異，得到該觀測面上的融合結(jié)果，利用位場延拓得到了z=0 m、z=-2000 m、z=-4000 m和z=-6000 m四個平面的重、磁力異常(圖16).可以看到其融合結(jié)果與各個計算面的理論重、磁力異常均吻合較好，重、磁力異常等值線光滑、異常變化趨勢一致；隨著延拓高度的增加，異常形態(tài)變化趨于平緩，符合位場數(shù)據(jù)隨高度變化的規(guī)律；融合計算結(jié)果與理論異常的誤差統(tǒng)計表(表7)顯示，其計算至不同高度平面結(jié)果均較穩(wěn)定，計算精度較高，驗證了該方法在多維多尺度重磁位場數(shù)據(jù)融合中的有效性.

圖14 單尺度和多尺度重、磁位場數(shù)據(jù)融合結(jié)果等值線圖(a) 單尺度重力異常融合結(jié)果； (b) 單尺度磁力異常融合結(jié)果； (c) 多尺度重力異常融合結(jié)果； (d) 多尺度磁力異常融合結(jié)果.Fig.14 The contour map of fusion result of single-scale and multi-scale gravity and magnetic anomalies(a) The fusion result of single-scale gravity anomalies; (b) The fusion result of single-scale magnetic anomalies; (c) The fusion result of multi-scale gravity anomalies; (d) The fusion result of multi-scale magnetic anomalies.

表6 理論模型多尺度重、磁力異常融合計算誤差統(tǒng)計表Table 6 The error table of multi-scale gravity or magnetic anomaly fusion of theoretical models

圖15 多維多尺度原始觀測面起伏Fig.15 The 3D surface map of original multi-dimension multi-scale observation surface

圖16 多維多尺度重、磁異常及融合計算結(jié)果立體圖(a) 觀測面重力異常； (b) 觀測面磁力異常； (c) 重力異常融合結(jié)果； (d) 磁力異常融合結(jié)果.Fig.16 The 3D surface map of multi-dimension multi-scale anomalies and fusion results of gravity and magnetic anomalies(a) The gravity anomalies on the observation surface; (b) The magnetic anomalies on the observation surface; (c) The fusion results of gravity anomalies; (d) The fusion results of magnetic anomalies.

表7 理論模型多維多尺度重、磁力異常融合誤差統(tǒng)計表Table 7 The error table of multi-dimension multi-scale gravity or magnetic anomaly fusion of theoretical models

2 實際資料處理試驗

為檢驗本文多維多尺度重、磁位場數(shù)據(jù)融合方法實用性，選取在中國某地航空和地面觀測的不同比例尺重力異常數(shù)據(jù)進(jìn)行多維多尺度位場數(shù)據(jù)融合試驗.

2.1 單維多尺度重力數(shù)據(jù)融合試驗

選取中國某地1∶50000地面實測重力異常(數(shù)據(jù)體一，圖17b)和1∶200000地面實測重力異常數(shù)據(jù)(數(shù)據(jù)體二，圖17c)檢驗單維多尺度重、磁位場數(shù)據(jù)融合方法的實用效果.兩個不同比例尺的地面數(shù)據(jù)的觀測面均為地形面(圖17a)，地形高度起伏變化范圍是747～2623 m，兩次重力觀測得到了不同比例尺的數(shù)據(jù)體一(幅值-8～100 mGal)和數(shù)據(jù)體二(幅值-43～96 mGal)，網(wǎng)格化后X和Y方向的網(wǎng)格間距都為150 m，其異常趨勢和幅值變化在相同位置上大體一致，但1∶50000的數(shù)據(jù)反映更多細(xì)節(jié)信息，而1∶200000數(shù)據(jù)的變化趨勢相對光滑.

利用單維多尺度重、磁位場數(shù)據(jù)融合方法將兩組不同比例尺的重力數(shù)據(jù)進(jìn)行融合，以高精度大比例尺數(shù)據(jù)為基準(zhǔn)，將小比例尺數(shù)據(jù)進(jìn)行數(shù)據(jù)基準(zhǔn)差校正，相關(guān)系數(shù)為0.883，回歸模型的斜率為0.996，截距(即數(shù)據(jù)基準(zhǔn)差)為5.636 mGal；校正后的小比例尺數(shù)據(jù)與大比例尺數(shù)據(jù)進(jìn)行偏差統(tǒng)計，其偏差平均值為0.017 mGal，均方根偏差為7.546 mGal，相對均方根偏差6.968%，其偏差較小能夠保證校正后數(shù)據(jù)與高精度數(shù)據(jù)的一致性；之后進(jìn)行數(shù)據(jù)融合得到單維多尺度重力數(shù)據(jù)融合結(jié)果(圖17d)，融合結(jié)果重力異常的形態(tài)和幅值(33～100 mGal)均與原觀測面基本一致，保留了地面高精度的大比例尺數(shù)據(jù)(圖17b)所包含的細(xì)節(jié)信息，同時用較低精度的小比例尺數(shù)據(jù)有效補充了高精度大比例尺數(shù)據(jù)的空白區(qū)，融合后的重力異常也符合位場的特點，驗證了該融合方法在不同尺度實測重力異常的融合中的實用性.

圖17 地形起伏、觀測重力異常及融合重力異常等值線圖(a) 地形起伏； (b) 數(shù)據(jù)體一1∶50000地面實測重力異常； (c) 數(shù)據(jù)體二1∶200000地面實測重力異常； (d) 單維多尺度融合重力異常；黑色點代表觀測點.Fig.17 The contour map of topography, measured gravity anomaly and the fusion result(a) The topography; (b) The gravity anomaly measured on the ground at scale 1∶50000; (c) The gravity anomaly measured on the ground at scale 1∶200000; (d) The fusion result of single-dimension multi-scale gravity anomaly; The black points represent the location of the observation points.

2.2 多維多尺度重力數(shù)據(jù)融合試驗

為檢驗多維多尺度重、磁位場數(shù)據(jù)融合方法的實用性，在地面不同比例尺重力異常的基礎(chǔ)上增加1∶50000航空實測重力異常數(shù)據(jù)，航空觀測面GPS高度為1638～2705 m(圖18a)，測線間隔為500 m，測點間隔為20 m，重力異常幅值為-43～96 mGal(圖18b)，X和Y方向的網(wǎng)格化尺寸為150 m，航空重力異常與地面重力異常變化趨勢大體一致，但總體幅值較低；兩者比較看出，地面數(shù)據(jù)顯示更多細(xì)節(jié)信息，航空數(shù)據(jù)觀測比例尺雖與地面1∶50000重力數(shù)據(jù)一致，但其觀測精度低于地面1∶50000數(shù)據(jù)，航空數(shù)據(jù)主要在覆蓋范圍和工作效率上更具有優(yōu)勢.

采用本文多維多尺度重、磁位場數(shù)據(jù)融合方法對數(shù)據(jù)體一(1∶50000地面重力異常，圖17b)、數(shù)據(jù)體二(1∶200000地面重力異常，圖17c)和數(shù)據(jù)體三(1∶50000航空重力異常，圖18b)三組航空和地面實測的不同精度重力異常數(shù)據(jù)進(jìn)行融合試驗.

數(shù)據(jù)體一和數(shù)據(jù)體三融合過程中，航空數(shù)據(jù)采用位場延拓計算至地面(維度統(tǒng)一)后進(jìn)行基準(zhǔn)差統(tǒng)一，與地面數(shù)據(jù)對比，其平均偏差小于10-5mGal，均方根偏差為7.379 mGal，相對均方根偏差為6.813%.可以看出，數(shù)據(jù)體三計算至地面后的精度略高于數(shù)據(jù)體二(相對均方根偏差為6.968%)，但精度相差不大；數(shù)據(jù)融合重力異常如圖19a所示，融合重力異常符合位場數(shù)據(jù)特征，形態(tài)和幅值均與原觀測面基本一致，幅值范圍為-44～124 mGal，保留了地面高精度數(shù)據(jù)(圖17b)包含的細(xì)節(jié)信息，同時用效率高、覆蓋范圍廣的航空重力數(shù)據(jù)有效地補充了地面重力數(shù)據(jù)的空白區(qū)，兩者數(shù)據(jù)相接的位置趨勢一致且光滑無梯級帶.

根據(jù)上述兩組數(shù)據(jù)融合的精度比較結(jié)果，該地區(qū)三組多維多尺度數(shù)據(jù)融合時，首先將計算至地面后精度相差不大的航空1∶50000和地面1∶200000重力異常數(shù)據(jù)融合，加權(quán)平均的權(quán)重為0.5，得到初步融合結(jié)果(圖19b)，再將其與精度最高的地面1∶50000重力數(shù)據(jù)進(jìn)行多尺度融合，即為此三組多維多尺度數(shù)據(jù)融合重力異常(圖19c).

數(shù)據(jù)體二和三的融合結(jié)果(圖19c)已經(jīng)結(jié)合了大范圍的航空重力數(shù)據(jù)與小比例尺小范圍地面數(shù)據(jù)的各自優(yōu)勢，經(jīng)校正后異常形態(tài)趨勢一致，最終融合結(jié)果(圖19c)相對于初步融合結(jié)果又包含更多地面高精度細(xì)節(jié)信息；統(tǒng)計初步融合結(jié)果校正至地面高精度數(shù)據(jù)后的偏差，平均值小于10-5mGal，均方根偏差為6.556 mGal，相對均方根偏差為6.054%，由此看出三組數(shù)據(jù)整體融合的效果優(yōu)于任意兩組數(shù)據(jù)，同樣驗證了該思路和方法在多維多尺度位場數(shù)據(jù)融合中的有效性.

圖18 航空觀測GPS高度和實測重力異常等值線圖(a) 觀測面GPS高度； (b) 1∶50000實測航空重力異常.Fig.18 The contour map of observed surface GPS altitude and gravity anomaly measured by aviation(a) The GPS altitude of observation surface; (b) Measured gravity anomaly at scale 1∶50000 by aviation.

圖19 多維多尺度融合重力異常等值線圖(a) 航空1∶50000與地面1∶50000重力異常融合； (b) 航空1∶50000與地面1∶200000重力異常融合； (c) 航空1∶50000與地面1∶200000、地面1∶50000重力異常融合； (d) 航空1∶50000與地面1∶200000、地面1∶50000重力異常當(dāng)作單維直接處理結(jié)果.Fig.19 The contour map of fusion results of multi-dimension multi-scale gravity anomaly(a) The fusion result of aviation 1∶50000 and ground 1∶50000 gravity anomaly; (b) The fusion result of aviation 1∶50000 and ground 1∶200000 gravity anomaly; (c) The fusion result of aviation 1∶50000, ground 1∶200000 and ground 1∶50000 gravity anomaly; (d) The processing result of aviation 1∶50000, ground 1∶200000 and ground 1∶50000 gravity anomaly regarded as single-dimension data.

另外，為了進(jìn)一步驗證和對比多維數(shù)據(jù)融合中維度統(tǒng)一的必要性和有效性，本文測試將多維數(shù)據(jù)當(dāng)作單維數(shù)據(jù)不進(jìn)行位場延拓直接“融合”(圖19d)，其結(jié)果顯示在不同數(shù)據(jù)的接邊處存在明顯梯級帶，且異常值的變化趨勢無法達(dá)到有效一致，不符合位場變化規(guī)律，由此也驗證不考慮觀測位置高度變化而將多維數(shù)據(jù)當(dāng)作單維數(shù)據(jù)直接“融合”的做法是錯誤的.本文多維多尺度重、磁位場數(shù)據(jù)融合方法原理及計算技術(shù)對于不同觀測比例尺的地面和航空數(shù)據(jù)進(jìn)行了有效融合，驗證了該方法在實測地面和航空數(shù)據(jù)融合中的實用性.

3 結(jié)論與建議

本文研究了多維多尺度重、磁位場數(shù)據(jù)融合問題，采用改進(jìn)的積分-迭代位場延拓方法實現(xiàn)維度統(tǒng)一，采用相關(guān)分析和回歸分析實現(xiàn)數(shù)據(jù)基準(zhǔn)統(tǒng)一，采用數(shù)據(jù)精度加權(quán)平均的方法實現(xiàn)尺度統(tǒng)一，解決了多維多尺度重、磁位場數(shù)據(jù)融合中維度、尺度和數(shù)據(jù)基準(zhǔn)差異的問題，提出了一套適用于多維多尺度重、磁位場數(shù)據(jù)融合的方法.

理論模型測試和實測資料試驗顯示，該方法適用于航空、地面或海洋觀測的重、磁位場數(shù)據(jù)融合，融合重、磁力異常能夠作為后續(xù)位場數(shù)據(jù)處理、反演和解釋的基礎(chǔ)數(shù)據(jù)，具有良好的實用價值.本文提出的重、磁位場數(shù)據(jù)融合的思路也可推廣至重、磁力異常各分量或張量數(shù)據(jù)的融合處理中.

致謝感謝論文評審專家提出的修改意見和編輯對論文的編輯、加工.