基于WLs的二維地球化學(xué)景觀多重分形建模

高 歆

(中國(guó)地質(zhì)大學(xué)(北京)地球科學(xué)與資源學(xué)院， 北京 100083)

利用多重分形理論建模構(gòu)成了過去數(shù)十年來客觀自然現(xiàn)象比較流行的一種尺度分析手段，廣泛的被應(yīng)用在包括地學(xué)、湍流和因特網(wǎng)等自然科學(xué)領(lǐng)域中。為了有效的指導(dǎo)地球化學(xué)勘探和礦產(chǎn)資源評(píng)價(jià)，Mandelbrot和Cheng[1-2]提出了多種多重分形模型來描述地球化學(xué)數(shù)據(jù)的空間尺度結(jié)構(gòu)，這些新的模型不僅考慮地殼中礦石元素的分布幾何形狀，而且將傳統(tǒng)的地質(zhì)統(tǒng)計(jì)學(xué)納入到了計(jì)算分維譜的框架之下，其中尤其是數(shù)據(jù)域的奇異性分析，即計(jì)算空間中每點(diǎn)的奇異值，一種高通濾波方式，可以有效的從由各種地質(zhì)過程生成的混合模式和特征中分解背景值和異常值。

對(duì)多重分形譜的估計(jì)是多重分形分析的重要任務(wù)。針對(duì)多重分形譜建模，有兩種最常用的建模方法，其中一種是基于分維譜f(α)，另外一種是基于余維數(shù)c(γ)。兩種模型之間的內(nèi)在聯(lián)系現(xiàn)已得到證明[3]。多重分形理論已廣泛用于各種學(xué)科領(lǐng)域，并因此產(chǎn)生了各種各樣的計(jì)算多重分形譜函數(shù)的方法，如矩方法[4]、直方圖法、小波方法[5],以及二次維矩法[6]等，其中應(yīng)用最廣泛的是將空間分成許多格子的矩統(tǒng)計(jì)方法。另外，多重分形模型將服從正態(tài)或者對(duì)數(shù)正態(tài)分布的微量元素的背景值與服從分形分布(冪律分布)的高低異常值統(tǒng)一到多重分形譜中，多重分形分析體現(xiàn)了統(tǒng)計(jì)手段上的進(jìn)步,目前地學(xué)數(shù)據(jù)處理中，最常用到的是基于盒子的多重分析方法。Cheng等[7]通過頻率-濃度、數(shù)量-濃度、濃度-濃度、濃度-周長(zhǎng)、S-A等分形模型及矩分析等多重分形方法，將地球化學(xué)異常從背景中分離出來；郭科等[8]運(yùn)用多重分形模型，通過研究Pb、Zn、Sb、Ag、Cu共五種元素的廣義分形維數(shù)Dq，認(rèn)為這些元素在空間上的分布服從多重分形分布特征，同時(shí)指出多重分形方法相比于傳統(tǒng)的化探方法來說在判斷元素的共生組合方面具有一定的優(yōu)越性。Zheng等[9]運(yùn)用周長(zhǎng)-面積等分形學(xué)研究方法研究了中國(guó)3個(gè)金礦床中黃鐵礦、雄黃等礦物表面As、S、Fe、Si等元素在成礦過程中的遷移和轉(zhuǎn)換規(guī)律。

自從Parish和Frisch[10]提出用湍流域中的速度增量來刻畫尺度函數(shù)：Sp(l)=<(δfl)p>～lξp，δfl(x)=f(x+l)-f(x)是尺度域中的增量。不幸的是，Muzy[11]發(fā)現(xiàn)了該種方法只能刻畫信號(hào)中比較強(qiáng)的奇異性，對(duì)于一些弱的奇異性和譜值左半部分顯得無能為力。Mallat[12]提出了利用小波極大模值捕捉信號(hào)的奇異性以來，Arneodo等人提出了WTMM代替直接的速度增量來計(jì)算尺度函數(shù)，從而避免了上述缺點(diǎn)。然而，WTMM需要連續(xù)小波變換，以及利用不同的半徑在不同的尺度上進(jìn)行搜索鄰域中的模極大值，復(fù)雜的計(jì)算過程，高昂的計(jì)算量都使的這種方法變的不很實(shí)用。為了避免上述所述WTMM涉及到的繁瑣的工作量，法國(guó)數(shù)據(jù)學(xué)Jaffard[13]等人提出了WLs的計(jì)算方法，不僅能夠描述譜值的整個(gè)區(qū)域，還避免了連續(xù)小波變換以及尋值等繁瑣步驟。多次研究證明，計(jì)算分形譜的方法嚴(yán)重依賴方法和數(shù)據(jù)本身，WLs的引入為這種情況提供了一種補(bǔ)充。Dojnow[14]通過基于小波系數(shù)的多重分形模型，考察了含有噪聲信號(hào)的多條腦電圖曲線，指出分析信號(hào)的相位攜帶了最多的尺度不變性信息；Khalil等[15]利用基于WTMM的多重分形模型考察了間歇染色體結(jié)構(gòu)，將染色體版圖從背景中分離出來，提出了一種新的邊界分割方法，闡述了染色體結(jié)構(gòu)、位置和不穩(wěn)定性之間的潛在關(guān)系。

