楊博，曹書錦,，毛雅靜，周勇，朱自強(qiáng)，陳新躍

1.湖南科技大學(xué) 頁(yè)巖氣資源利用湖南重點(diǎn)實(shí)驗(yàn)室，湖南 湘潭 411201；2.湖南科技大學(xué) 資源環(huán)境與安全工程學(xué)院，湖南 湘潭 411201；3.中南大學(xué) 地球科學(xué)與信息物理學(xué)院，長(zhǎng)沙 410083；4.湖南文理學(xué)院 資源環(huán)境與旅游學(xué)院，湖南 常德 415000

0 引言

歐拉反褶積法能使用少量先驗(yàn)信息且不需要準(zhǔn)確的解釋模型，即可快速準(zhǔn)確地圈定異常源的基本輪廓或標(biāo)識(shí)異常源中心，特別適合大尺度重磁位場(chǎng)數(shù)據(jù)的快速解釋。歐拉反褶積法一般通過(guò)一定大小的滑動(dòng)窗口在整個(gè)勘探區(qū)域移動(dòng)，根據(jù)人工經(jīng)驗(yàn)或枚舉構(gòu)造指數(shù)，在滑動(dòng)窗口內(nèi)基于歐拉齊次方程構(gòu)建線性方程組[1--7]。該線性方程組由人為設(shè)定大小的滑動(dòng)窗口內(nèi)的觀測(cè)數(shù)據(jù)構(gòu)建，其條件數(shù)極大，這導(dǎo)致歐拉反褶積易于受到觀測(cè)數(shù)據(jù)噪聲、異常源形態(tài)及其空間位置和滑動(dòng)窗口尺寸等因素的影響。

在實(shí)際地質(zhì)情況中，地質(zhì)構(gòu)造一般為多個(gè)不同形態(tài)異常源的組合，傳統(tǒng)歐拉反褶積解釋工作流程中僅使用單一構(gòu)造指數(shù)易于導(dǎo)致歐拉解發(fā)散，而通過(guò)枚舉構(gòu)造指數(shù)易于導(dǎo)致歐拉解集過(guò)于龐大[8--14]，單一過(guò)濾標(biāo)準(zhǔn)又將忽視其他異常源[1]。張季生等[15]針對(duì)歐拉反褶積方法易受高頻干擾的問(wèn)題，提出低階的歐拉反褶積與解析信號(hào)相結(jié)合的位場(chǎng)反演方法，將高頻噪聲的干擾減少到可以控制的水平[16]。此外，將重力梯度張量解析信號(hào)引入歐拉反褶積，或?qū)⒅亓μ荻葟埩繑?shù)據(jù)進(jìn)行聯(lián)合歐拉三維反演，以提高歐拉反褶積解的收斂性[17--18]。

對(duì)于歐拉反褶積在具體實(shí)施過(guò)程中最優(yōu)化窗口選擇的問(wèn)題，周文月探索了小窗口多數(shù)據(jù)對(duì)歐拉反褶積方法的影響，在一定程度上解決了采用大窗口時(shí)解的平滑問(wèn)題，以及不同場(chǎng)源混疊相互干擾引起的解的不精確問(wèn)題[19]。Krahenbuhl et al.[20]指出過(guò)小的滑動(dòng)窗口難以獲得地下異常源信息，而過(guò)大的滑動(dòng)窗口易于引入相鄰異常源的干擾，但沒(méi)有給出如何選擇最優(yōu)滑動(dòng)窗口大小的策略。研究最優(yōu)滑動(dòng)窗口大小確定方法，可以有效地捕捉地下異常源信息、減少歐拉反褶積解集大小和提高工作效率[1, 21--22]。

