，,3，

(1.中國地質(zhì)大學(xué)地質(zhì)過程與礦產(chǎn)資源國家重點(diǎn)實(shí)驗(yàn)室，湖北武漢430074； 2.中國地質(zhì)大學(xué)(武漢)資源學(xué)院，湖北武漢430074； 3.York大學(xué)地球空間科學(xué)系和地理系，加拿大多倫多M3J1P3)

基于柵格數(shù)據(jù)的局部奇異性分析迭代方法及其軟件實(shí)現(xiàn)

陳志軍1,2，成秋明1,2,3，陳建國1,2

(1.中國地質(zhì)大學(xué)地質(zhì)過程與礦產(chǎn)資源國家重點(diǎn)實(shí)驗(yàn)室，湖北武漢430074； 2.中國地質(zhì)大學(xué)(武漢)資源學(xué)院，湖北武漢430074； 3.York大學(xué)地球空間科學(xué)系和地理系，加拿大多倫多M3J1P3)

近些年，弱緩化探異常識別已成為成礦預(yù)測和勘查評價(jià)中十分關(guān)鍵的問題。多重分形理論的局部奇異性分析因其嶄新的思路、簡便的方法、優(yōu)良的效果而在弱緩異常識別中受到廣泛關(guān)注。在深入剖析局部奇異性分析基本思想及計(jì)算方法的基礎(chǔ)上，引入正規(guī)化尺度參數(shù)L，提出了迭代方法來改進(jìn)局部奇異性指數(shù)的估值。給出了奇異性迭代算法并用C++編程實(shí)現(xiàn)，軟件功能強(qiáng)勁，操作靈活，不僅適用于各向同性情形，還適用于各向異性尺度的奇異性計(jì)算和方向性奇異性計(jì)算。軟件的動態(tài)鏈接庫版本已在GeoDAS礦產(chǎn)資源定量預(yù)測專業(yè)軟件中應(yīng)用。

局部奇異性分析; 柵格數(shù)據(jù); 多重分形; 迭代方法; 正規(guī)化尺度參數(shù)

0 引 言

分形與多重分形在勘查地球化學(xué)數(shù)據(jù)、勘查地球物理數(shù)據(jù)、礦石品位數(shù)據(jù)等地質(zhì)勘查數(shù)據(jù)中的普遍存在性為成礦過程奇異性與礦產(chǎn)預(yù)測定量化的新理論與新方法研究提供了嶄新的思路和新型的數(shù)學(xué)工具(成秋明, 2007; Agterberg, 2003)。不均勻性、自相似性(自放射性、廣義自相似性)、奇異性等通常被用來描述地球化學(xué)場特征、礦床空間分布、熱液成礦系統(tǒng)礦床品位-噸位分形模型的非線性性質(zhì)。一種五參數(shù)的多重級聯(lián)模型最近被提出，該模型對de Wijs’s鋅數(shù)據(jù)等多重分形經(jīng)典數(shù)據(jù)的非對稱奇異譜進(jìn)行了成功建模，為現(xiàn)實(shí)世界中地球化學(xué)場多重分形性的存在性提供了理論依據(jù)。分形維數(shù)、多重分形對稱度、多重分形奇異譜等參數(shù)從非線性角度對描述成礦物質(zhì)富集與貧化的統(tǒng)計(jì)規(guī)律發(fā)揮了重要作用(成秋明, 2000, 2007)，然而，對于礦產(chǎn)資源定量預(yù)測中的“定位”(礦產(chǎn)資源空間分布位置)預(yù)測功能則有所不足，需要將全局性統(tǒng)計(jì)轉(zhuǎn)變?yōu)榫植炕Ｐ??；诙嘀胤中卫碚摻⒌木植科娈愋苑治瞿Ｐ?，試圖通過局部異常的有效識別來達(dá)到空間預(yù)測目的，該模型正逐漸在覆蓋區(qū)勘查地球化學(xué)弱緩異常識別中表現(xiàn)出巨大的應(yīng)用潛力，成為非線性礦產(chǎn)預(yù)測的有效方法之一(成秋明, 2012)。應(yīng)用的需求也促進(jìn)了局部奇異性分析自身模型的發(fā)展，出現(xiàn)了多種形式的擴(kuò)展模型(陳志軍, 2007)。主要從提高奇異性指數(shù)計(jì)算精度方面介紹局部奇異性分析迭代分析方法，設(shè)計(jì)算法并編程實(shí)現(xiàn)，同時(shí)指明準(zhǔn)確計(jì)算奇異性指數(shù)需要注意的一些關(guān)鍵參數(shù)設(shè)置。

1 局部奇異性分析方法的基本思想

奇異性理論指出：奇異性是指在很小的時(shí)間-空間范圍內(nèi)具有巨大能量釋放或巨量物質(zhì)形成的現(xiàn)象(成秋明, 2006; Cheng et al, 2009)。假定在E維(拓?fù)?空間中有一由質(zhì)點(diǎn)凝聚成分維數(shù)為D的分形(如礦物的富集)，L尺度上質(zhì)點(diǎn)集合的質(zhì)量M(L)與集合中質(zhì)點(diǎn)的數(shù)目N(L)成正比，則該分形體的密度可表示為(陳颙等, 2005)：

(1)

由式(1)可知，當(dāng)D

該模型的基本內(nèi)容：記B(x,ε)為任意空間位置x，尺度大小為ε的盒子(可以是方形、矩形、或其他更復(fù)雜的形狀)，在盒子B(x,ε)覆蓋區(qū)域內(nèi)的金屬量為μ(B(x,ε))，則有：

μ(B(x,ε))=c(x,ε)εα(x)

(2)

為便于理解與制圖，令：

Δα(x) =E-α(x)

(3)

設(shè)ε大小鄰域上的平均密度為〈ρ(B(x,ε))〉，〈〉表示統(tǒng)計(jì)期望，式(2)于是寫成密度形式，有：

〈ρ(B(x,ε))〉 =μ(B(x,ε)) /εE=c(x,ε)ε-Δα(x)

(4)