1 基于盒子的多重分形分析

(1)

這里的非整數(shù)αi稱為coarse holder指數(shù)[16]。把在分形上具有相同α值的小盒子數(shù)目記為Nα(ε)，它與ε的大小有關(guān)，并且可以寫成：

(2)

將上式與N(ε)～ε-D的簡(jiǎn)單分形公式對(duì)比，可見f(α)的物理意義是表示具有相同α值的子集的分形維數(shù)，或者說它描述了ε→0+直方圖N(ε)的變化，f(α)定義為具有相同α值的子集的Hausdorff維數(shù)，稱為多重分形譜或奇異性譜。在多重分形譜的眾多計(jì)算方法當(dāng)中，矩方法是最常用的方法之一。為了解f(α)的分布特性，定義分割函數(shù)：

(3)

可以獲得質(zhì)量指數(shù)函數(shù)τ(q)，τ(q)是關(guān)于q的凸函數(shù)。當(dāng)q=1時(shí)達(dá)到最大曲率，該點(diǎn)的曲率可以用于判斷是否具有多重分形性。通過勒讓德變換，可以得到α(q)和f(α)如下的兩個(gè)函數(shù)：

(4)

三個(gè)函數(shù)τ(q)、α(q)和f(α)是相互關(guān)聯(lián)的。

2 基于WLs的多重分形分析

(5)

繼續(xù)設(shè)H0(k)和G0(k)為兩個(gè)具有有限沖擊響應(yīng)的濾波器，引入二進(jìn)區(qū)間以及并集：

