高維強，史朝洋，張利明，馮旭亮

(1.陜西省礦產(chǎn)地質(zhì)調(diào)查中心，陜西 西安 710068；2.西安石油大學 地球科學與工程學院，陜西 西安 710065)

0 引言

實測重力數(shù)據(jù)經(jīng)過高度校正、緯度校正、地形校正和中間層校正，即可得到布格重力異常，其為重力資料解釋的基礎。然而，對于地形變化較大的山區(qū)，布格重力異常往往與地形高程呈鏡像或同像變化，即山形異常[1-2]，給地質(zhì)構(gòu)造解釋推斷帶來極大的困難。

引起山形異常的主要原因為地形和中間層改正時中間層密度不準確[3]，可利用實測重力數(shù)據(jù)與高程散點圖擬合的中間層校正系數(shù)C=2πGσ反算中間層密度[2,4]。但地形起伏較大的山區(qū)，其地形物質(zhì)規(guī)?？s小，所產(chǎn)生的重力效應小于布格板重力效應，據(jù)此反算的中間層密度偏小[1]。Thorarinsson and Magnusson[5]將分形分析引入重力勘探領域，在冰島西南部取得了良好的地質(zhì)效果，但該技術(shù)必須對復雜地表按同一密度層分區(qū)，區(qū)內(nèi)巖石密度越單一，該方法的效果越好[6]，因此實際應用中較難做到[7]。

相關分析法是選取最佳中間層密度的傳統(tǒng)方法，一般在重力剖面上選用一系列中間層密度值進行校正，選用與地形無關或最小相關的布格重力異常值對應的密度值作為最佳中間層密度[8-9]。在同一測區(qū)，如果表層密度較為均勻，則圖解剖面法是有效的[10]，否則圖解法所求密度只是剖面位置處的合適密度，不能代表整個測區(qū)密度[10-11]，可對全區(qū)地改后的布格重力異常與地形高程進行相關分析，選取相關性最小的地層密度作為最佳地改密度[11]。當?shù)乇砻芏葯M向變化較大時，可分區(qū)選用最佳密度值并拼接形成全區(qū)變化的中間層密度[12]，也可在常密度改正結(jié)果的基礎上進行補償校正[13-15]。然而，當研究區(qū)面積較大且需要分區(qū)計算時，計算量較大，并且只能得到與最佳中間層密度最為接近的密度而并非最“真實”的密度。

蔡尚中[2]建立了自由空間重力異常與高程之間的函數(shù)關系，之后根據(jù)每個測點實際高程計算出回歸異常并從實測自由空間重力異常中消除，得到回歸剩余異常，有利于定性解釋。當研究區(qū)較大時，可采用滑動窗口回歸分析方法，并可利用變差函數(shù)給出窗口大小選取的量化參考值[16]。然而，回歸分析方法僅考慮了研究區(qū)范圍之內(nèi)地形和重力異常的關系，得到的回歸重力值僅為與地形起伏呈整體區(qū)域性變化的重力異常，并不能消除局部異常(特別是研究區(qū)以外地形的重力影響)，此外，該方法消除的區(qū)域異常中可能包含了深部地質(zhì)信息引起的有用的重力異常。

若有研究區(qū)地表地質(zhì)露頭實測的巖石密度，可采用插值和圓滑方法構(gòu)建地表密度分布模型，并建立浮動基準面進行地形和中間層校正[7]。然而當研究區(qū)面積較大時，有限的露頭采樣和鉆井測量資料并不能全面反映測區(qū)地表巖石密度的變化情況。牛源源等[17]針對這一問題，從一維自由空氣重力異常數(shù)據(jù)出發(fā)，采用貝葉斯方法估算近地表巖石密度，其與地表地質(zhì)及區(qū)域地層密度資料較為吻合，但目前僅用于剖面重力異常分析，并未見到面積性重力異常的相關應用。

重力數(shù)據(jù)同時包含了構(gòu)造和密度信息，已有研究表明利用重力數(shù)據(jù)可以較準確地估算近地表巖石密度[17]或密度界面上下的密度差[18]。Florio[19-20]進行沉積盆地基底深度重力反演時，根據(jù)回歸分析獲得盆地沉積層與基底的密度差，經(jīng)過幾次迭代便可獲得滿意的結(jié)果。受這一思路的啟發(fā)，本文在處理唐昭陵所在的九嵕山地區(qū)的高精度重力資料時，利用逐次回歸技術(shù)分析最佳地改密度，將其應用于重力資料外部校正之中，取得了較好的應用效果。

