一種用于重力場分離的約束插值切割方法

孟慶奎 舒 晴 高 維* 周堅鑫 張文志 李 勇

(①自然資源部航空地球物理與遙感地質重點實驗室，北京 100083； ②中國自然資源航空物探遙感中心，北京100083； ③中國地質大學(北京)地球物理與信息技術學院，北京 100083； ④中國地質科學院地球物理地球化學勘查研究所，河北廊坊 065000； ⑤自然資源部地球物理電磁法探測技術重點實驗室，河北廊坊 065000)

0 引言

重力場分離是重力數據反演、地質解釋與推斷的前提[1-3]，是業內公認的重點研究課題和技術難題。早在1952年，Dobrin等[4]就提出重力解釋中兩個最重要的問題：一個是區域場與局部場的分離； 另一個是基于分離得到的異常獲取目標地質體的空間位置和密度信息。程方道等[5]指出，與磁異常相比，重力異常普遍存在強大的區域背景場，往往會造成局部異常發生畸變，甚至難于辨識，要對重力資料進行有效的地質解釋，首先必須分離重力區域場與局部場； 管志寧等[6]強調，正確劃分不同規模的區域場與局部場是面積性測量資料解釋能否獲得良好地質效果的一個重要因素； 曾華霖[7]認為，在密度界面反演中，只要平均深度取值合理，界面之間的密度差異比較接近實際情況，無論采用何種反演方法，都能得到達到一定的精度及相似的結果，但前提是用于反演的異常必須單純由目標體引起； 徐世浙等[8]認為重力場三維反演是個難度極大的課題，提出異于常規解決方案的開拓性思路，主要包括三個方面：對平面上的重力場進行不同深度的分離； 用大深度位場向下延拓的迭代法將各個切割層在地面的重力場下延至相應深度； 將頂面的重力場反演為視密度。現今，位場大深度的向下延拓技術和層源異常的反演技術已經得到了很好發展，位場分離方法是目前影響反演方法效果的關鍵技術，繼續加強對異常分離方法的研究具有重要意義。

重力場分離方法大致可分為空間域方法和頻率域(波數域)方法。頻率域方法主要包括維納濾波法[9]、匹配濾波法[10-12]、小波多尺度分析法[13-17]、優選濾波法[18-19]和優化濾波法[20]等，這些方法各有特點和優勢，但同時又有不同的局限性。維納濾波法是維納濾波器在信號與干擾不相關條件下的一個特例，適用于白噪聲的去除，不適合于場的分離，目前相關研究與應用不多； 匹配濾波法對于埋深和尺寸相差較大的場源分離效果較好，但是若異常的功率譜不具明顯的線性特征時，如何更加有效地確定深、淺源場的頻段分布，是該方法在實際資料處理過程中的一個難點； 小波多尺度分析法是一種新興的位場分離方法，在區域性深部地殼結構研究方面取得了一些成果，但是存在小波母函數選擇和階次選擇較難等問題； 優選濾波法實現了指定頻段或尺度的局部異常的分離，但實際應用中延拓高度往往很難正確選擇； 優化濾波法異常分離不受延拓高度的困擾，但其前提是目標層與其他深度層場源信息互不相關。實踐表明，空間域方法有一定的技術優勢。管志寧等[6]基于多年的重磁數據處理經驗認為，由于受測區范圍有限、數據離散化，以及不同埋深、不同位置場源體的重磁異常頻譜特征部分重疊等諸多因素的影響，在波數域變換過程中會產生明顯的誤差和畸變，只能用作初步的定性估計； Spector[21]也認為，波數域方法只能對深源(區域場)和淺源(局部場)異常作定性的劃分，不能用于定量研究，無法取代空間域中異常曲線擬合類方法。因而，本文以空間域方法為研究對象。