式(4)中，E為研究對象的空間拓?fù)渚S數(shù)(E=1,2,3，分別對應(yīng)一維、二維或三維問題)，待定系數(shù)α(x)稱為空間位置x處的局部奇異性指數(shù)，它表征了奇異性的強(qiáng)弱程度，是與尺度無關(guān)的量，無量綱；若ε趨近于0且c(x,ε)存在，則此時(shí)的c值被認(rèn)為是位置x處的“α(x)維分形密度”。若ε用最大尺度進(jìn)行了歸一化，則c(x,ε)值是對該最大尺度的滑動平均的逼近。對地球化學(xué)采樣數(shù)據(jù)，測量所得數(shù)據(jù)通常表示的是濃度含量，可視為一種廣義的“密度”數(shù)據(jù)。因此，利用式(4)可直接對化探測量原始數(shù)據(jù)建模，對結(jié)果解釋更直觀，故而在分析中式(4)比式(3)更常用。需要注意的是，式(3)也是非常有用的，特別是在奇異性指數(shù)的算法設(shè)計(jì)方面。式(3)在計(jì)算中可直接獲得α值并且其擬合優(yōu)度高于用式(4)直接獲得的Δα。這是由于普通最小二乘(OLS)模型只考慮縱向誤差，而α的均值水平趨近于E、Δα的均值水平常趨近于0的緣故。(-Δα)作為log-log坐標(biāo)系中回歸直線的斜率，若其結(jié)果近似為0，意味著自變量和應(yīng)變量不存在相關(guān)關(guān)系；同樣的數(shù)據(jù)用α來表示(此時(shí)近似為E)，自變量和應(yīng)變量的相關(guān)關(guān)系將變得顯著。無論用式(3)還是式(4)計(jì)算，獲得的奇異性指數(shù)和局部系數(shù)是一致的，不同的是log-log坐標(biāo)系下對應(yīng)的回歸方程具有不同的顯著性水平。要解決此問題，對式(4)和式(3)在log-log坐標(biāo)系下可選用標(biāo)準(zhǔn)加權(quán)最小二乘法(SWLS)，權(quán)重為1/(1+b2)(其中b為擬合直線斜率)，以觀測點(diǎn)到回歸直線的垂直距離平方和最小化來確定直線參數(shù)，當(dāng)然該計(jì)算過程比普通最小二乘模型略顯復(fù)雜。

奇異性指數(shù)α值(或Δα)有著優(yōu)良的異常指示功能，在實(shí)際應(yīng)用中常用來度量不同空間位置物質(zhì)(或能量)較其周圍相對富集或虧損的程度。α值與E偏差越大，說明局部富集或虧損越越強(qiáng)烈，也即奇異性越強(qiáng)。以二維化探異常研究為例，此時(shí)E取值為2，當(dāng)α<2(即Δα>0)表示局部正異常，元素相對局部富集；α>2(即Δα<0)表示局部負(fù)異常，元素相對局部虧損；而α= 2(即Δα=0)表示無異常，屬于背景。通常在全圖范圍內(nèi)計(jì)算各個(gè)位置的α值，對α值(或Δα值)進(jìn)行制圖表達(dá)，從而來識別不同空間位置的局部化探異常，并進(jìn)一步篩選出有用異常并在空間上圈定出來。沒有奇異性的區(qū)域(α=E)所構(gòu)成子集的分形維數(shù)接近于分形盒維數(shù)E，對于地球化學(xué)背景區(qū)域，它們的空間分布通常占全圖的絕大部分，其統(tǒng)計(jì)分布規(guī)律通常符合正態(tài)分布或?qū)?shù)正態(tài)分布；而奇異性區(qū)域(α≠E)所構(gòu)成子集的分形維數(shù)有很多個(gè)，可以由f(α)

利用奇異性指數(shù)識別化探異常的效果經(jīng)常被用來與襯值異常進(jìn)行比較。陳志軍等(2009)探討了α值和滑動襯值(簡記為CV)之間的定量關(guān)系。當(dāng)無標(biāo)度區(qū)間上符合嚴(yán)格的冪律關(guān)系時(shí)，α的估值可任選2個(gè)尺度進(jìn)行計(jì)算，有如下關(guān)系成立：

(5)

式(5)中，la和lb分別表示較小和較大2個(gè)不同尺度的窗口；Δα和log(CV)相比較而非直接比較CV，在統(tǒng)計(jì)上更能區(qū)分辨別二者的異常能力。由于對數(shù)函數(shù)是典型的單調(diào)函數(shù)，log(CV)和CV具有相同的排序值，可以通過比較Δα和CV的排序值(或用百分位數(shù))來判別異常識別能力(陳志軍, 2007)。

局部奇異性指數(shù)的計(jì)算通常要考察多個(gè)尺度的趨勢性變化，由統(tǒng)計(jì)分形關(guān)系來確定，在log-log雙對數(shù)圖中采用回歸分析技術(shù)來估值，相對僅由2個(gè)窗口度量值所確定的滑動襯值來說，Δα值對局部異常強(qiáng)度的指示作用更穩(wěn)健，更能代表局部范圍內(nèi)異常的總體水平。滑動襯值則常因2個(gè)窗口選擇的隨意性而易受波動。通常，局部奇異性指數(shù)相比滑動襯值異常識別能力更佳(陳志軍等, 2009)。

2 基于柵格數(shù)據(jù)的算法剖析

2.1 柵格數(shù)據(jù)模型

從計(jì)算簡便性考慮，二維情形下基于柵格數(shù)據(jù)模型的計(jì)算方案效率最高。在柵格形式表達(dá)中，雖然對于柵格的形狀沒有限制(如等三角形、等六邊形等)，然而正方形網(wǎng)格的簡單性、運(yùn)算方便性、直觀性使得這種形式的柵格數(shù)據(jù)模型最為普遍。這些形狀一樣的正方形網(wǎng)格作為基本單元(單元格或稱像元)組成了形如數(shù)學(xué)中規(guī)則排列的矩陣。像元大小、空間位置及像元值構(gòu)成了柵格數(shù)據(jù)的基本要素。像元大小表示數(shù)據(jù)的最小尺度，或稱空間分辨率；像元的取值即為空間對象的測量值；每個(gè)像元有唯一的行索引和列索引，測量值的空間位置隱含于柵格行列位置，也即可根據(jù)行列號轉(zhuǎn)換成相應(yīng)的空間坐標(biāo)。一個(gè)柵格數(shù)據(jù)集描述了某個(gè)專題內(nèi)容在區(qū)域的位置、特性和空間上的相對位置。