為此，針對(duì)超定的歐拉齊次方程組易于受到多種因素的干擾，筆者引入截?cái)嗾`差對(duì)歐拉反褶積解的優(yōu)良性進(jìn)行評(píng)價(jià)，并引入模糊聚類對(duì)歐拉反褶積解集進(jìn)行劃分[3]，在此基礎(chǔ)上，當(dāng)異常源尺寸一定時(shí)，以不同深度的異常源、不同尺寸的滑動(dòng)窗口和歐拉解集中優(yōu)解的比例等三者為研究對(duì)象，分析滑動(dòng)窗口尺寸和異常源尺寸間的比例關(guān)系，以此選擇最優(yōu)化窗口。最后，數(shù)值算例進(jìn)一步表明了方法的有效性。

1 歐拉反褶積

重力場(chǎng)歐拉齊次方程可以寫為[23]：

fx(x-x0)+fy(y-y0)+fz(z-z0)=Nf(x-x0,y-y0,z-z0)+NB

(1)

式中：f為位于點(diǎn)(x，y，z)的異常源在觀測(cè)點(diǎn)(x0，y0，z0)處產(chǎn)生的異常響應(yīng)；fx、fy和fz為f沿x、y和z軸向的梯度，可通過(guò)離散傅立葉變換或離散余弦變換獲得；N為異常幅值隨場(chǎng)源深度衰減變化“陡緩”的量度，即構(gòu)造指數(shù)；B為正常場(chǎng)存在的偏差，主要是為消除區(qū)域背景場(chǎng)的影響而引入的一個(gè)代表背景場(chǎng)的參數(shù)。在傳統(tǒng)歐拉反褶積法中，一般通過(guò)人工經(jīng)驗(yàn)選擇或枚舉構(gòu)造指數(shù)N求解式(1)，但單一構(gòu)造指數(shù)N易于導(dǎo)致歐拉解發(fā)散或無(wú)解；而在一定范圍內(nèi)，如0～3，以某一步長(zhǎng)枚舉構(gòu)造指數(shù)N易使歐拉解解集過(guò)于龐大，導(dǎo)致解釋工作量劇增。為此，將構(gòu)造指數(shù)N作為待求參數(shù)，將式(1)重寫為：

xfx+yfy+zfz+Nf=x0fx+y0fy+z0fz

(2)

由于B數(shù)量級(jí)很小，故此，在式(2)中已將B略去。在歐拉反褶積中，一般利用長(zhǎng)寬均為wx的滑動(dòng)窗口在勘探區(qū)內(nèi)移動(dòng)。對(duì)于每個(gè)滑動(dòng)窗口而言，以滑動(dòng)窗口內(nèi)的n=wx×wx個(gè)觀測(cè)值根據(jù)式(2)構(gòu)建線性方程組，并求解該方程組以獲得歐拉解(x，y，z，N)。

將構(gòu)造指數(shù)N當(dāng)作未知數(shù)求解可在一定程度上避免歐拉解發(fā)散，但式(2)本質(zhì)上仍然是求解超定的歐拉齊次方程，因而歐拉解仍然存在發(fā)散的趨勢(shì)。為此，以截?cái)嗥娈愔捣纸夥ǖ玫降恼`差來(lái)估計(jì)歐拉解的優(yōu)良性，以聚類分析來(lái)劃分歐拉解解集，并評(píng)價(jià)歐拉解解集中各簇的優(yōu)良性。在此基礎(chǔ)上，分析滑動(dòng)窗口的尺寸與地下異常源尺寸的關(guān)系，研究最優(yōu)滑動(dòng)窗口大小選擇標(biāo)準(zhǔn)，以滿足在反演地下不同尺寸異常源的同時(shí)，減少相鄰干擾信息的侵入，控制歐拉解集的大小，降低歐拉反演的后續(xù)處理工作量。