空間域方法主要包括趨勢面擬合法[22-24]、經驗模態分解法[25]、最小曲率法[26-28]、曲率濾波法[29-30]、有限單元法[31-37]、卡爾曼濾波法[38-39]、矩陣底秩分解法[40-41]和插值切割法[42-53]等。目前，空間域方法中較為前沿且應用較廣的重力場分離方法是插值切割法，該方法基于多次切割法發展而來。程方道等[5]首次提出多次切割法，其基本原理是以四點圓周平均值為切割算子，通過連續切割得到切割區域場，并將其從總異常中減去，便得到切割局部場。為了壓制多次切割法中的切割發散效應，文百紅等[42]基于新型算子首次提出了插值切割法； 之后，秦葆瑚[43]、段本春等[44]、徐世浙等[45]、邢怡等[46]、劉東甲[47]、葛粲等[48-49]圍繞優化切割算子、探尋切割半徑與異常埋深的關系展開了系統的研究工作，形成了以插值切割法為核心技術的方法系列，汪炳柱等[50]、楊金玉等[51]、萬勝榮等[52]和趙文舉等[53]將該方法廣泛應用于油氣勘探、礦產勘查、基礎地學研究等諸多領域。

以上提及的插值切割法及其改進算法很少涉及重力異常產生的物理屬性，多針對異常曲線本身進行插值切割計算。如果忽視異常產生的地質屬性，往往會造成分析問題不夠全面，甚至產生一定的計算誤差。故本文在分析插值切割法基本原理基礎上，證明了插值切割法的切割算子僅考慮了局部異常與區域異常彎曲方向相同的疊合情況，不能完全體現實際復雜的地下地質情況。因此，本文結合重力局部異常與區域異常的四種疊合關系，提出了一種基于異常約束條件的插值切割法，并將其命名為約束插值切割法。通過典型理論模型對比試算和實際測量數據應用，驗證了約束插值切割法的準確性和優越性，可為大區域的基礎地學研究提供技術支撐，在具備局部高精度重力資料的情況下，也可為油氣遠景區的預測提供一定技術參考。

1 插值切割法

插值切割法是一種簡單、有效的空間域非線性濾波方法，以當前計算點場值與周圍多點平均值的插值運算為切割算子，通過連續切割得到切割區域場，并將其從位場異常中減去，便得到切割局部場[42]。根據插值切割算子的不同，形成了基于四點圓周插值切割算子[42]、八點窗口插值切割算子[46]和動態改進型插值切割算子[49]的插值切割方法。

1.1 計算思路

插值切割法雖衍生出不同的切割算子，但其計算原理是一致的。下面以經典的四點圓周插值切割算子為例介紹插值切割法的計算思路。

設重力異常Z(x,y)由區域場R(x,y)和局部場L(x,y)組成，即

Z(x,y)=R(x,y)+L(x,y)

(1)

令M(x,y)表示與點(x,y)距離為r的4個點的重力異常平均值，即

Z(x,y+r)+Z(x,y-r)]

(2)

式中r也稱為切割半徑。區域場R(x,y)是Z(x,y)與M(x,y)的加權平均，即

R(x,y)=aM(x,y)+bZ(x,y)

(3)

式中a和b是加權系數，且a+b=1。令

(4)

(5)

ΔHx=Z(x+r,y)+Z(x-r,y)

(6)

ΔHy=Z(x,y+r)+Z(x,y-r)

(7)

(8)

(9)

e=c+d(0≤e≤2)

(10)

(11)

用上述方法得到的區域場稱為第1次切割的區域場，用R1(x,y)表示； 對R1(x,y)重復使用以上方法，得到第2次切割的區域場R2(x,y)； 依次迭代，當

(12)

有

(13)

上式得到的區域場即為切割半徑的區域場，將其與重力異常Z(x,y)相減，得到局部場

L(x,y)=Z(x,y)-R(x,y)

(14)

一般地，切割半徑r越大，得到的局部場所反映的密度體的深度和規模就越大。根據前人大量模擬實驗和實測資料[42-53]分析，以切割半徑r得到的局部場可近似代表深度r以上地層的重力效應，但尚未得出切割半徑r與密度體埋深的確切計算公式。由于切割算子中的加權系數引入了半二階差分量，它們與負二階水平導數即曲率成正比，因此插值切割對不同非線性異常有不同的作用，可提高對不同特征異常場的分離效果。

1.2 誤差傳遞分析

插值切割法采用當前計算點值與四點圓周平均值之間的插值關系作為切割算子G，將切割算子G作用于重力場Z(x,y)，即可得到區域場

(15)

設原始觀測數據各點的均方誤差δZ(x,y)=ε0，根據誤差傳遞公式[54]得到區域場的誤差傳遞關系為

[δM(x,y)]2

(16)

那么，第1次切割后區域場的誤差傳遞系數k1有