眾多GIS軟件都支持柵格數(shù)據(jù)類型，并且有眾多數(shù)據(jù)格式。表1展示了ArcGIS和Surfer 2種主流的柵格數(shù)據(jù)文本文件格式，由此可以直觀地了解柵格數(shù)據(jù)的編碼規(guī)則。在對這兩類柵格數(shù)據(jù)文件進(jìn)行相互轉(zhuǎn)化中，需要注意ArcGIS ASC文件所定義的左下角位置是左下角柵格單元的左下邊界的角點(diǎn)，而Surfer 6 Text Grid柵格數(shù)據(jù)所定義的X最小值和Y最小值位于柵格單元的中心位置，因此兩者偏差半個(gè)像元大小(表1示例文件中該偏差值為50)。ArcGIS ASC文件文件頭信息可自定義缺失數(shù)據(jù)標(biāo)記值，如-99 999；而在Surfer 6 Text Grid文件中，缺失數(shù)據(jù)默認(rèn)標(biāo)記為1.701 41E+038。此外，兩者數(shù)據(jù)陣列的存放方案是不一樣的，具有上下對稱的關(guān)系。

表1 ArcGIS和Surfer柵格數(shù)據(jù)文件格式示例

2.2 局部奇異性分析算法中的量綱問題

局部奇異性分析算法通?；跂鸥駭?shù)據(jù)的多窗口鄰域統(tǒng)計(jì)來開展計(jì)算。利用柵格數(shù)據(jù)在鄰域分析方面的簡便性，可以快速獲得各個(gè)空間位置上不同尺度的滑動平均值，進(jìn)而在log-log坐標(biāo)下進(jìn)一步計(jì)算奇異性指數(shù)(空間維數(shù)E=2)。柵格數(shù)據(jù)模型具有統(tǒng)一的幾何支撐，意味著像元分辨率范圍內(nèi)研究對象具有同質(zhì)性。對于某個(gè)柵格數(shù)據(jù)來說，空間分辨率一旦確定，局部奇異性特征尺度區(qū)間的最小值也就隨之確定。因此，獲取的局部奇異性指數(shù)是在一定尺度區(qū)間上估計(jì)出來的，揭示的是統(tǒng)計(jì)分形規(guī)律，在現(xiàn)實(shí)中特征尺度不能無限趨近于0。

局部奇異性分析方法算法的一些關(guān)鍵計(jì)算技術(shù)如：空間覆蓋系列盒子的各向同性、各向異性窗口的多種構(gòu)造方式、空間加權(quán)處理、掩模矩陣構(gòu)置、等效尺度的確定及l(fā)og-log最小二乘線性擬合、數(shù)據(jù)邊緣區(qū)擴(kuò)邊、多種測度參見文獻(xiàn)(陳志軍，2007)。為便于推導(dǎo)迭代算法，這里有必要分析一下尺度ε的處理問題。ε的量綱是否及如何影響計(jì)算結(jié)果？假定窗口由小到大依次遞增序列為：lmin=l1

2.2.1ε=l具有長度量綱 此時(shí)，c值在理論上的量綱是[ρ][l]Δα，估計(jì)Δα值(或α值)和c值的回歸方程可表示為：

log〈ρ(x,li)〉 = log(c(x))-Δα(x)×log(li),i=1,2,3，…,m

(6)

用α(x)表示，也即

log〈ρ(x,li)×li2〉=log(c(x))+α(x)×log(li)，i=1,2,3，…,m

(7)

其中，Δα(x)和log(c(x))為待定系數(shù)。式(6)與(7)相比，式(6)在形式上更為簡單，但是由于Δα(x)均值水平接近于0而使得對普通最小二乘(OLS)準(zhǔn)則下回歸方程顯著性水平偏低，作者建議先用式(7)形式計(jì)算α(x)，再由式(3)轉(zhuǎn)換成Δα(x)。此時(shí)，回歸直線段斜率即代表α(x)，通常與柵格數(shù)據(jù)模型的維數(shù)2比較接近，然后由式(4)計(jì)算Δα(x)。當(dāng)li=1(帶長度單位)時(shí)，可以獲得擬合值即為c(x)值。這里，需要注意到c(x,l1=1)、c(x,l1=lmin)、c(x,l1=lcellsize)三者可以是完全不同的，c(x,l1=lcellsize)才表示原始柵格分辨率意義下對原始測量值的估值。

2.2.2ε=l/L=N(l)/N(L) 消除量綱 引入正規(guī)化尺度參數(shù)L(為一定值)，通過比值l/L消除ε的量綱。在柵格數(shù)據(jù)模型中，可以直接用窗口邊長所占像元個(gè)數(shù)的比值：N(l)/N(L)等效表達(dá)物理長度比值l/L，這里，N(·)操作符表示取長度覆蓋的柵格像元數(shù)，允許取小數(shù)表示等效個(gè)數(shù)。對非方形窗口，借鑒面積與邊長的關(guān)系，N(l)的大小可由窗口覆蓋像元數(shù)的平方根來等效表達(dá)。此時(shí)，Δα(x)值與c值均無量綱，Δα(x)值大小與不消除量綱情形獲得的結(jié)果相同，而c值則不同。式(3)可轉(zhuǎn)變?yōu)椋?/p>

〈ρ(x,l)〉=c(x,l)(l/L)-Δα(x)

(8)

對x位置處奇異性指數(shù)及其系數(shù)的估計(jì)可由log(l/L)-log(ρ(x,l)×(l/L)2)圖中回歸直線來確定，對應(yīng)回歸方程可表示為：

log〈ρ(x,li)×(li/L)2)〉=log(c(x))+α(x)×log(li/L)，i=1,2,3，…,m

(9)