λj,k1,k2={[k12j,(k1+1)2j],[k22j,

(6)

(7)

理論上的WLs定義為(圖1)

(8)

圖1 WL的定義示意圖

WLs正是通過其局部奇異性刻畫能力構(gòu)成了了尺度分析模型，如果X(t)在t0點(diǎn)存在奇異性指數(shù)α(t0)≥0，消失矩Nψ≥α(t0)，WLs與尺度之間的冪律關(guān)系為：

(9)

通過計(jì)算位于尺度2j上的WLs的q階矩的平均值，可以得到結(jié)構(gòu)函數(shù)和質(zhì)量指數(shù)函數(shù)。

(10)

(11)

通過勒讓德變換，奇異值函數(shù)和分維譜函數(shù)定義為：

(12)

同樣，這三個(gè)函數(shù)τL(q)、αL(q)和fL(α)是相互聯(lián)系的。

3 方法驗(yàn)證

3.1 De Wijs模擬數(shù)據(jù)

為了對(duì)比這兩種方法的分析效果，本文采用De Wijs模擬數(shù)據(jù)[18]驗(yàn)證包括N=64和N=1028兩種不同尺度數(shù)據(jù)，生成算子為1.96、0.84、0.84、0.36。圖2是尺度指數(shù)和多重分形譜的估計(jì)值，對(duì)于兩種不同數(shù)據(jù)盒子估計(jì)法均獲得了比較的擬合，但是在譜值的首末端有點(diǎn)偏差，主要是由于邊緣效應(yīng)造成的。結(jié)果表明，對(duì)于小尺度的模擬數(shù)據(jù)，圖2A中與理論曲線相比，質(zhì)量指數(shù)曲線在左邊和右邊分別出現(xiàn)了上移和下移的情況，相應(yīng)的譜值左邊和右邊偏離標(biāo)準(zhǔn)的譜域，因?yàn)榛赪Ls的計(jì)算中對(duì)原小波系數(shù)進(jìn)行挑選導(dǎo)致的發(fā)散，而基于盒子法的估計(jì)比較接近理論值，也顯示了盒子法對(duì)于小尺度的測(cè)度數(shù)據(jù)的穩(wěn)定性和準(zhǔn)確性。同樣，值得注意的是，對(duì)于n=1的消失矩WLs估計(jì)，原理上等同于盒子法(盒子法：測(cè)度上的求和或者積分，一階消失矩：加權(quán)求和或者積分)，譜值的左邊現(xiàn)實(shí)了良好的逼近，即模擬數(shù)據(jù)中的小值模擬準(zhǔn)確，而對(duì)于大值出現(xiàn)了偏差。對(duì)于大尺度多重分形分析，盒子估計(jì)法和WLs的基本都和理論曲線實(shí)現(xiàn)了基本重疊，表明對(duì)于大尺度數(shù)據(jù)，WLs和盒子估計(jì)都可以用來多重分形分析。總之，對(duì)于小尺度的逼近，由于分辨率有限，邊緣效應(yīng)，以及高階消失矩的過分平滑都會(huì)導(dǎo)致估計(jì)發(fā)散，而對(duì)于大尺度的數(shù)據(jù)，兩者都可以作為分析手段，對(duì)于消失矩，可以去2、3、4階的均值來近似估計(jì)。

3.2 地球化學(xué)數(shù)據(jù)驗(yàn)證

本文地球化學(xué)數(shù)據(jù)驗(yàn)證來自云南個(gè)舊錫礦水系沉積物元素Sn、Cu和Pb(圖3)，研究區(qū)域?yàn)?9400 km2，，采樣間隔為2km，分辨率為2km×2km，N=76×50，分析小波為Daubechies，消失矩為Nψ=1,2,3,4。如圖4所示，對(duì)于Sn元素譜值估計(jì)，基于WLs不同的消失矩的計(jì)算，差別明顯，顯示了WLs估計(jì)得不穩(wěn)定性，而Cu和Pb的估計(jì)，在左邊還表示了部分的重疊，正則性方面強(qiáng)于Sn元素。另外，利用經(jīng)典級(jí)聯(lián)過程驗(yàn)證多重分形數(shù)據(jù)的正則性，也可獲得同樣的結(jié)果，本文限于篇幅，另文闡述。通過地球化學(xué)元素?cái)?shù)據(jù)的驗(yàn)證，盒子估計(jì)法展現(xiàn)出來了優(yōu)良的穩(wěn)定性和準(zhǔn)確性，對(duì)于小尺度數(shù)據(jù)WLs估計(jì)需慎重。從另外的角度來說，利用基于函數(shù)奇異性的這套框架某種程度上不符合以測(cè)度刻畫的地學(xué)數(shù)據(jù)，或許WLs優(yōu)勢(shì)在于震蕩的一維信號(hào)。

圖2 De Wijs模擬數(shù)據(jù)的質(zhì)量指數(shù)和分形譜(盒子法: “+”；WL:消失矩N=1,2,3,4，符號(hào)：“o”;“×”; “*”; “矩形”)

圖3 云南個(gè)舊水系沉積物Sn、Cu和Pb元素IDW插值圖(A-C: 分別為Sn、Cu和Pb元素濃度值)

圖4 地球化學(xué)元素?cái)?shù)據(jù)的質(zhì)量指數(shù)和分形譜(盒子法: “+”；WL:消失矩N=1,2,3,4，符號(hào)：“o”; “×”; “*”; “矩形”)

4 結(jié)論

對(duì)于分辨率有限、大小且不規(guī)則的現(xiàn)實(shí)地學(xué)數(shù)據(jù)，基于WLs的尺度分析，在準(zhǔn)確性和穩(wěn)定性方面都弱于盒子法，但是有時(shí)候尺度分析并不涉及到其多重分形譜的計(jì)算，僅僅利用系數(shù)與尺度之間的冪律關(guān)系，需要找到一個(gè)拐點(diǎn)，因此，小波系數(shù)不失為找不同尺度閾值一個(gè)好方法。對(duì)于小波消失矩的選擇，消失矩過小，捕捉奇異性的強(qiáng)度不夠，過大會(huì)導(dǎo)致分析的數(shù)據(jù)濾波過于平滑，從而導(dǎo)致譜值的兩端發(fā)散，因此，一般取消失矩的比較大的加權(quán)平均較合適。另外，基于WLs的估計(jì)涉及到小波變換的多分辨率分析，在圖像分割過程中，利用小波刻畫的奇異性進(jìn)行分割有著內(nèi)在的優(yōu)越性，延伸繼續(xù)刻畫多重分形譜，正如前文引用，是一個(gè)有效的方法。