(17)

第i次切割后區域場的誤差傳遞系數ki有

(18)

由上式可見，多次切割后區域場的傳遞誤差不會放大。

2 改進的插值切割法

分析插值切割法基本理論公式，該方法主要對異常曲線本身進行插值切割計算，而對重力異常產生的地質屬性考慮不足。忽略產生異常的地質屬性得到的計算結果難免產生一定的誤差。下面首先在理論層面分析插值切割法的不足之處，再從實際地質屬性出發，結合重力場局部異常與區域異常的四種疊合關系，對插值切割法加以改進，提出一種基于異常約束條件的插值切割法，即約束插值切割法。

2.1 插值切割法不足的理論分析

根據插值切割法的基本理論可知，切割區域場計算公式(式(15))中，0≤e≤2，不難證明R(x,y)介于Z(x,y)與M(x,y)之間，這就意味著經插值切割法計算得到的區域場值不能超出總場異常值與周圍點平均值的取值區間。這種情形僅代表局部異常與區域異常彎曲方向相同的疊合情況，并不能完全代表復雜多變的實際地下地質情況，下面進一步分析局部異常與區域異常不同疊合關系的情形。

2.2 局部異常與區域異常疊合關系分類

在勘探地球物理領域，一般認為區域異常是由埋藏較深、水平延伸范圍較大的地質體或構造所產生的，在異常圖中表現出形態寬緩、幅值較大、范圍較廣的特征； 局部異常則是由埋藏較淺、水平延伸范圍較小的地質體所產生的，在異常圖中表現出形態陡窄、幅值較小、異常范圍小的特征。在遵循這一基本原則的前提下，本文針對四種實際地質情形，總結了相應的局部異常與區域異常的四種疊合關系分類情形，詳見表1和圖1。不難發現，式(11)僅適用于分類情形1和2，即區域異常介于總場異常值與周圍點平均值的取值區間。然而，對于區域異常與局部異常曲線彎曲方向相反的情形，即分類情形3和4，如仍使用式(11)計算切割區域場，計算結果雖可無限向周圍點平均值逼近，但依然忽略了圖1c和圖1d中黑色豎線所示部分的異常，這部分異常代表計算得到切割區域場與理論區域場的差異。

表1 局部異常與區域異常疊合關系分類表

圖1 對應表1中分類情形1～4(a～d)的局部異常與區域異常疊合關系分類

2.3 算法改進方案

為彌補傳統插值切割法的不足，本文提出了一種基于異常約束條件的改進方案，其基本思路是：首先，給定切割半徑r，分析區域異常和局部異常曲線彎曲方向； 其次，對區域異常與局部異常曲線彎曲方向相同的情形(即分類情形1和2)，采用式(11)計算切割區域場，對彎曲方向相反的情形(即分類情形3和4)，在式(11)基礎上，增加切割區域異常修正量ΔZ(x,y)(圖2)，即

圖2 切割區域異常修正量示意圖

(19)

首先考慮x方向異常剖面，假設局部異常附近的區域異常是曲率半徑為N、曲率中心為O的曲率圓的部分弧段(圖2中紅色曲線與藍色曲線重合部分)，則x方向切割區域的異常修正量為

(20)

(21)

(22)

式中Z′D，x、Z″D，x分別為Z(x,y)在D點對x的一階導數和二階導數。用一階差分替代一階導數，并忽略一階余量，得到

(23)

(24)

同理，可得到y方向切割區域的異常修正量ΔZ(y)，并取ΔZ(x,y)=[ΔZ(x)+ΔZ(y)]/2。

約束插值切割法詳細計算步驟見附錄A。以上算法基于MATLAB(R2008a)編程平臺實現程序和界面編制。

3 理論模型試算

3.1 理論模型設計

本文理論模型設計參考了程方道等[5]、Kaftan等[31]、Zhu等[40]及劉東甲[47]的設計思路，即深部大尺度密度體產生的異常視為區域異常，淺部小尺度密度體產生的異常視為局部異常，不同之處是在同一個理論模型中綜合考慮了局部異常與區域異常的四種疊合關系。鑒于2．2部分中分類情形1、3分別與2、4呈鏡像對稱關系，可合并處理，故本文設計了圖3所示的理論模型。三個不同深度的淺部球體(編號①、②、③)產生的是局部重力異常，其中淺部球體①、②產生局部重力高，淺部球體③產生局部重力低； 三個不同深度的深部球體(編號④、⑤、⑥)產生的是區域重力異常，其中球體④、⑥引起的是區域重力高異常，球體⑤引起的是重力低異常。