L參數(shù)加入到模型中，增加了模型處理獲得c值的靈活性，以下介紹的迭代方法與L參數(shù)有緊密聯(lián)系。L參數(shù)的取值可以是多樣化的，下面列舉4種典型的情形。

(1)L=lcellsize(柵格數(shù)據(jù)空間分辨率)，此時(shí)，ε的大小就是柵格像元個(gè)數(shù)，也即ε=l/L=N(l)。對柵格數(shù)據(jù)計(jì)算窗口邊長通過數(shù)柵格像元個(gè)數(shù)在實(shí)際中操作簡單，僅需對計(jì)算窗口的行列號進(jìn)行操作，不涉及具體空間距離計(jì)算，因此在奇異性指數(shù)計(jì)算中使用頻繁，此種情形c(x)值為原始數(shù)據(jù)的非線性擬合值。

(2)L=lmax(尺度窗口序列的最大值)，這種取值方式在許多分形計(jì)算場合常被采用。log(c(x))值在式(9)中為log-log圖中截距(此時(shí)li=L)，相當(dāng)于在lmax×lmax窗口滑動平均對數(shù)值的擬合值。在log-log擬合中，當(dāng)L取lmax時(shí)，log(c(x))值在橫軸的最大值位置獲得，從回歸分析誤差角度考慮，由于偏離擬合數(shù)據(jù)的重心而有較大的誤差。

(3)L取值為尺度窗口序列{li} (i=1,2,3，…,m)的幾何平均值：

或表示成等效窗口數(shù)：

(9)

此種情形克服了式(2)中的不足，log-log回歸分析中平移縱軸使之經(jīng)過數(shù)據(jù)重心，從而使得c(x)估值的誤差最小化。用式(9)確定的L進(jìn)行奇異性分析時(shí)，用c(x)值擬合相應(yīng)窗口大小的滑動平均值具有最佳的預(yù)測可信度。

(4)L取值小于柵格數(shù)據(jù)的分辨率lcellsize。例如L=lcellsize/ 3，此時(shí)c(x)估計(jì)值表示了x位置處(lcellsize/ 3)空間分辨率意義下的分形插值結(jié)果。

綜上可見，ε的量綱對于Δα值(或α值)的估算沒有影響，而c值的大小及其量綱有差異，它的大小取決于正規(guī)化尺度參數(shù)L的取值方式。某些作者在開展局部奇異性分析模型時(shí)，對尺度指標(biāo)ε的取值張冠李戴，不加區(qū)分：在實(shí)際處理中采用了像元個(gè)數(shù)N(l)來進(jìn)行計(jì)算，卻聲稱窗口計(jì)算參數(shù)為空間距離某某千米變化到某某千米。盡管對于α(x)的估值不受影響，但嚴(yán)格來說這是不妥當(dāng)?shù)?，忽視了c(x)值的數(shù)量水平、量綱及其應(yīng)有的意義。

從模型計(jì)算的通用性角度考慮，筆者建議用形式上更為一般化的式(8)作為基于柵格數(shù)據(jù)模型的常規(guī)局部奇異性的算法?？梢栽谶B續(xù)區(qū)間[lmin,lmax]上內(nèi)插預(yù)測及適度外推預(yù)測任意窗口大小的滑動平均結(jié)果，或者獲得更小分辨率的分形插值結(jié)果。本研究所介紹的迭代方法是在L>lcellsize的情形下作出推導(dǎo)的，此時(shí)c(x)值與原始數(shù)據(jù)具有相同的量綱，且c(x)值代表了L尺度上的光滑成分。

3 迭代方法及算法設(shè)計(jì)

3.1 迭代方法的基本思想

局部奇異性分析方法的特點(diǎn)在于采用局部化的α估值與空間拓?fù)渚S數(shù)E值的差異性來實(shí)現(xiàn)異常識別，于是合理估計(jì)α值(或Δα值)就成為一個(gè)關(guān)鍵問題。利用一系列方形窗口基于柵格模型的局部奇異性分析算法簡單，運(yùn)算快捷。然而，嵌套的這一系列尺度窗口一次性計(jì)算所得的α值是最好的嗎？最優(yōu)的α值應(yīng)該符合怎樣的數(shù)學(xué)條件？筆者注意到與α值所伴隨的c值的特性有其重要作用。從式(4)可見，局部奇異性分析可將原始數(shù)據(jù)分離出兩大因子Δα成分和c成分。如果各個(gè)位置處的Δα值完全刻畫了奇異性的強(qiáng)度，那么由全體c值構(gòu)成的c集就應(yīng)成為一個(gè)非奇異性的成分；否則，若c集仍帶有奇異性，那么對應(yīng)的Δα集便與理想的奇異性指數(shù)不對應(yīng)，c集還可以繼續(xù)分離出新的Δα成分和c成分。Δα成分和c成分應(yīng)該分別代表奇異性因子和非奇異性因子，c成分應(yīng)足夠光滑，沒有奇異性，具有良好的連續(xù)性和可微性。在實(shí)際計(jì)算中，獲取Δα成分的log-log擬合過程中總會存在剩余誤差，它并沒有完全保證c成分不再帶有奇異性信息。設(shè)法去除c集中的奇異性成分，使之具有足夠高的光滑性，而將奇異性信息都?xì)w屬于Δα集，這成為筆者研究局部奇異性分析迭代方法來優(yōu)化α估值的關(guān)鍵思路(Chen et al, 2005, 2007)。

3.2 迭代方法的數(shù)學(xué)表達(dá)

記α*(x)和c*(x)分別為期望獲取的一種最優(yōu)的局部奇異性指數(shù)和局部系數(shù)，仍有

ρ(x)=c*(x)εα*(x)-E=c*(x)ε-Δa*(x)

(10)