1 研究區(qū)概況

九嵕山位于鄂爾多斯盆地南緣，西鄰祁連—河西走廊(六盤山)—賀蘭山構(gòu)造帶，南隔渭河地塹與秦嶺造山帶相依，地處穩(wěn)定的鄂爾多斯地塊與活動的祁連—秦嶺造山帶之間(圖1)，具有長期復雜的地質(zhì)構(gòu)造背景。研究區(qū)一帶位于鄂爾多斯南緣翹起帶，為一強烈的近東西向褶皺帶，總體構(gòu)造形態(tài)為一復背斜。次級背、向斜東西成帶狀展布，褶皺緊密，局部向南倒轉(zhuǎn)，兩翼斷層多見。唐王陵向斜槽部為上奧陶統(tǒng)唐王陵組，南翼向北傾，350°-10°∠20°-30°，北翼向南傾，傾角較陡，一般產(chǎn)狀為170°-190°∠50°-60°。由于北翼北斷層破壞，地層出露不全，總體為一不對稱向斜。

圖1 九嵕山及鄰區(qū)無人機激光雷達測量的地形

研究區(qū)位于唐王陵向斜的南翼，主要出露奧陶系唐王陵組，按巖性可分為兩段：下段(O3t1)為雜色(黃褐、灰綠、紫紅、深灰)含礫頁巖，底部夾厚層狀巖屑砂巖，并發(fā)育燧石條帶白云巖大漂礫；上段(O3t2)為灰色、淺灰色巨厚層狀復成分礫巖、角礫巖，偶夾礫屑砂巖，上下兩段間為整合接觸。

2 地球物理特征

研究區(qū)共進行了1 722個點的重力測量，其中線距為10 m，點距為5 m。對原始重力測量值進行高度校正和正常場校正，得到自由空間重力異常如圖2所示，其與地形(圖3)呈明顯的同相變化關系，在不到200 m的海拔高度變化范圍內(nèi)，重力異常的局部變化接近10×10-5m/s2。

圖2 研究區(qū)自由空間重力異常

圖3 研究區(qū)地形

研究區(qū)大面積出露奧陶系唐王陵組，第四系非常薄，且僅分布于局部地區(qū)，最大厚度不超過1 m。未獲得研究區(qū)地形改正密度，在研究區(qū)不同位置實測了156塊標本的密度(位置如圖1中藍色五角星所示)，統(tǒng)計結(jié)果見表1，對于標本數(shù)量較少的巖性，僅給出平均密度。研究區(qū)各類巖石密度整體較大，大多數(shù)巖石密度為2.7×103kg/m3左右，砂礫巖和雜色礫巖的密度超過了2.77×103kg/m3。可見，盡管研究區(qū)主要出露奧陶系唐王陵組，但由于地層巖性變化較大，使得地層密度也并非統(tǒng)一的密度，且具有較為明顯的變化。

表1 研究區(qū)各類巖石密度統(tǒng)計

3 常規(guī)地形改正結(jié)果

嚴良俊等[7]根據(jù)實測的地表密度資料建立橫向的密度變化并用于地形校正，取得較好的地質(zhì)效果。根據(jù)此方法，本文對除唐瓦和唐磚之外實測的153個點上的巖石密度值進行網(wǎng)格化以建立表層變密度模型用于地形校正。為避免密度不連續(xù)引起的布格重力異常的畸變，對原始密度網(wǎng)格數(shù)據(jù)進行光滑，并根據(jù)地表地質(zhì)資料對實測密度外圍的巖石密度進行推斷補充及擴邊，最終形成與圖1所示的地形范圍一致的地表橫向變密度數(shù)據(jù)，作為地形改正的依據(jù)。