圖3 理論模型

3.2 理論模型正演

理論模型的重力異常為各球體剩余質量的引力位沿重力方向(定義為z方向)的導數之和[55]，即

(25)

基于式(25)，通過正演計算可得模型3的正演重力異常(圖4、圖5)。由圖4可見，局部異常基本淹沒在區域異常之中。由圖5可見，1、2、3號異常分別對應2．2部分中分類情形1、4、3。

圖4 圖3模型重力異常正演結果

圖5 圖3模型重力異常正演剖面(y=0)

3.3 理論模型重力場分離

分別采用插值切割法、約束插值切割法、小波多尺度分解法[13]對理論模型正演結果進行重力場分離，分離得到的區域場和局部異常見圖6。基于切割半徑與密度體埋深近似相等的對應關系，計算過程中以r=400m為初始值，以100m為增減量，采用多個切割半徑r，經對比篩選，確定取r=500m。根據小波多尺度分解的尺度場源深度轉換公式[14]，因網格間距為100m，故采用1～4階小波細節代表淺部球體產生的局部重力異常，用4階小波逼近代表深部球體產生的區域重力異常。對比約束插值切割法和小波多尺度分解法得到的局部異常，可見前者更接近于理論局部異常(圖4c)。由于區域異常幅值較大，基本淹沒了局部異常，故通過三種方法分離出的區域異常(圖6上)與理論區域異常(圖4b)基本一致。為了更直觀地分析本文方法的優越性，繪制了沿y=0剖面的三種方法分離得到的局部異常曲線(圖7)。由圖可見，由于1號異常為區域重力高疊加局部重力高(分類情形1)，三種方法的結果基本一致，都能較準確地體現分類情形1所產生的局部重力異常； 關于2號和3號異常，本文方法得到的分離結果明顯改善了插值切割法的計算結果，且與小波多尺度分解的結果吻合。三種方法均在局部異常兩翼附近產生了小幅度假異常，這是因為方法本身存在一定的計算誤差。以上分析證明了本文提出的約束插值切割法基本理論的正確性和分離效果的正確性。

圖6 圖3模型重力場不同方法分離得到的區域異常場(上)和局部異常場(下)

圖7 不同方法分離的局部重力異常場(y=0)

為了定量對比分析插值切割法、約束插值切割法和小波多尺度分解法的重力場分離效果，用下式計算局部異常ΔgⅠ與理論局部異常ΔgⅡ的相關系數

(26)

式中：Cor(·)表示協方差； Var(·)表示均方差。經計算，插值切割法、約束插值切割法和小波多尺度分解法得到的局部異常與理論局部異常的相關系數分別為0.42、0.75、0.70，可見本文提出的約束插值切割法得到的局部異常與理論局部異常具有最強的相關性。

另外，為了評價約束插值切割法對插值切割法的改善效果，引入評價航磁動態補償效果的“改善率”(Improvement Ratio，IR)[56]評價改善效果，其計算公式為

(27)

式中：Δgcz表示插值切割法得到的局部異常； Δgys表示約束插值切割法得到的局部異常； Std(·)表示求標準差。

因為本文提出的新方法僅對2號、3號異常進行了優化，故僅對2號、3號局部異常的改善率進行評估。經計算，約束插值切割法對插值切割法的改善率IR為1.48，可見相比于插值切割法，新方法計算的局部異常效果改善了48%，優化效果明顯。

4 應用實例

4.1 中國大陸及周邊區域重力數據分離

本文以中國大陸及周邊區域的重力數據為試驗對象，范圍為東經63°～136°，北緯15°～58°。重力數據來自美國國家地理空間情報局地球重力場研發小組發布的超高階地球重力場模型EGM2008(Earth Gravitational Model 2008)，該模型的階次完全至2159次，空間分辨率約為5′，可用于小比例尺重力編圖和構造研究[57]。本實例使用的布格重力異常的網格間距為4km。