成立。這里，ε=l/L且L>lcellsize，盒子B(x,ε) 內(nèi)的平均密度簡記為ρ(x)，ρ(x)在所有空間位置的取值構(gòu)成的數(shù)據(jù)集合稱為ρ集，c*(x)在所有空間位置的取值構(gòu)成的數(shù)據(jù)集合稱為c*集。注意到，c*集與ρ集量綱相同，c*集反映了比ρ集更大尺度上的滑動平均非線性擬合信息，c*集比ρ集具有更小的離散程度。c*集可看作是ρ集的一個(gè)粗尺度逼近，其細(xì)節(jié)(粗糙度)信息保存在Δα*圖中，理論上的c*集應(yīng)不再具有多重分形特征，這樣ρ集的奇異性信息都應(yīng)集中反映在α*集中。若對c*集再做奇異性分析，則其計(jì)算所得的各個(gè)α值都趨近于同一個(gè)常數(shù)，即空間拓?fù)渚S數(shù)E。

基于同尺度的局部系數(shù)圖迭代方法可以不斷校正奇異性指數(shù)，將具有奇異性的c(x)通過迭代不斷演化成非奇異的(Chen et al, 2007)：首先將ρ(x)分解成c0(x)和α0(x)，被估計(jì)出來的c0(x)可以進(jìn)一步估計(jì)其局部奇異性，即把它當(dāng)作新的ρ(x)，再次進(jìn)行分解得到新的c1(x)，(下標(biāo)1表示進(jìn)行第1次迭代數(shù)，下同)，可知c1(x)的離散程度比c0(x)小，比ρ(x)更小，不斷重復(fù)這樣的過程，可以預(yù)料，最終獲取的局部系數(shù)圖將趨近于一種平穩(wěn)狀態(tài)，最終將不再具有奇異性，此時(shí)它的局部奇異性指數(shù)α趨近于自身的歐氏空間維數(shù)E，即Δα趨近于0。迭代方法與常規(guī)的非迭代方法的區(qū)別在于前者考慮了局部系數(shù)c(x)的作用。

令c-1(x) =ρ(x)，于是有如下迭代形式：

ck-1(x)=ck(x)ε-Δαk(x),k=0,1,2,...

(11)

c集不斷視為新的輸入數(shù)據(jù)集，再次進(jìn)行同樣的局部奇異性計(jì)算，如此循環(huán)，新產(chǎn)生的c集將變得越來越光滑，趨近于一種平穩(wěn)狀態(tài)……于是最終可獲得所期望的c*集。迭代分解的示意見圖1。

圖1 局部奇異性分析迭代方法示意圖

假定進(jìn)行了n次迭代時(shí)達(dá)到理論上的迭代終止條件：αn→0時(shí)，研究獲取包括第0次(非迭代)和n次迭代的結(jié)果，有式(12)成立。

對比式(10)和式(12)，于是可得最終局部奇異性指數(shù)及系數(shù)的估值，見式(13)。

(13)

從中可見，當(dāng)n=0時(shí)，迭代方法退化為常規(guī)方法；n>1時(shí)，迭代方法獲取的最終奇異性指數(shù)α*的結(jié)果中則多了一個(gè)附加校正項(xiàng)，Δα*的表達(dá)則更為簡練，是歷次迭代所產(chǎn)生的Δαk的累加?？梢?，迭代方法具有對常規(guī)方法結(jié)果進(jìn)行校正的能力，且其形式簡潔。奇異性迭代方法得到了學(xué)者的關(guān)注，α*被稱為最終奇異性指數(shù)(Agterberg, 2012)。

3.3 迭代終止條件與邊緣效應(yīng)

理論上的迭代終止條件要求Δαn中處處取值為0，在實(shí)際計(jì)算中這一條件過于嚴(yán)格，可用均方根誤差(RMSE)作為判別條件，而均方根誤差在數(shù)值上與Δαn的標(biāo)準(zhǔn)方差等價(jià)，于是，可直接用s(Δαk) < eps(eps為精度要求，例如0.01)迭代約束條件。

當(dāng)精度要求過高時(shí)，將導(dǎo)致迭代過程無法在期望時(shí)間內(nèi)結(jié)束，因此與一般的迭代算法一樣，可人為限制迭代最多次數(shù)N0，使之定能在有限時(shí)間內(nèi)完成計(jì)算。當(dāng)?shù)?jì)數(shù)器k>N0時(shí)，跳出循環(huán)，終止迭代。作者對地質(zhì)勘查數(shù)據(jù)的大量計(jì)算表明，經(jīng)過數(shù)次迭代(通常3～5次)就能大幅度地降低c集的不光滑度，獲得較理想的奇異性指數(shù)。當(dāng)N0取值為0時(shí)，迭代方法就退化為常規(guī)方法。

當(dāng)多尺度窗口滑動到柵格數(shù)據(jù)邊緣位置時(shí)，一些窗口部分位置因覆蓋于數(shù)據(jù)范圍之外而造成數(shù)據(jù)獲取的不完整，從而產(chǎn)生邊緣效應(yīng)，窗口尺度越大且迭代次數(shù)越多，邊緣效應(yīng)會越來越嚴(yán)重。對此問題，可以采用傅里葉變換類似的手段來擴(kuò)邊，補(bǔ)齊數(shù)據(jù)范圍之外的像元取值；還可以采用缺失數(shù)據(jù)處理手段，僅對有效數(shù)據(jù)作鄰域統(tǒng)計(jì)量計(jì)算，獲得等效的平均密度。

3.4 基于柵格數(shù)據(jù)的迭代算法設(shè)計(jì)

局部奇異性分析的迭代算法可設(shè)計(jì)如下：

設(shè)置多尺度窗口參數(shù){li},i=1,2,3，…,m

設(shè)置L參數(shù)且L>lcellsize

設(shè)置迭代精度eps與最多次數(shù)N0

輸入數(shù)據(jù)初始化為原始柵格數(shù)據(jù)

最終的局部奇異性指數(shù)Δα*所有元素初始化為0

for(intk= 0,i<=N0;k++) { //約束迭代次數(shù)不超過N0

for(inti= 0,i

for(int j=0;j

計(jì)算[i][j]柵格焦元位置m個(gè)不同

大小窗口的鄰域平均值 //**

雙對數(shù)log-log回歸分析獲得奇異性

指數(shù)αk[i][j]和系數(shù)ck[i][j]

}

}

Δαk=E-αk; //默認(rèn)二維情形E=2；退化為方向性奇異性，則E=1

Δα*+=Δαk; //最終局部奇異性指數(shù)的累計(jì)計(jì)算

if(s(Δαk)< eps) break; //迭代達(dá)到精度要求，結(jié)束迭代

else輸入數(shù)據(jù)賦值更新為計(jì)算所得的ck集

}