自由空間重力異常經(jīng)地形改正和中間層改正之后可得到布格重力異常。傳統(tǒng)的做法是先將任一測點周圍的地形“削平補齊”，然后利用布格板公式計算各測點與基點之間的無限大物質(zhì)層的重力影響，即先做地形改正，再做中間層改正。在計算測點的地形改正值時，先計算測點附近4個節(jié)點的地形改正值，然后將4個節(jié)點的地形改正值內(nèi)插到測點位置上來作為測點的地形改正值[21]。本文在進行地形改正時，根據(jù)實測的地形起伏直接計算地形在研究區(qū)各重力測點引起的重力變化，并從自由空間重力值中消除，得到布格重力異常。這一過程相當于將傳統(tǒng)的地形改正和中間層改正合并為一項進行，簡化了計算過程，并且避免了地形改正值內(nèi)插的誤差。因此，如不特別說明，本文中的地形改正均指的是傳統(tǒng)的地形改正和中間層改正之和。

如果重力測量僅用于研究范圍較小的局部異常，地形校正范圍可以采用最大半徑為20 km進行校正，把遠區(qū)地形校正當做區(qū)域異常將其消除即可[22]。九嵕山研究區(qū)重力測量的目的是推斷昭陵墓道位置及墓室結(jié)構(gòu)，其為局部淺地表目標體。因此，本文僅計算了圖1所示范圍內(nèi)地形的影響而忽略了該區(qū)外圍地形的影響。計算地形影響時選擇最低測點所在的位置為基點位置，并考慮測點之外地形的整體變化趨勢，最終選擇1 000 m為基點所在平面，正演計算實測地形表面與該平面之間的物質(zhì)在各測點所引起的重力異常，如圖4a所示，從自由空間重力異常中消除地形影響，得到布格重力異常如圖4b所示。

圖4b所示的布格重力異常與地形呈鏡像關系，即呈明顯的山形異常，二者之間相關系數(shù)為-0.989，可見布格重力異常受地形的影響非常嚴重。從圖4a來看，地形影響值變化形態(tài)與自由空間重力異常形態(tài)非常相似，二者的相關系數(shù)為0.996，但前者相對變化明顯大于后者，其原因為利用實測數(shù)據(jù)建立的地形改正密度偏大。通過野外踏勘也發(fā)現(xiàn)，研究區(qū)地表巖石風化比較嚴重，導致巖石的實際等效密度遠小于實測密度。圖5為不同測點處地表巖石露頭照片，可以看出巖石裂隙非常發(fā)育且不同位置裂隙規(guī)模差異較大，并且裂隙在縱向上延伸范圍也比較大。此外，部分裂隙中充填了黃土，部分裂隙中無充填物。研究區(qū)各類巖石裂隙規(guī)模的變化及充填物的變化導致最佳地形改正密度無法利用實測巖石密度及裂隙的發(fā)育程度去估計，只能采用間接的方法，即根據(jù)實際地形與重力異常的關系進行估計。

圖4 研究區(qū)由實測密度計算的地形影響值及相應的布格重力異常

圖5 研究區(qū)部分測點處地表露頭俯視(a～e)和遠景(f)照片

4 逐次回歸地形改正密度選取方法及應用

在地形起伏較大的山區(qū)利用實測重力數(shù)據(jù)與高程散點圖擬合的中間層校正系數(shù)，進而根據(jù)布格板公式反算的地形改正密度偏小[1]，可對該密度進行迭代調(diào)整而得到最佳地改密度。利用重力資料反演密度界面深度時，常利用布格板重力公式構(gòu)建密度界面深度修改量的迭代計算公式，并逐次迭代消除剩余重力異常值則可反演得到密度界面深度[23-24]。受迭代法的思想啟發(fā)，本文在選取最佳地形改正密度時采用以下步驟：

1)根據(jù)自由空間重力異常與高程的散點圖，利用公式進行擬合，并計算第一次地形改正密度σ(1)=a/(1.6πG)，其中G為萬有引力常量。需要說明的是，由于三維條件下，實際地形引起的重力變化小于相同密度時布格板的重力異常，因此本文將布格板公式中的系數(shù)2πG改為1.6πG以加速迭代收斂過程。