分別利用插值切割法和本文提出的約束插值切割法對布格重力異常進行分離。為了展示新方法的正確性和改進效果，以葛粲等[48]的研究結果為參照。根據徐世浙等[8]、葛粲等[48]提出的“切割半徑應近似等于異常體埋深”的層切割思想，本例的切割半徑r以網格間距4km為基準，并以其整數倍的形式逐步遞增，每次切割得到的區域場作為下次切割的輸入，并結合前人經驗總結的“切割半徑r得到的局部場可以近似代表深度r以上的地層重力效應”，當切割半徑為19倍網格間距(r=76km)時，計算得到的局部異常可視為72～76km參考深度層的重力異常值，并將其視為上地幔頂部附近介質產生的重力異常。插值切割法和約束插值切割法計算得到的上地幔頂部附近重力異常分別見圖8和圖9。

對比圖8與葛粲等[48]的計算結果，可見重力異常分布特征基本一致，驗證了本文編程實現的插值切割法的正確性。對比插值切割法(圖8)與約束插值切割法(圖9)得到的上地幔頂部附近重力異常，可見整體上重力異常形態相似，但約束插值切割法得到的正值重力異常幅度和范圍相對更大，說明在負值背景下疊加正異常的情況下，約束插值切割法對異常分布進行了合理調整，這與約束插值切割法的基本原理和理論模型結論相一致。

圖8 插值切割法得到的中國大陸及周邊地區上地幔頂部附近重力異常場

圖9 約束插值切割法得到的中國大陸及周邊地區上地幔頂部附近重力異常場

為了定量對比約束插值切割法和插值切割法得到的重力異常，計算了兩者的差值(圖10)，即前文所述的修正量。由圖可見，在準噶爾盆地、喜馬拉雅地塊等局部地區，本文方法的修正達50%以上。

圖10 圖9與圖8的差值

根據圖9可得到初步的地質認識，即在喜馬拉雅地塊邊界局部地區(圖中綠色虛線橢圓區域)，深部高密度支撐物質強烈北移侵入，而柴達木盆地和準噶爾盆地正高值范圍的擴大(圖中白色虛線所示)，則意味著穩定盆地下方可能存在高密度物質。

楊文采等[58]利用小波多尺度分解法對中國陸域及周邊布格重力異常進行了分離，獲取了包括上地幔頂部異常(D7+D8)在內的4個深度界面所產生的重力異常。對比約束插值切割法與小波多尺度分解法對上地幔頂部異常的分離結果，可見在青藏高原等區域的異常分布特征與楊文采等[58]的研究成果相近，但約束插值切割法得到的異常細節更加豐富，主要是因為小波多尺度分解法僅將布格重力異常分離為4層，而約束插值切割法可分離為20層，對異常的分辨率更高[48]。

由于篇幅限制，結合地震、地質等資料的綜合地質解釋將另文詳細闡述。據本例分析可見，約束插值切割法在實際應用中實現了對插值切割法計算結果的修正，可得到更準確的重力異常、更明確的重力異常梯度帶。

4.2 中國南海萬安盆地及周邊區域重力數據分離

為了進一步驗證約束插值切割法的實際應用效果，選取中國南海萬安盆地及周邊區域布格重力異常數據為研究對象，范圍為東經104°～111°，北緯5°～11°，數據來源同4.1部分重力數據。陳玲等[59]認為南海萬安盆地新生界沉積層厚度最大為12.5km，馮旭亮等[60]的研究結果顯示該區域最大厚度大于10km，故本例約束插值切割法中選取切割半徑為3倍網格間距(12km)，經計算得到萬安盆地及周邊區域新生界沉積層的重力異常(圖11a)。圖11b為選取的一條地震剖面所對應的新生界沉積層重力異常，圖11c為該地震剖面構造解釋剖面[60]。

為了說明本文方法計算結果的可靠性，與馮旭亮等[60]采用最小曲率位場分離方法得到的計算結果(文獻[60]中的圖7、圖8)進行了對比，可以發現，圖11a中重力異常的整體分布特征與文獻[60]中的圖7較吻合，圖11b中的重力異常曲線走勢與文獻[60]中的圖8a大致相同，但本文的計算結果更能反映沉積基底的起伏變化細節，特別是更加直觀地顯示了F1斷層右側沉積基底的局部凸起(圖11c中的虛線圈所示)。因此，基于本文方法分離得到的局部重力異常，可更清晰地解釋新生界沉積層底界面的深度變化，為該區域新生界深度反演研究提供了較可靠的數據基礎，同時也驗證了約束插值切割法的實用性。