c*=ck; //獲得最終的局部系數(shù)

注：**步驟可以采用積分圖快速算法。

4 軟件功能實(shí)現(xiàn)及其效果

利用C++編程語言在Windows平臺下開發(fā)了RasterDataMining 1.0軟件，實(shí)現(xiàn)了基于柵格數(shù)據(jù)模型的局部奇異性分析專業(yè)分析功能。以主流文本柵格格式ArcGIS ASC Grid、Surfer 6 Text Grid來讀寫文件，制圖功能則可在ArcGIS或Surfer平臺下完成。該軟件從柵格數(shù)據(jù)局部奇異性分析的深度應(yīng)用需求出發(fā)，具有如下功能(圖2)。

圖2 局部奇異性程序功能界面(圖上帶圈數(shù)字編號說明內(nèi)容見正文功能介紹)

具體功能按數(shù)字編號依次介紹如下。

(1) 實(shí)現(xiàn)局部奇異性常規(guī)方法和迭代方法，并支持缺失數(shù)據(jù)計(jì)算。

(2) 支持各向同性及各向異性鄰域形狀并自動構(gòu)建多尺度窗口。

在柵格模型下直接以像元空間分辨率為單位來表達(dá)尺度大小，尺度大小為窗口邊長所占的像元數(shù)，正規(guī)化尺度參數(shù)L也采用此方案。鄰域形狀可選方形、圓形、矩形、橢圓，有寬度和高度2個(gè)維度。窗口數(shù)、最內(nèi)窗口、半徑增量、旋轉(zhuǎn)角度幾項(xiàng)關(guān)鍵參數(shù)來創(chuàng)建多尺度鄰域形狀，其中，最內(nèi)窗口設(shè)置值約定取奇數(shù)，默認(rèn)為1，也可以設(shè)置成3、5等其他奇數(shù)值，此時(shí)的最小尺度大于像元的空間分辨率。多尺度窗口序列有2種自動構(gòu)建方式：線性等間距(Linear space)，形成等差數(shù)列；對數(shù)等間距(Log space)，形成等比數(shù)列。在自動生成的基礎(chǔ)上，用戶可以自由修改尺度序列。

程序支持各向同性情形的奇異性分析，通過設(shè)置非等寬高的窗口參數(shù)(還可包括旋轉(zhuǎn)角度)來實(shí)現(xiàn)。開發(fā)了某個(gè)像元是否在1個(gè)橢圓內(nèi)(帶旋轉(zhuǎn)角度)的快速判斷算法。

程序支持方向性奇異性分析計(jì)算。程序由二維自動退化為一維情形，即E=1。選擇矩形窗口(可設(shè)置旋轉(zhuǎn)角度)，固定其中一個(gè)維度的窗口大小不變化。例如，在高度方向窗口數(shù)設(shè)置為1或者3等定值(長度單位為柵格空間分辨率)且尺度半徑增量固定為0，于是奇異性指數(shù)沿著東西方向進(jìn)行計(jì)算。若高度窗口數(shù)固定為1，相當(dāng)于原始柵格數(shù)據(jù)被分成多條東西向的一維數(shù)據(jù)序列；若高度窗口固定為3，則等效于原始柵格數(shù)據(jù)作南北向3個(gè)窗口的滑動平均后的多條一維數(shù)據(jù)序列。此時(shí)的奇異性指數(shù)圖可突出南北向的差異性。

(3) 支持自定義計(jì)算位置，包括2個(gè)方面。圖2中3a：自定義感興趣區(qū)域進(jìn)行計(jì)算，有3種選擇：數(shù)據(jù)格網(wǎng)所有位置、數(shù)據(jù)格網(wǎng)有效數(shù)據(jù)為準(zhǔn)、用戶文件自定義位置；圖2中3b：指定某一空間位置進(jìn)行奇異性探查與預(yù)測(圖3)。

圖3 局部奇異性分析特定計(jì)算位置奇異性探查與預(yù)測

(4) 支持批量計(jì)算，包括2個(gè)方面。圖2中4a：對多套不同的多尺度窗口進(jìn)行批量計(jì)算；圖2中4b：對多個(gè)輸入柵格數(shù)據(jù)進(jìn)行批量計(jì)算。

(5) 輸出結(jié)果多樣化，滿足進(jìn)一步分析的需要。結(jié)果文件保存不局限于奇異性指數(shù)，還包括系數(shù)c集，擬合優(yōu)度等多種指標(biāo)，可以追蹤歷次迭代的詳細(xì)結(jié)果，也可以只保留最終結(jié)果。

此外，程序還提供了友好的柵格數(shù)據(jù)管理功能(支持ArcGIS和Surfer柵格數(shù)據(jù)批量相互轉(zhuǎn)換)以及地圖代數(shù)常見分析功能。

為幫助用戶優(yōu)化參數(shù)設(shè)置，程序提供了單個(gè)位置的計(jì)算功能，用戶可以獲得該位置的感興趣信息：數(shù)據(jù)子集提取、鄰域窗口統(tǒng)計(jì)、雙對數(shù)回歸分析、奇異性模型計(jì)算結(jié)果、因變量預(yù)測。圖3顯示W(wǎng)X數(shù)據(jù)集在行列坐標(biāo)(53,75)位置處的奇異性計(jì)算信息，該位置計(jì)算所得Δα值為0.67(對應(yīng)α值為1.33)，表明具有較強(qiáng)的正奇異性，擬合優(yōu)度系數(shù)達(dá)到95.6%，可以通過顯著性檢驗(yàn)。按照所建立的奇異性模型，當(dāng)窗口大小為lcellsize/2和lcellsize/3時(shí)，預(yù)測結(jié)果分別為22.67，29.70，明顯高于原始柵格像元值20.41，與該位置的正奇異性特征一致，這提供了一種非線性插值方法。對應(yīng)窗口大小為3.5個(gè)像元單位時(shí)，也可以通過奇異性公式獲得預(yù)測值，為6.20，該預(yù)測值為分?jǐn)?shù)維窗口內(nèi)的滑動平均估值。