2)以σ(1)為地形改正密度，正演計算地形在測點處的重力影響值gt(1)，并根據(jù)gB(1)=gf-gt(1)計算布格重力異常。根據(jù)自由空間重力異常與第一次計算的地形在測點處的重力影響值的散點圖，利用公式gf=c·gt(1)+d進行擬合，并根據(jù)系數(shù)c的大小判斷：當c明顯大于1時，說明計算的地形影響值偏小，地形改正密度σ(1)偏小，繼續(xù)下一步驟；當c明顯小于1時，說明計算的地形影響值偏大，地形改正密度σ(1)偏大，繼續(xù)下一步驟；當c≈1時，說明地改密度σ(1)取值合適，此時停止迭代，將布格重力異常作為最終結(jié)果。

3)根據(jù)第一次計算得到的布格重力異常與高程的散點圖，利用公式gB(1)=e·h+f進行擬合，計算第一次地形改正密度修正值δ(1)=e/(1.6πG)。

4)根據(jù)σ(2)=σ(1)+δσ(1)計算第二次地形改正密度，并令σ(1)=σ(2)轉(zhuǎn)回步驟2)。

根據(jù)以上步驟，本文計算了九嵕山地區(qū)布格重力異常。研究區(qū)自由空間重力異常與測點高程的散點圖及擬合公式如圖6所示。由擬合一次項系數(shù)得到第一次地改密度值σ(1)=1.489×103kg/m3，據(jù)此計算得到地形重力影響值及相應的布格重力異常如圖7a、b所示。此時自由空間重力異常與圖7a所示的地形影響的散點圖(圖8)中，線性擬合的一次項系數(shù)為1.4843，其明顯大于1，說明計算的地形影響值偏小。根據(jù)布格重力異常與高程的散點圖(圖9)，由擬合一次項系數(shù)得到第一次地形改正的修正值δ(1)=0.463×103kg/m3，則第二次的地形改正密度σ(2)=1.952×103kg/m3。

圖6 研究區(qū)自由空間重力異常與高程散點

圖7 研究區(qū)由第一次回歸密度計算的地形影響值及相應的布格重力異常

圖8 研究區(qū)自由空間重力異常與第一次計算的地形影響值散點

圖9 研究區(qū)第一次計算的布格重力異常與高程散點

利用此方法，經(jīng)過5次迭代，得到最佳地形改正密度σ(5)=2.154×103kg/m3。地改密度隨迭代次數(shù)的變化如圖10a所示，自由空間重力異常與計算的地形在測點處的重力影響值的相關系數(shù)c隨迭代次數(shù)的變化如圖10b所示。迭代初始，密度隨迭代次數(shù)變化較大，隨著迭代次數(shù)增加，密度值變化不明顯，相關系數(shù)c也逐漸趨于1，證明此事已經(jīng)得到了最佳密度，再增加迭代次數(shù)意義不大。

圖10 密度及相關系數(shù)隨迭代次數(shù)的變化

假設研究區(qū)地表巖石平均密度為2.7×103kg/m3，則該最佳密度值近似于未風化巖石密度的80%，從地表巖石露頭的裂隙發(fā)育來看，估算裂隙度10%～30%，考慮裂隙中存在部分黃土充填，因此本文估算的最佳地形改正密度應該是合理的。根據(jù)該密度計算的布格重力異常如圖11a所示，其與地形的相關系數(shù)僅為0.0329，可見布格重力異常與地形完全無關。

圖11 研究區(qū)布格重力異常

研究區(qū)有4個大小不等的石硐(位置如圖11中黑色方框所示)，尺寸約為4 m(寬)×3 m(高)×5.5 m(深)，石硐頂距地表約8～10 m。若按照上述估算的中間層平均密度，則石硐與圍巖的密度差為-2.154×103kg/m3，正演計算4個石硐引起的重力異常，其為局部重力異常低，最大幅值可達-0.046×10-5m/s2。利用位場分離方法求取剩余布格重力異常如圖11b所示，4個石硐處表現(xiàn)為明顯的局部重力異常低，相同參數(shù)下對圖4b所示的布格重力異常求取的剩余異常則無法反映石硐的重力異常特征，且剩余重力異常走向與地形呈明顯相關關系?？梢?，本文利用逐次回歸方法得到的中間層密度較為合理，能較好地消除與地形無關的虛假異常。

5 討論

