一種精妙的數(shù)據(jù)重建方法：網(wǎng)函數(shù)插值法

陳啟宏， 張兆田， 邱 鈞

(1.上海財(cái)經(jīng)大學(xué)， 上海 200433; 2.國家自然科學(xué)基金委員會(huì)， 北京 100085; 3.北京信息科技大學(xué)， 北京 100101)

0 引言

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

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

1 網(wǎng)函數(shù)插值法

簡而言之，網(wǎng)函數(shù)插值法是由數(shù)學(xué)上直紋面構(gòu)造曲面的一種方法。它由不同的三個(gè)直紋面構(gòu)成，其中兩個(gè)是作為近似曲面的直紋面，第三個(gè)是作為全局調(diào)整的直紋面。不失一般性，以二元網(wǎng)函數(shù)插值為例，闡述具體構(gòu)造算法。

1.1 二元網(wǎng)函數(shù)插值及其余項(xiàng)

在xOy平面上一閉矩形區(qū)域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)

假設(shè)函數(shù)f(x,y)在矩形區(qū)域D上是連續(xù)的。可用如下方法，由函數(shù)f(x,y)在矩形區(qū)域D邊界上的值，計(jì)算函數(shù)f在D內(nèi)任一點(diǎn)的近似值。

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

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

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

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

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

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

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

最后，由三張直紋面構(gòu)造的網(wǎng)函數(shù)插值曲面為：z=F(x,y)。其中

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

(5)

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

由此可見，區(qū)域D面積A越小，函數(shù)f與網(wǎng)函數(shù)插值法函數(shù)F的逼近程度越高。偏導(dǎo)數(shù)最大值M與函數(shù)f在區(qū)域D上的起伏大小有關(guān)。對于在區(qū)域D上變化較平緩的函數(shù)，網(wǎng)函數(shù)插值函數(shù)與f的逼近程度也就較高。

1.2 矩形網(wǎng)函數(shù)插值的算子表示

上述近似曲面構(gòu)造過程中,F3由F1和F2計(jì)算得來的，進(jìn)而網(wǎng)插值函數(shù)F可用算子簡潔表示：

F=F1?F2=F1+F2-F1F2

(6)

式中，?為布爾和。

1.3 相應(yīng)的偏微分方程邊值問題的解空間

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

(7)

這里函數(shù)F的解空間可表示為：

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

(8)

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

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

1.4 矩形網(wǎng)函數(shù)插值的面積加權(quán)簡潔計(jì)算格式

利用上述三張直紋面構(gòu)造的矩形網(wǎng)函數(shù)插值，可得到其與面積有關(guān)的簡潔計(jì)算格式。設(shè)D是xOy的平面上的一矩形，D的四個(gè)角點(diǎn)記作Pi(i=1,2,3,4)。假設(shè)f是定義在D上的連續(xù)函數(shù)，過Q點(diǎn)作平行于D的四條邊的兩條平行線，交四條邊于四個(gè)點(diǎn)記作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)

因此，網(wǎng)函數(shù)插值法表示為：

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

(12)

從計(jì)算格式上看出，區(qū)域D上任一點(diǎn)Q的計(jì)算值，均由在區(qū)域邊界上的相應(yīng)8個(gè)已知值線性表示。其中，線性系數(shù)是相應(yīng)小矩形面積與大矩形面積的比值，取決于點(diǎn)Q的幾何位置，其與坐標(biāo)系的選擇無關(guān)[1]。因此，該算法簡單，易于計(jì)算。

2 網(wǎng)函數(shù)插值法應(yīng)用案例

2.1 網(wǎng)函數(shù)插值法在重力勘查中的應(yīng)用

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

網(wǎng)函數(shù)插值法能夠成功地應(yīng)用于地質(zhì)勘探領(lǐng)域[2-3]，主要得益于兩點(diǎn)：其一是該算法結(jié)構(gòu)十分簡潔，便于計(jì)算機(jī)實(shí)現(xiàn)；其二是該算法適用于由良好的解空間包含的函數(shù)類，適合于分離大小不一、形態(tài)各異的局部異常。改進(jìn)后的方法為分離更高階次的區(qū)域場創(chuàng)造了可能性，便于分離大局部異常中更次級的小局部異常。該方法保證精度的關(guān)鍵是，使四條邊界線上待分離局部異常和旁側(cè)局部異常的影響最小。

2.2 網(wǎng)函數(shù)插值法在植物群落種群分布格局中的應(yīng)用

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

3 相關(guān)研究進(jìn)展

(1)鑒于地貌環(huán)境復(fù)雜、勘探測線多樣性，邱佩璋先生研究團(tuán)隊(duì)把矩形網(wǎng)函數(shù)插值法擴(kuò)展到三角網(wǎng)函數(shù)插值法，更靈活適應(yīng)勘查重力場計(jì)算[7]。

(2)基于函數(shù)類特點(diǎn)和目標(biāo)函數(shù)先驗(yàn)知識(shí)，發(fā)展了基于數(shù)據(jù)驅(qū)動(dòng)和模型驅(qū)動(dòng)的網(wǎng)函數(shù)插值方法，即改進(jìn)的網(wǎng)函數(shù)插值法[8]。若預(yù)先測知函數(shù)f的形態(tài)，但其數(shù)學(xué)表達(dá)式不符合網(wǎng)函數(shù)法所限定的函數(shù)類型。如果可通過數(shù)學(xué)變換T，使得T(f)成為(近似)符合網(wǎng)函數(shù)插值法所限定的函數(shù)類型，則網(wǎng)函數(shù)F計(jì)算的值T(f)是(近似)精確的，再對F進(jìn)一步做數(shù)學(xué)逆變換T-1就可得到待插函數(shù)f。

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

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

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

(6)CT是由系列探測網(wǎng)格點(diǎn)上的探測數(shù)據(jù)，利用切片定理重建目標(biāo)圖像。如果探測數(shù)據(jù)不完全，在一定條件下，也可通過網(wǎng)函數(shù)插值方法進(jìn)行數(shù)據(jù)重建，重構(gòu)出近似的目標(biāo)圖像。

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

致謝

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