在前述參數(shù)設(shè)置條件下，增加迭代計(jì)算約束條件：標(biāo)準(zhǔn)誤差s(Δαk)≤0.01，最多迭代次為25次，對整個(gè)數(shù)據(jù)集計(jì)算全圖范圍的局部奇異性(圖4)。程序運(yùn)行結(jié)果表明，當(dāng)?shù)降?0次時(shí)，s(Δαk)=0.009 7，于是終止迭代。圖4展示了20次迭代過程中計(jì)算結(jié)果空間模式和常見統(tǒng)計(jì)量的變化趨勢。圖4曲線展示了s(Δαk)隨迭代次數(shù)k的變化，可見由s(Δα0)到s(Δα1)具有最大的變化性，迭代4～5次后，s(Δαk)的下降速度越來越緩慢，趨近于0。與之相對應(yīng)，局部系數(shù)c圖也變得越來越光滑。原始數(shù)據(jù)是一個(gè)典型的多重分形數(shù)據(jù)集，具有強(qiáng)烈的不均勻性，對原始數(shù)據(jù)制圖，僅能分辨有限的若干高值(顯性異常，形如“大型山峰”)，而更多局部鄰域上相對高值(弱緩異常，形如“小型山峰”)難以識別，將數(shù)據(jù)轉(zhuǎn)換成對數(shù)后在一定程度上達(dá)到了異常增強(qiáng)效果，可以識別較多的局部異常。通過局部奇異性常規(guī)分析所獲得的Δα0圖可見，原始數(shù)據(jù)還有更多空間位置具有局部弱緩異常存在。通過迭代方法，可以獲得任意迭代次數(shù)后的最終奇異性指數(shù)，由圖4下面所展示的3個(gè)結(jié)果來看，異常的總體面貌雷同，但清晰程度各有差別。實(shí)際上，欲獲得較好的奇異性指數(shù)，結(jié)合s(Δαk)曲線的變化性，對該數(shù)據(jù)集迭代4～5次就可得到較優(yōu)的結(jié)果。

此外，對于批量處理等功能，限于篇幅，在此不一一展示。綜上所述，所研發(fā)的基于柵格數(shù)據(jù)的局部奇異性分析程序計(jì)算效率高，功能全面，實(shí)用性強(qiáng)。核心算法代碼還編譯生成了動態(tài)鏈接庫，提供給第三方使用。軟件GeoDAS就采用了本程序的動態(tài)鏈接庫，將局部奇異性計(jì)算與可視化結(jié)合到1個(gè)環(huán)境中，促進(jìn)了方法的推廣并方便了用戶的使用。

圖4 局部奇異性迭代分析結(jié)果空間模式和常見統(tǒng)計(jì)量趨勢(枚舉了迭代次數(shù)為0, 1,10,20的c圖和Δα圖，圖件由Surfer軟件生成，圖件渲染顏色僅表示該圖自身取值的相對高低，不同圖件相同顏色值所對應(yīng)的數(shù)據(jù)值無對比性，常見統(tǒng)計(jì)量見左統(tǒng)計(jì)表)

5 結(jié) 論

在深入在剖析局部奇異性分析基本思想及計(jì)算方法的基礎(chǔ)上，引入正規(guī)化尺度參數(shù)L，提出了迭代方法來改進(jìn)局部奇異性指數(shù)的估值。迭代方法使得在計(jì)算中通常被忽視的局部系數(shù)c值的作用被挖掘并利用。迭代分析方法的成立是建立在正規(guī)化尺度參數(shù)L大于柵格空間分辨率的基礎(chǔ)上，通過迭代，可以獲得越來越光滑的c圖，與此相伴，在Δα中將積累更多的奇異性信息，因此迭代分析方法對Δα估值有其良好的特性。不同的L參數(shù)取值，可以獲得不同的c值，增加了奇異性模型的靈活性。基于建立的奇異性模型，可以開展非線性空間插值、滑動平均估值。所設(shè)計(jì)的軟件功能強(qiáng)勁，為用戶應(yīng)用迭代分析新方法提供了高效便捷的計(jì)算工具。不僅支持各向同性情形，還適用于各向異性尺度的奇異性計(jì)算和方向性奇異性計(jì)算。軟件的動態(tài)鏈接庫版本也在GeoDAS軟件中被推廣應(yīng)用，促進(jìn)了方法的推廣。

本研究提出的迭代方法的數(shù)學(xué)公式同樣適用于矢量數(shù)據(jù)模型。矢量數(shù)據(jù)作為最主要的空間數(shù)據(jù)類型之一，由于矢量數(shù)據(jù)空間位置的不規(guī)則性，相比柵格數(shù)據(jù)模型，需要優(yōu)化算法提高時(shí)間效率，這將是后續(xù)研究的方向之一。目前研究的迭代策略是全圖幅同步進(jìn)行的，迭代使用的鄰域形狀也是全區(qū)一樣，這些方面也可能進(jìn)行局部化而優(yōu)化迭代分析的結(jié)果。迭代方法中局部系數(shù)c集在迭代過程中間接地獲取更大范圍上的數(shù)據(jù)信息，具有加權(quán)滑動平均的效應(yīng)。隨著迭代次數(shù)的增加，各空間位置的ck(x)值被磨光的速率可能是不一致的，全圖幅同步迭代將導(dǎo)致一些位置被過度計(jì)算。迭代次數(shù)越多邊緣效應(yīng)越嚴(yán)重，這在解釋數(shù)據(jù)時(shí)需要引起注意。迭代過程中使用固定的鄰域形狀、固定的尺度變化范圍，簡化了計(jì)算，但忽視了局部系數(shù)在迭代過程中自身所發(fā)生的變化。本程序提供的批處理功能方便了用戶從多種計(jì)算方案中優(yōu)選結(jié)果，但這還不夠智能化，更好的方案是從空間數(shù)據(jù)自身特點(diǎn)著手來發(fā)現(xiàn)最合適的窗口形狀及尺度變化范圍，例如借助空間U統(tǒng)計(jì)量方法(Cheng et al, 1996, 2007)。這些問題、難題的解決將有助于基于多重分形的局部奇異性方法的進(jìn)一步發(fā)展與完善。