本文根據(jù)自由空間重力異常與地形的相關關系，利用逐次回歸方法確定了山區(qū)地形改正的密度，九嵕山地區(qū)的實測重力異常校正結(jié)果說明該方法的正確性。對于山區(qū)重力數(shù)據(jù)的地形改正問題，最常用的方法是選取一系列地形改正密度值進行試算，將與地形相關性最小的布格重力異常確定為最佳重力異常，對應的密度為最佳地改密度[11]。使用該方法時，密度間隔取值太小，計算量過大；反之密度間隔取值太大時，可能不能獲取最佳地形改正密度。常用的選取方法是以0.05×103kg/m3為步長進行計算，即便如此，若在(2.0～2.7)×103kg/m3的范圍內(nèi)進行試算，也需15次地形改正計算，工作量非常大。在給定的區(qū)間內(nèi)采用二分法選取密度值進行地形改正計算[10]能在一定程度上減小工作量，但若想得到較準確的最佳地改密度，仍然需要多次計算才行。而本文提出的方法只需有限的幾次計算即可得到最佳地改密度，工作量較小。

九嵕山實測重力資料校正時，全區(qū)采用的統(tǒng)一的密度。其原因在于地表地質(zhì)調(diào)查結(jié)果表明研究區(qū)出露的地層較為單一，可以用統(tǒng)一的密度去等效實際地層的密度變化。此外，自由空間重力異常與測點高程的回歸分析結(jié)果可以較好地用線性函數(shù)擬合，也反映出該區(qū)地形改正密度無明顯的橫向變化。實際工作中，若測區(qū)面積較大，可利用自由空間重力異常與測點高程的散點圖進行判斷，若線性函數(shù)擬合誤差較大，則可采用分區(qū)處理的方式，最終再將校正結(jié)果拼接即可。

本文將布格重力異常與地形起伏變化無關作為判斷布格重力異常校正是否合理的標準，事實上其為山區(qū)地形改正時的普遍原則。然而，實際工作中，必須針對勘探的具體目標進行具體分析。例如區(qū)域重力調(diào)查中，通常會采用統(tǒng)一的中間層密度(一般取為地殼平均密度2.67×103kg/m3)進行校正以方便數(shù)據(jù)拼接及對比，此時布格重力異常通常與地形呈現(xiàn)鏡像關系，區(qū)域性的重力變化反映了地殼厚度的變化特征。在局部重力測量時，布格重力異常也可能表現(xiàn)為與地形同像的特征，例如盆山結(jié)合部位，造山帶通常由變質(zhì)巖、火成巖等組成，而盆地內(nèi)部沉積層較厚，會使得重力異常呈現(xiàn)山區(qū)高、盆地低的特征，這種同像的重力異常是符合地質(zhì)意義的，若研究盆山構(gòu)造特征，則不能簡單的以布格重力異常與地形無關而判斷。而本文中的九嵕山地區(qū)，勘探目標為昭陵墓道口及墓室，其為近地表局部異常，因此整個山體引起的重力變化為非目標異常，布格重力異常的計算過程也相當于消除了區(qū)域背景。因此，實際工作中應將地形校正的過程與布格重力異常的后續(xù)處理結(jié)合起來，以提取勘探對象引起的重力異常為最終目的，若在地形校正時采用地殼平均密度等常密度值進行校正，則在重力資料處理時需要采用一定的措施消除或減弱與勘探目標無關的異常；相反，若在地形校正時采用了分區(qū)校正、變密度校正等措施已在一定程度上減弱了非勘探目標的重力影響，則在重力資料處理中可針對性的進行處理即可。

6 結(jié)論

針對山區(qū)重力測量時地形改正密度選取不準確引起的山形異常問題，本文以九嵕山重力異常為例，基于回歸分析方法選取了最佳的地形改正密度，改正后的布格重力異常明顯不受地形影響，證明了這一方法的正確性。本文使用回歸分析方法時，采取了迭代計算的措施，與常規(guī)的采用一系列密度值進行試算的方法相比，本文方法可以快速獲得最佳地形改正密度。從九嵕山地區(qū)的重力測量工作來看，現(xiàn)有技術(shù)條件高精度地形數(shù)據(jù)獲取已不是難點，以實際地質(zhì)構(gòu)造為約束的改正方法及匹配的數(shù)據(jù)處理方法，是改進重力改正效果及提高解釋精度的關鍵之一。