2 截?cái)嗾`差估計(jì)

將滑動(dòng)窗口內(nèi)n個(gè)觀測(cè)值f代入到式(2)所得到的線性方程組寫為矩陣的形式：

Gm=b

(3)

式中：G為[fxfyfzf]；b為x0fx+y0fy+z0fz構(gòu)成的向量；m=[x，y，z，N]T為歐拉解。

式(3)為超定方程組，因而在勘探區(qū)內(nèi)移動(dòng)滑動(dòng)窗口時(shí)，歐拉反褶積可能因?yàn)槭?3)條件數(shù)過(guò)大而導(dǎo)致無(wú)解(即解為非數(shù)的情況)，或有解但不為正確解/優(yōu)解(如歐拉解的空間位置不在勘探區(qū)內(nèi)、構(gòu)造指數(shù)不在合理范圍內(nèi)等)。為此，這里采用截?cái)嗥娈愔捣纸夥▽?duì)式(3)的誤差進(jìn)行估計(jì)，并在此基礎(chǔ)上，構(gòu)建一個(gè)估計(jì)函數(shù)εr以過(guò)濾歐拉解集中的繆解：

(4)

式中：ε為截?cái)嗥娈愔捣纸夥▽?duì)式(3)的誤差估計(jì)值。

3 模糊聚類原理

為判定不同大小的滑動(dòng)窗口反演地下異常源的能力，特引入聚類分析方法將歐拉解集按相似性分類，以標(biāo)示地下多個(gè)不同的異常源。在多種聚類分析方法中，基于劃分的模糊聚類分析不需要訓(xùn)練樣本，僅通過(guò)無(wú)監(jiān)督學(xué)習(xí)達(dá)到自動(dòng)分類的目的，最終能夠獲得比較客觀地反映地質(zhì)構(gòu)造的簇。模糊聚類不但能從原始數(shù)據(jù)中提取特征，而且能在提取特征之后，進(jìn)一步提供最近鄰原型分類器和提取空間劃分特性等。為把p個(gè)向量xi(i=1，2，…，p)分類成c個(gè)簇，并使得目標(biāo)函數(shù)最小，求得簇的分類信息和簇的中心，以模糊聚類算法構(gòu)建如下目標(biāo)函數(shù)：

(5)

(6)

(7)

由模糊聚類得到的歐拉解簇的聚集程度可由下式計(jì)算得到：

(8)

式中：ni為第i個(gè)歐拉解簇中元素的個(gè)數(shù)；ni,εr為第i個(gè)歐拉解簇中優(yōu)解的個(gè)數(shù)。

4 模型試驗(yàn)

4.1 算法驗(yàn)證

為驗(yàn)證本文歐拉反褶積算法和基于截?cái)嗥娈愔捣纸夥ㄔO(shè)計(jì)的估計(jì)函數(shù)的正確性，特基于有正演解析解，且構(gòu)造指數(shù)易于確定的立方體模型，構(gòu)建如下地球物理模型I：在地下均勻半空間內(nèi)，存在一個(gè)大小為500 m×500 m×500 m、頂板埋深為200 m、底板埋深為700 m且剩余密度為0.3 g/cm3的異常體。以大小50 m×50 m×50 m的網(wǎng)格離散該地球物理模型，網(wǎng)格沿著笛卡爾坐標(biāo)系三軸向的剖分?jǐn)?shù)nx、ny和nz分別為32、32和16；觀測(cè)點(diǎn)間距為50 m×50 m，觀測(cè)點(diǎn)個(gè)數(shù)為nx×ny=32×32=1 024，即勘探區(qū)尺寸為1 550 m×1 550 m，觀測(cè)點(diǎn)高度為地面上25 m。

首先，通過(guò)解析解計(jì)算立方體異常源的重力異常響應(yīng)f，并在f中添加3%的高斯噪聲作為觀測(cè)數(shù)據(jù)，繼而利用快速傅立葉變換計(jì)算得到fx、fy和fz。在此基礎(chǔ)上，應(yīng)用式(3)計(jì)算異常源的空間位置及其構(gòu)造指數(shù)，利用聚類分析對(duì)歐拉解集進(jìn)行分類，并采用估計(jì)函數(shù)εr對(duì)歐拉解的優(yōu)良性進(jìn)行判定。

圖1為孤立異常源的歐拉解分布示意圖，滑動(dòng)窗口尺寸wx為5。從圖中可以看出構(gòu)造指數(shù)趨近于2(立方體異常源構(gòu)造指數(shù)的理論值為2)的歐拉解聚集于異常源中心，這表明了本文算法的正確性。

圖1 孤立異常源歐拉解分布示意圖Fig.1 Scatter plot of Euler solutions for a single anomalous sources

進(jìn)一步，以εr>0.25為標(biāo)準(zhǔn)過(guò)濾孤立異常源的歐拉解，其分布如圖2所示。通過(guò)與圖1進(jìn)行對(duì)比，可以發(fā)現(xiàn)基于εr過(guò)濾后的歐拉解能有效的標(biāo)示異常源的空間位置，而且盡可能的保留了一些構(gòu)造指數(shù)偏離理論值的歐拉解，如N≤1.5的歐拉解。這表明本研究提出的估計(jì)函數(shù)的正確性。

圖2 孤立異常源的歐拉解以εr>0.25為標(biāo)準(zhǔn)過(guò)濾后的分布示意圖Fig.2 Scatter plot of Euler solutions for a single anomalous sources filtered by εr>0.25

4.2 不同尺寸滑動(dòng)窗口提取不同深度孤立異常源的能力分析

為研究不同尺寸滑動(dòng)窗口提取不同深度孤立異常源的能力，在其他參數(shù)與地球物理模型I保持一致的前提下，將地球物理模型I以步長(zhǎng)Δz沿Z軸向下移動(dòng)從而生成一系列不同深度的孤立異常源模型II，其中勘探區(qū)尺寸與4.1小節(jié)中一致。

解率為歐拉解集中不為非數(shù)的解與歐拉解個(gè)數(shù)的百分比。在歐拉反褶積中，滑動(dòng)窗口wx最小為3。在勘探區(qū)域內(nèi)，滑動(dòng)窗口沿著x和y軸移動(dòng)的最大步數(shù)分別為nx-wx和ny-wx，滑動(dòng)窗口移動(dòng)次數(shù)最大為(nx-wx)×(ny-wx)，且歐拉解集中解的個(gè)數(shù)小于或等于這一數(shù)字。因而，對(duì)于固定大小的勘探區(qū)域，隨著滑動(dòng)窗口尺寸不斷增大，解的數(shù)目在不斷的減少。

由圖3可見(jiàn)，當(dāng)滑動(dòng)窗口較小時(shí)，即滑動(dòng)窗口尺寸小于勘探區(qū)域的1/4時(shí)，淺部異常源對(duì)應(yīng)的有解率比深部的要大；當(dāng)滑動(dòng)窗口尺寸約為勘探區(qū)域的1/4～1/2時(shí)，深部異常源對(duì)應(yīng)的有解率比淺部的要大；而當(dāng)滑動(dòng)窗口尺寸大于勘探區(qū)域尺寸1/2時(shí)，不同深度的異常源對(duì)應(yīng)的有解率幾乎一致。

圖3 不同尺寸滑動(dòng)窗口下不同深度孤立異常源的有解率曲線Fig.3 Different moving window sizes for Euler solution rate of isolated anomalous sources with various offset Δz

優(yōu)解率為歐拉解集中優(yōu)解(排除空間位置不在勘探區(qū)內(nèi)、構(gòu)造指數(shù)不在合理范圍內(nèi)的歐拉解)與歐拉解集中不為非數(shù)的解的個(gè)數(shù)的百分比。圖4為利用估計(jì)函數(shù)εr對(duì)不同深度孤立異常源的歐拉解進(jìn)行過(guò)濾后的優(yōu)解率曲線。當(dāng)滑動(dòng)窗口很小時(shí)，淺部異常源對(duì)應(yīng)的歐拉優(yōu)解率占優(yōu)；隨著滑動(dòng)窗口的增大，這一關(guān)系仍然得以保持；當(dāng)滑動(dòng)窗口尺寸大于勘探區(qū)域的1/2時(shí)，不同深度異常源對(duì)應(yīng)的優(yōu)解率趨于一致。

圖4 不同尺寸滑動(dòng)窗口下不同深度孤立異常源的優(yōu)解率曲線(3%高斯噪聲)Fig.4 Different moving window sizes for optimal Euler solution rate of isolated anomalous sources with various offset Δz, survey data contaminated with 3% Gaussian noise

在圖4中試驗(yàn)的基礎(chǔ)上，進(jìn)一步在觀測(cè)數(shù)據(jù)中添加7%的高斯噪聲，以對(duì)比不同噪聲情況下，不同尺寸滑動(dòng)窗口提取不同深度孤立異常源的能力(圖 5)。當(dāng)滑動(dòng)窗口很小時(shí)，幾乎難于提取深部異常源的有用信息；隨著滑動(dòng)窗口的不斷增大，歐拉反褶積優(yōu)解率曲線的斜率先增大后放緩，且這種趨勢(shì)在深部異常上更為明顯。相比于圖4，隨著孤立異常源的觀測(cè)數(shù)據(jù)中噪聲的增加，將導(dǎo)致基于小滑動(dòng)窗口的歐拉反褶積難于提取到有用信息。

圖5 不同尺寸滑動(dòng)窗口下不同深度孤立異常源的優(yōu)解率曲線(7%高斯噪聲)Fig.5 Different moving window sizes for optimal Euler solution rate of isolated anomalous sources with various offset Δz, survey data contaminated with 7% Gaussian noise

4.3 不同尺寸滑動(dòng)窗口提取不同深度組合異常源的能力分析

在上述研究的基礎(chǔ)上，進(jìn)一步研究不同尺寸滑動(dòng)窗口提取不同深度組合異常源的能力，特設(shè)置如下地球物理模型III：由兩個(gè)剩余密度均為0.3 g/cm3且大小均為500 m×500 m×500 m的異常源組成，兩者質(zhì)心分別為(0 m，-1 500 m，500 m)與(0 m，1 500 m，700 m)。以大小50 m×50 m×50 m的物性網(wǎng)格離散該地球物理模型，網(wǎng)格沿著笛卡爾坐標(biāo)系三軸向的剖分?jǐn)?shù)nx、ny和nz分別為101、101和50；觀測(cè)點(diǎn)間距為50 m×50 m，觀測(cè)點(diǎn)個(gè)數(shù)為nx×ny=101×101=10 201，即勘探區(qū)尺寸為5 000 m×5 000 m。觀測(cè)點(diǎn)高度為地面上25 m。以步長(zhǎng)Δz沿Z軸向下移動(dòng)從而生成一系列的地球物理模型IV。

本小節(jié)中觀測(cè)數(shù)據(jù)的計(jì)算與4.1小節(jié)中相同，通過(guò)式(3)獲得組合異常源的歐拉解集后，利用模糊聚類對(duì)歐拉解集進(jìn)行分簇，通過(guò)式(8)計(jì)算歐拉解中的優(yōu)解率來(lái)判定最優(yōu)的滑動(dòng)窗口大小。這里優(yōu)解率定義為模糊聚類的簇中歐拉優(yōu)解的個(gè)數(shù)與歐拉解的個(gè)數(shù)的百分比。

圖6為組合異常源中深部異常源在不同尺寸滑動(dòng)窗口下的優(yōu)解率曲線，總的趨勢(shì)為隨著異常源深度的增加，優(yōu)解率在不斷地降低。對(duì)于深部異常源而言，隨著滑動(dòng)窗口的增大，優(yōu)解率先增大后逐漸減?。划?dāng)滑動(dòng)窗口尺寸約為異常源尺寸的3倍時(shí)優(yōu)解率最大，且當(dāng)滑動(dòng)窗口尺寸大于勘探區(qū)尺寸的1/2時(shí)，將導(dǎo)致歐拉反褶積難于提取到有用信息。

圖6 深部異常源不同尺寸滑動(dòng)窗口下的優(yōu)解率曲線Fig.6 Different moving window sizes for optimal Euler solution rate of deep anomalous sources with various offset Δz

圖7為組合異常源中淺部異常源在不同尺寸滑動(dòng)窗口下的優(yōu)解率曲線，總的趨勢(shì)與圖 6類似，隨著異常源深度的增加，優(yōu)解率在不斷地降低。當(dāng)Δz=100和Δz=200時(shí)，隨著滑動(dòng)窗口的增大，優(yōu)解率先增大后變緩，并因受到深部異常源的干擾，出現(xiàn)一定的波動(dòng)。當(dāng)Δz=300、Δz=400和Δz=500時(shí)，隨著滑動(dòng)窗口的增大，優(yōu)解率先增大后逐漸減??；當(dāng)滑動(dòng)窗口尺寸約為異常源尺寸的3倍時(shí)獲得最佳的優(yōu)解率。與圖6相比，由于深部異常屬于弱異常源，其位場(chǎng)總場(chǎng)效應(yīng)相對(duì)較小，因而滑動(dòng)窗口尺寸過(guò)大，易于導(dǎo)致歐拉反褶積難以提取有效異常。與圖4相比，對(duì)于淺部異常源而言，當(dāng)勘探區(qū)域內(nèi)僅為孤立異常源時(shí)，滑動(dòng)窗口尺寸相對(duì)較大有利于歐拉反褶積提取異常源信息；但與圖5相比，對(duì)于淺部異常源而言，滑動(dòng)窗口尺寸相對(duì)較大易于導(dǎo)致式(3)引入大量干擾信息(如觀測(cè)數(shù)據(jù)噪聲、相鄰弱異常源干擾等)。綜上，當(dāng)測(cè)量點(diǎn)距為異常源尺寸的1/10時(shí)，最優(yōu)化窗口尺寸應(yīng)為異常源尺寸的3倍且小于勘探區(qū)尺寸的1/2。

圖7 淺部異常源不同尺寸滑動(dòng)窗口下的優(yōu)解率曲線Fig.7 Different moving windows sizes for optimal Euler solution rate of shallow anomalous sources with various offset Δz

5 結(jié)論

(1)使用單一構(gòu)造指數(shù)過(guò)濾歐拉繆解，易于遺漏有用異常信息。采用估計(jì)函數(shù)εr過(guò)濾能有效的保留構(gòu)造指數(shù)相對(duì)較小的歐拉解。

(2)當(dāng)勘探區(qū)域大小固定時(shí)，隨著滑動(dòng)窗口的增大，有解率先增大后變小，當(dāng)滑動(dòng)窗口尺寸約為勘探區(qū)域尺寸的1/3～1/2時(shí)，有解率最大。

(3)對(duì)于孤立異常源而言，隨著滑動(dòng)窗口的不斷增大，歐拉反褶積的優(yōu)解率先穩(wěn)步增大后趨于穩(wěn)定；當(dāng)滑動(dòng)窗口很小時(shí)，淺部異常源對(duì)應(yīng)的優(yōu)解率占優(yōu)；隨著滑動(dòng)窗口的不斷增大，這一關(guān)系仍然得以保持；當(dāng)滑動(dòng)窗口尺寸大于勘探區(qū)域尺寸的1/2時(shí)，不同深度異常源對(duì)應(yīng)的優(yōu)解率趨于一致。

(4)對(duì)于有不同深度異常源的組合模型而言，隨著異常源深度的增加，優(yōu)解率在不斷地降低；由于深部異常屬于弱異常源，其位場(chǎng)總場(chǎng)效應(yīng)相對(duì)較小，因而滑動(dòng)窗口尺寸過(guò)大易于導(dǎo)致歐拉反褶積難以提取有效異常。對(duì)于淺部異常源而言，滑動(dòng)窗口尺寸相對(duì)較大易于引入大量干擾信息。因而在歐拉反褶積中，當(dāng)測(cè)量點(diǎn)距為異常源尺寸的1/10時(shí)，最優(yōu)化窗口的尺寸應(yīng)為異常源尺寸的3倍且小于勘探區(qū)尺寸的1/2。