一種精妙的數據重建方法：網函數插值法

陳啟宏， 張兆田， 邱 鈞

(1.上海財經大學， 上海 200433; 2.國家自然科學基金委員會， 北京 100085; 3.北京信息科技大學， 北京 100101)

0 引言

數據重建通常是由函數在子流形的相關信息，重建其在高維流形的函數值。CT(Computed Tomography)就是由函數在系列直線上的積分構成的低維投影數據集，重建高維函數。網函數插值是由函數在高維流形邊界上的探測數據，重建高維流形上的函數(近似)值。

網函數插值問題是由勘查重力場數據尋找礦藏實際需求驅動的。重力場通常是由大范圍的背景區域重力場和由局部地質體(含礦體)引起的局部異常重力場構成。重力勘查不可缺少的步驟之一就是，依據野外獲取的大量網格點上重力探測數據，分離區域重力場，進而獲得由局部地質體(含對尋找油氣和煤田有重要意義的地質構造和規模巨大的礦體)引起的局部異常重力場分布[1]。

1 網函數插值法

簡而言之，網函數插值法是由數學上直紋面構造曲面的一種方法。它由不同的三個直紋面構成，其中兩個是作為近似曲面的直紋面，第三個是作為全局調整的直紋面。不失一般性，以二元網函數插值為例，闡述具體構造算法。

1.1 二元網函數插值及其余項

在xOy平面上一閉矩形區域D={(x,y)|x0≤x≤x1,y0≤y≤y1}上，其四條邊界曲線記為：

u0={(x0,y),f(x0,y)|y0≤y≤y1}

(1)

u1={(x1,y),f(x1,y)|y0≤y≤y1}

(2)

λ0={(x,y0),f(x,y0)|x0≤x≤x1}

(3)

λ1={(x,y1),f(x,y1)|x0≤x≤x1}

(4)

假設函數f(x,y)在矩形區域D上是連續的?？捎萌缦路椒?，由函數f(x,y)在矩形區域D邊界上的值，計算函數f在D內任一點的近似值。

圖1 曲面的四條邊界線z=F(x,y)：μ0,μ1,λ0,λ1

圖2 直紋面：z=F1(x,y)

這里從幾何角度，通過三張直紋面在矩形區域D上構造一曲面F(x,y)，作為函數f(x,y)的近似(以下函數和曲面視為相同含義)。

第一步：設平面y=yq與曲線u0、u1相交于兩點，將這兩點用直線段聯結起來。對于[y0,y1]中任意一數yq，均如上方公式作直線段，這樣得到第一張直紋面，記作z=F1(x,y)。

第二步：對于[x0,x1]中任意一點xq，設平面x=xq與曲線λ0、λ1相交于兩點，將這兩點用直線段聯結起來，這樣構成第二張直紋面，記作z=F2(x,y)。

圖3 直紋面z=F2(x,y)

圖4 直紋面z=F3(x,y)

最后，由三張直紋面構造的網函數插值曲面為：z=F(x,y)。其中

F(x,y)=F1(x,y)+F2(x,y)-F3(x,y)

(5)

不難看出，近似曲面F(x,y)在區域D的邊界上與f(x,y)取值相等，而直紋面F1(x,y)、F2(x,y)保留了在區域D的邊界上與f(x,y)比較接近的部分，用直紋面F3(x,y)在區域D上整體調整F1(x,y)+F2(x,y)與f(x,y)偏差較大的部分。

由此可見，區域D面積A越小，函數f與網函數插值法函數F的逼近程度越高。偏導數最大值M與函數f在區域D上的起伏大小有關。對于在區域D上變化較平緩的函數，網函數插值函數與f的逼近程度也就較高。

1.2 矩形網函數插值的算子表示

上述近似曲面構造過程中,F3由F1和F2計算得來的，進而網插值函數F可用算子簡潔表示：

F=F1?F2=F1+F2-F1F2

(6)

式中，?為布爾和。

1.3 相應的偏微分方程邊值問題的解空間

假設F(x,y)在區域D的邊界上二次連續可微，則由三張直紋面構造的函數F滿足偏微分方程：

(7)

這里函數F的解空間可表示為：

F(x,y)=φ(x)+yφ1(x)+ψ(y)+xψ1(y)

(8)

φ(x)、φ1(x)在區間[x0,x1]，ψ(y)、ψ1(y)在區間[y0,y1]上均為二次可微函數[1]。從上式可見，偏微分方程邊值問題解F的解空間具有廣泛的選擇性，可解釋為具有良好的函數逼近能力。

若函數f屬于網函數插值法所限定的函數類型，則網插值函數F計算的值是精確的。例如，f為三次多項式函數，則F=f。正是由于重力區域異常通?？捎刹桓哂谌蔚亩囗検綌M合，符合網函數插值法的應用條件，因而在重力勘查數據處理中取得了理想效果。

1.4 矩形網函數插值的面積加權簡潔計算格式

利用上述三張直紋面構造的矩形網函數插值，可得到其與面積有關的簡潔計算格式。設D是xOy的平面上的一矩形，D的四個角點記作Pi(i=1,2,3,4)。假設f是定義在D上的連續函數，過Q點作平行于D的四條邊的兩條平行線，交四條邊于四個點記作Qj(j=1,2,3,4)。

過Q的兩條平行線將面積為A的矩形分成四塊小矩形，小矩形面積記作Ai(i=1,2,3,4)，如圖5。