圖11 中國南海萬安盆地及周邊地區重力數據及地震剖面

5 結論

本文通過基本理論推導，并結合實際重力異常的疊合分類情形，實現了對插值切割法的優化和改進，提出了基于異常約束條件的插值切割法(約束插值切割法)，并開展了理論模型試驗和衛星重力數據處理、分析。得出以下幾點認識。

(1)通過對插值切割法切割算子的公式推導，結合四種情形的異常疊合關系，證明了插值切割法的切割算子僅考慮了局部異常與區域異常彎曲方向相同的疊合情況。基于此，給出了插值切割算子的補償方法，提出了基于異常約束的插值切割法，制定了新方法的計算流程，并基于Matlab平臺編程實現。

(2)理論模型試算結果顯示，插值切割法、約束插值切割法、小波多尺度分解法中約束插值切割法得到的局部異常與理論局部異常具有最強的相關性(相關系數為0.75)，與插值切割法相比，對局部異常效果的改善達48%，優化效果明顯。

(3)對中國陸域及周邊區域的布格重力異常進行了分離，獲取了72～76km參考深度層產生的重力異常，約束插值切割法獲取了該層位更真實可靠的重力異常值信息，在準噶爾盆地、喜馬拉雅地塊等局部地區的修正達50%以上，可為上地幔頂部構造特征研究提供參考。

(4)利用約束插值切割法分離得到了萬安盆地及周邊區域新生界沉積層重力異常，與地震剖面顯示的基底起伏特征基本一致，說明該方法可用于基底深度變化的反演研究。

總之，通過理論模型數據和實際重力數據分析，約束插值切割法計算結果均明顯優化了插值切割法的計算結果，驗證了前者理論的正確性及實用性，可為實際重力數據的地質解釋提供更準確、可靠的重力異常信息。本文提出的新方法也可用于磁場數據的異常分離。

在構思本文的初期，與插值切割法的原創者——中國石油勘探開發研究院文百紅老師進行了有益探討，特別是在插值切割法的細節上得到了文百紅老師的悉心指導，在此深表謝意； 感謝中國自然資源航空物探遙感中心駱遙教授為文章的修改提出了建設性的意見！

附錄A 約束插值切割法計算步驟

(1)根據目標層位深度信息，對重力異常Z(x,y)設置切割半徑r；

(2)利用式(2)計算與點(x,y)相距r的某4個點重力異常的平均值M(x,y)；

(3)利用式(10)計算加權參數e;

(4)比較區域異常與局部異常曲線彎曲方向，具體可細分為2步：

①計算下列參數

②當比較參數不能同時滿足下列條件時，表明區域異常和局部異常彎曲方向相同(即分類情形1和2)

Ix1Iy1>0

Iy1Iy2<0

Ix1Ix2<0

Iy2Ix2>0：

當比較參數同時滿足下列條件時，表明區域異常和局部異常彎曲方向相反(即分類情形3和4)

Ix1Iy1>0

Iy1Iy2<0

Ix1Ix2<0

Iy2Ix2>0：

(5)對于彎曲方向相同的情形，區域場等于Z(x,y)與M(x,y)的加權平均，采用式(11)計算切割區域場；

(6)對于彎曲方向相反的情形，區域場為Z(x,y)與M(x,y)的加權平均的基礎上增加切割區域異常修正量ΔZ(x,y)，此種情況可細分為以下4步：

①首先考慮x方向異常剖面，假設局部異常附近的區域異常為曲率半徑為N、曲率中心為O的曲率圓的部分弧段，利用式(20)計算x方向切割區域異常修正量ΔZ(x)。

②同理可得y方向切割區域異常修正量ΔZ(y)；

③取ΔZ(x,y)=[ΔZ(x)+ΔZ(y)]/2；

④計算修正量ΔZ(x,y)后，采用式(19)計算切割區域場；

(7)用上述方法得到第1次切割的區域場R1(x,y)，對R1(x,y)重復使用步驟(1)～(6)，得到第2次切割的區域場R2(x,y)，依次迭代，直到滿足式(12)，迭代終止，得到切割區域場；

(8)利用式(14)計算最終的局部場R(x，y)。