這篇文章主要討論了以積分為主的測(cè)度和以累積量為主的尺度分析，以經(jīng)典礦物分割模型de Wijs模型模擬數(shù)據(jù)為例，結(jié)果表明，對(duì)于盒子法的尺度分析，不管是小尺度還是大尺度，都表現(xiàn)出了更優(yōu)的結(jié)果，基于WLs的尺度分析對(duì)于大尺度的數(shù)據(jù)也表現(xiàn)出了良好的結(jié)果。此外，基于小波系數(shù)的WLs不僅可以為分析數(shù)據(jù)提供多分辨率分析，還可以檢測(cè)信號(hào)的振蕩器異性等行為。總之，一種是基于概率測(cè)度下的奇異性，一種基于連續(xù)函數(shù)的連續(xù)性，分別來源于不同的數(shù)據(jù)應(yīng)用，對(duì)于大分辨率的數(shù)據(jù)，兩者的刻畫能力相同，而對(duì)于小尺度的地學(xué)數(shù)據(jù)，盒子法的效果更優(yōu)。

[1] Mandelbrot BB. Intermittent turbulence in self similar cascades: Divergence of high moments and dimension of the carrier [J]. Fluid. Mech., 1974, 62: 331.

[2] Cheng Q. Multifractality and spatial statistics [J]. Computers & Geosciences, 1999, 25(9): 949-961.

[3] Cheng Q, Agterberg FP. Comparison between two types of multifractal modeling [J]. Mathmatical Geoglogy, 1996, 28(8):1001-1016.

[4] Halsey TC, Jensen MH, Kadanoff LP. Fractal measures and their singularities: The characterization of strange sets [J]. Phys. Rev., 1986, 33(2): 1141-1151

[5] Arneodo A, Bacry E, Muzy JF. The thermodynamics of fractals revisited with wavelets [J]. Physica A, 1995, 213: 232-275.

[6] Schertzer D, Lovejoy S, eds. Non-linear variability in geophysics [M]. The Netherlands: Kluwcr Acad. Publ., Dordrecht, 1991.

[7] Parisi G., Frisch U. On the singularity structure of fully developed turbulence, in Turbulence and Predictability in Geophysical Fluid Dynamics [M], edited by M. Ghil, R. Benzi, and G. Parisi, New York: Elsevier,1985.

[8] Muzy JF, Bacry E, Arneodo A.. Wavelets and multifractal formalism for singular signals-application to turbulent data [J]. Physical Review Letters, 1991, 67(25):3515-3518.

[9] Mallat S, Hwang WL. Singularity detection and processing with wavelets [J]. IEEE Trans. Inf. Theory, 1992, 38, 617-643.

[10] Jaffard S. Multifractal formalism for functions [J]. SIAM J. Math. Anal., 1997, 28(4): 944-998.

[11] Cheng QM, Agterberg FP, Ballantyne SB. the separation of geochemical anomalies from background by fractal methods [J]. Journal of geochemical exploration, 1994, 51(2):109-130.

[12] 郭科，施澤興，唐菊興等.用多重分形研究元素的共生組合.電子科技大學(xué)學(xué)報(bào),2004,33(2)：221-224.

[13] Zheng Z, Mao H, Cheng QM. Fractal geometry of element distribution on mineral surfaces [J]. mathematical of geology, 2001, 33(2): 217-228.

[14] Dojnow P. Multifractal analysis of amplitude and phase components of EEG-based analytical signals [J]. Comptes Rendus De L Academic Bulgare Des Sciences, 2007, 60(10):1071-1076.

[15] Khalil A, Grant GL, Caddle Lb, et al. Chromosome territories have a highly nonspherical morphology and nonrandom positioning [J]. Chromosome Research, 2007, 15(7):899-916.

[16] Evertsz CJG., Mandelbrot BB. Multifractal measures, In: Peitgen, H.O., Jurgens, H., Saupe, D. (Eds.), Chaos and Fractals. Springer [M], New York: Springer, 1992.

[17] Wendt H, Roux SG., Jaffard S, et al. Wavelet leaders and bootstrap for multifractal analysis of images [J]. Signal Processing, 2008, 89(6): 1100-1114.

[18] De Wijs, HJ, Statistics of ore distribution [J], Geologic en Mijnbouw. 1951, 13 (8), 365-375.