圖5 矩形剖分

則上述三張直紋面可分別表示為：

F1(Q)=[(A2+A3)f(Q2)+(A4+A1)f(Q4)]/A

(9)

F2(Q)=[(A1+A2)f(Q1)+(A3+A4)f(Q3)]/A

(10)

F3(Q)=[A1f(P1)+A2f(P2)+A3f(P3)+A4f(P4)]/A

(11)

因此，網函數插值法表示為：

F(Q)=F1(Q)+F2(Q)-F3(Q)

(12)

從計算格式上看出，區域D上任一點Q的計算值，均由在區域邊界上的相應8個已知值線性表示。其中，線性系數是相應小矩形面積與大矩形面積的比值，取決于點Q的幾何位置，其與坐標系的選擇無關[1]。因此，該算法簡單，易于計算。

2 網函數插值法應用案例

2.1 網函數插值法在重力勘查中的應用

20世紀 70 年代中期至 80 年代初，針對內蒙古地質局劉士毅先生提及關于勘查重力場數據處理的實際需求，內蒙古大學邱佩璋先生帶領團隊與內蒙古物探隊劉士毅、郭積忠等緊密合作，首先將網函數插值法用于重力勘查區域場校正。當時對鄂爾多斯含油氣盆地北部 750 km2的 1∶200 000 重力普查資料，用網函數插值法等作區域校正。與通常采用的圓周法和趨勢分析法相比，網函數插值法計算結果表明，共 35 個大小不同的局部異常全都得到了突出，異常特征更清晰、形態更完整。相關理論計算結果后期被實地勘探所證實。網函數插值法的成功創新應用，為重力勘查做出了重要貢獻，相關工作被授予地質礦產部科技進步獎。

網函數插值法能夠成功地應用于地質勘探領域[2-3]，主要得益于兩點：其一是該算法結構十分簡潔，便于計算機實現；其二是該算法適用于由良好的解空間包含的函數類，適合于分離大小不一、形態各異的局部異常。改進后的方法為分離更高階次的區域場創造了可能性，便于分離大局部異常中更次級的小局部異常。該方法保證精度的關鍵是，使四條邊界線上待分離局部異常和旁側局部異常的影響最小。

2.2 網函數插值法在植物群落種群分布格局中的應用

20世紀80年代中期，內蒙古大學數學和生物學研究人員楊在中、郝敦元和楊持等將網函數插值法應用于生態研究中，給出一種研究植物群落種群分布格局的新方法[4]。實踐表明，網函數插值法可有效地應用于生物群落種群分布格局的研究[4-5]。新方法避開了英國著名生態學家 Greig-Smith 提出“鄰接格子樣方法”面臨的數據統計量大、難于計算等問題[6]，在實踐中得到推廣應用[3]。

3 相關研究進展

(1)鑒于地貌環境復雜、勘探測線多樣性，邱佩璋先生研究團隊把矩形網函數插值法擴展到三角網函數插值法，更靈活適應勘查重力場計算[7]。

(2)基于函數類特點和目標函數先驗知識，發展了基于數據驅動和模型驅動的網函數插值方法，即改進的網函數插值法[8]。若預先測知函數f的形態，但其數學表達式不符合網函數法所限定的函數類型。如果可通過數學變換T，使得T(f)成為(近似)符合網函數插值法所限定的函數類型，則網函數F計算的值T(f)是(近似)精確的，再對F進一步做數學逆變換T-1就可得到待插函數f。

(3)網函數插值法可推廣到n維空間(m)網[9]。例如，一維網函數插值與二點線性插值、三點拋物插值等密切相關。這給出一種構造曲線或曲面的新視角。

(4)利用網函數插值生成 Coons 型分形曲面，曲面形態自然并可保持分形特征[10]。

(5)網函數插值方法可有效地重建探測過程中的不完全數據或在數據傳輸過程中的局部信息缺損[11-12]。

(6)CT是由系列探測網格點上的探測數據，利用切片定理重建目標圖像。如果探測數據不完全，在一定條件下，也可通過網函數插值方法進行數據重建，重構出近似的目標圖像。

注：邱佩璋先生1951年畢業于上海交通大學數學系，進入中科院數學所，師從吳新謀先生，研究偏微分方程基本解、Hadamard理論等。1957年響應黨和國家號召支援邊疆，參與組建內蒙古大學數學系。網函數插值思想萌發在70年代，在發現不完全Huyghens現象，揭示了基本解升維結構，及其展開系數對稱性引起的奇特性質，進而形成了結點網產生器的原始思路和方法。當年由重力勘查數據處理的實際探尋礦藏這個重大需求而牽引，邱佩璋先生領導的理論計算團隊與劉士毅先生帶領的內蒙物探實踐團隊經過多年的合作研究與迭代演進，網函數插值法在80年代逐步完善，并在區域場校正等等實踐中發揮了重要作用。鑒于網函數插值法在重力勘查中的成功運用，1980年，邱佩璋先生等獲得了內蒙古自治區的科技進步三等獎，并榮獲內蒙古自治區先進科技工作者稱號。在邱佩璋先生逝世周年之際，謹以此文緬懷先生教誨。

致謝

感謝中國地質調查局發展研究中心劉士毅研究員在本文撰寫過程中的指導、修訂和支持，感謝北京信息科技大學劉暢老師等對該文稿的修訂整理。