6 致 謝

在寫作過程中，筆者與Agterberg博士進(jìn)行了有益的探討，在此表示由衷的感謝！

成秋明.2000.多重分形理論和地球化學(xué)因素分布規(guī)律[J].地球科學(xué)：中國地質(zhì)大學(xué)學(xué)報(bào),25(3):311-318.

成秋明.2004.空間模式的廣義自相似性分析和礦產(chǎn)資源評價(jià)[J].地球科學(xué)：中國地質(zhì)大學(xué)學(xué)報(bào),29(6):733-744.

陳颙,陳凌.2005.分形幾何學(xué)[M].2版.北京:地震出版社.

成秋明.2006.非線性成礦預(yù)測理論:多重分形奇異性-廣義自相似性-分形譜系模型與方法[J].地球科學(xué)：中國地質(zhì)大學(xué)學(xué)報(bào),31(3):337-348.

成秋明.2007.成礦過程奇異性與礦產(chǎn)預(yù)測定量化的新理論與新方法[J].地學(xué)前緣,14(5):42-53.

陳志軍.2007.多重分形局部奇異性分析方法及其在礦產(chǎn)資源信息提取中的應(yīng)用[D].武漢:中國地質(zhì)大學(xué)(武漢).

陳志軍,成秋明,陳建國.2009.利用樣本排序方法比較化探異常識別模型的效果[J].地球科學(xué)：中國地質(zhì)大學(xué)學(xué)報(bào),34(2):353-364.

成秋明.2012.覆蓋區(qū)礦產(chǎn)綜合預(yù)測思路與方法[J].地球科學(xué)：中國地質(zhì)大學(xué)學(xué)報(bào)，37(6):1109-1125.

AGTERBERG F P. 2001. Multifractal simulation of geochemical map patterns[J]. Journal of China University of Geosciences, 12(1): 31-39.

AGTERBERG F P. 2003. Past and future of mathematical Geology[J]. Journal of China University of Geosciences, 14(3): 191-198.

AGTERBERG F P. 2012. Sampling and analysis of chemical element concentration distribution in rock units and orebodies[J]. Nonlin Processes Geophys, 19: 23-44.

CHENG QIUMING, AGGTERBERG F P, BONHAM-CARTER G F. 1996.Spatial analysis method for geochemical anomaly separation[J]. Journal of Geochemical Exploration, 56(3): 183-195.

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

CHEN ZHIJUN, CHENG QIUMING, CHEN JIANGUO. 2005. Significance of fractal measure in local singularity analysis of multifractal model[C]//CENG QIUMING, GRAEME BONHAM-CATER. Proceedings ofIAMG'05: GIS and Spatial Analysis. International Association for Mathematical Geology. Wuhan, China: China University of Geosciences Printing House, (1): 475-480.

CHENG QIUMING. 2006. GIS-based fractal/multifractal anomaly analysis for modeling and prediction of mineralization and mineral deposits[C]//HARRIS J R. GIS Application in the Earth Sciences. St John’s, NL, Canada:Geological Association of Canada, Special book: 289-300.CHEN ZHIJUN, CHENG QIUMING. 2007. Generalized local singularity analysis and its application[J]. Journal of China University of Geosciences, 18(Special issue): 348-350.

CHEN ZHIJUN, CHENG QIUMING, XIE SHUYUAN, et al. 2007. A novel iterative approach for mapping local singularities from geochemical data[J]. Nonlin Processes Geophys, 14: 317-324.

CHENG QIUMING, AGTERBERG F P. 2009. Singularity analysis of ore-mineral and toxic trace elements in stream sediments[J]. Computers & Geosciences, 35(2): 234-244.EVERTZ C J G, MANDEBROT B B. 1992. Multifractal measures: Appendix B[M]//PEITGEN H O, JURGENS H, SAUPE D. Chaos and Fractals. New York: Springer Verlag, 922-953.

Iterative approach to local singularity analysis and software implementation based on raster data

CHENZhi-jun1,2,CHENGQiu-ming1,2,3,CHENJian-guo1,2

(1. State Key Laboratory of Geological Processes and Mineral Resources, China University of Geosciences, Wuhan 430074, China; 2. Faculty of Earth Resources, China University of Geosciences (Wuhan), Wuhan 430074, China; 3. Department of Earth and Space Science, Department of Geography, York University, Toronto M3J1P3, Canada)

Weak anomaly identification was a challenging problem in the metallogeny prognosis and mineral resource assessment in recent years. Great attentions were drawn to the local singularity analysis based on multifractal theory for its new research ideas, simple method and beneficial effects to anomaly recognition. The authors introduced the singularity principles and its computational methods. The authors took into account the regularization scale parameter L and brought forward the iterative approach for the singularity analysis to improve the estimation of the singularity index. The iterative algorithm and software were realized and implemented using C++ language. This iterative approach could work not only for the isotropic case, but also could apply to the anisotropic scaling case and directional singularity index estimation. The dynamic-link library version of this software had been applied in GeoDAS which was professional GIS software for quantitative prediction of mineral resources.

Local singularity analysis; Raster data; Multifractals; Iterative approach; Regularization scale parameter

10.3969/j.issn.1674-3636.2014.04.613

2014-08-18；編輯：陸李萍

國家自然科學(xué)基金項(xiàng)目(41272361, 40802081),中國地質(zhì)調(diào)查局工作項(xiàng)目(12120113089300),中央高?；究蒲袠I(yè)務(wù)費(fèi)專項(xiàng)資金資助項(xiàng)目(CUG120102)

陳志軍(1978— )，男，副教授，博士，主要從事數(shù)學(xué)地質(zhì)的教學(xué)與科研工作, E-mail: chenzhijunCS@163.com

P628

：A

：1674-3636(2014)04-0613-10