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

，,3，

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

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

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

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

近些年，弱緩化探異常識別已成為成礦預(yù)測和勘查評價中十分關(guān)鍵的問題。多重分形理論的局部奇異性分析因其嶄新的思路、簡便的方法、優(yōu)良的效果而在弱緩異常識別中受到廣泛關(guān)注。在深入剖析局部奇異性分析基本思想及計算方法的基礎(chǔ)上，引入正規(guī)化尺度參數(shù)L，提出了迭代方法來改進局部奇異性指數(shù)的估值。給出了奇異性迭代算法并用C++編程實現(xiàn)，軟件功能強勁，操作靈活，不僅適用于各向同性情形，還適用于各向異性尺度的奇異性計算和方向性奇異性計算。軟件的動態(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ù)的非對稱奇異譜進行了成功建模，為現(xiàn)實世界中地球化學(xué)場多重分形性的存在性提供了理論依據(jù)。分形維數(shù)、多重分形對稱度、多重分形奇異譜等參數(shù)從非線性角度對描述成礦物質(zhì)富集與貧化的統(tǒng)計規(guī)律發(fā)揮了重要作用(成秋明, 2000, 2007)，然而，對于礦產(chǎn)資源定量預(yù)測中的“定位”(礦產(chǎn)資源空間分布位置)預(yù)測功能則有所不足，需要將全局性統(tǒng)計轉(zhuǎn)變?yōu)榫植炕Ｐ??；诙嘀胤中卫碚摻⒌木植科娈愋苑治瞿Ｐ?，試圖通過局部異常的有效識別來達到空間預(yù)測目的，該模型正逐漸在覆蓋區(qū)勘查地球化學(xué)弱緩異常識別中表現(xiàn)出巨大的應(yīng)用潛力，成為非線性礦產(chǎn)預(yù)測的有效方法之一(成秋明, 2012)。應(yīng)用的需求也促進了局部奇異性分析自身模型的發(fā)展，出現(xiàn)了多種形式的擴展模型(陳志軍, 2007)。主要從提高奇異性指數(shù)計算精度方面介紹局部奇異性分析迭代分析方法，設(shè)計算法并編程實現(xiàn)，同時指明準(zhǔn)確計算奇異性指數(shù)需要注意的一些關(guān)鍵參數(shù)設(shè)置。

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

奇異性理論指出：奇異性是指在很小的時間-空間范圍內(nèi)具有巨大能量釋放或巨量物質(zhì)形成的現(xiàn)象(成秋明, 2006; Cheng et al, 2009)。假定在E維(拓撲)空間中有一由質(zhì)點凝聚成分維數(shù)為D的分形(如礦物的富集)，L尺度上質(zhì)點集合的質(zhì)量M(L)與集合中質(zhì)點的數(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)計期望，式(2)于是寫成密度形式，有：

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

(4)

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

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

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

(5)

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

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

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

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

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

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

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

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

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

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

2.2.1ε=l具有長度量綱 此時，c值在理論上的量綱是[ρ][l]Δα，估計Δα值(或α值)和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)形式計算α(x)，再由式(3)轉(zhuǎn)換成Δα(x)。此時，回歸直線段斜率即代表α(x)，通常與柵格數(shù)據(jù)模型的維數(shù)2比較接近，然后由式(4)計算Δα(x)。當(dāng)li=1(帶長度單位)時，可以獲得擬合值即為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ù)模型中，可以直接用窗口邊長所占像元個數(shù)的比值：N(l)/N(L)等效表達物理長度比值l/L，這里，N(·)操作符表示取長度覆蓋的柵格像元數(shù)，允許取小數(shù)表示等效個數(shù)。對非方形窗口，借鑒面積與邊長的關(guān)系，N(l)的大小可由窗口覆蓋像元數(shù)的平方根來等效表達。此時，Δα(x)值與c值均無量綱，Δα(x)值大小與不消除量綱情形獲得的結(jié)果相同，而c值則不同。式(3)可轉(zhuǎn)變?yōu)椋?/p>

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

(8)

對x位置處奇異性指數(shù)及其系數(shù)的估計可由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ù)，也即ε=l/L=N(l)。對柵格數(shù)據(jù)計算窗口邊長通過數(shù)柵格像元個數(shù)在實際中操作簡單，僅需對計算窗口的行列號進行操作，不涉及具體空間距離計算，因此在奇異性指數(shù)計算中使用頻繁，此種情形c(x)值為原始數(shù)據(jù)的非線性擬合值。

(2)L=lmax(尺度窗口序列的最大值)，這種取值方式在許多分形計算場合常被采用。log(c(x))值在式(9)中為log-log圖中截距(此時li=L)，相當(dāng)于在lmax×lmax窗口滑動平均對數(shù)值的擬合值。在log-log擬合中，當(dāng)L取lmax時，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進行奇異性分析時，用c(x)值擬合相應(yīng)窗口大小的滑動平均值具有最佳的預(yù)測可信度。

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

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

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

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

3.1 迭代方法的基本思想

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

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

記α*(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*集可看作是ρ集的一個粗尺度逼近，其細節(jié)(粗糙度)信息保存在Δα*圖中，理論上的c*集應(yīng)不再具有多重分形特征，這樣ρ集的奇異性信息都應(yīng)集中反映在α*集中。若對c*集再做奇異性分析，則其計算所得的各個α值都趨近于同一個常數(shù)，即空間拓撲維數(shù)E。

基于同尺度的局部系數(shù)圖迭代方法可以不斷校正奇異性指數(shù)，將具有奇異性的c(x)通過迭代不斷演化成非奇異的(Chen et al, 2007)：首先將ρ(x)分解成c0(x)和α0(x)，被估計出來的c0(x)可以進一步估計其局部奇異性，即把它當(dāng)作新的ρ(x)，再次進行分解得到新的c1(x)，(下標(biāo)1表示進行第1次迭代數(shù)，下同)，可知c1(x)的離散程度比c0(x)小，比ρ(x)更小，不斷重復(fù)這樣的過程，可以預(yù)料，最終獲取的局部系數(shù)圖將趨近于一種平穩(wěn)狀態(tài)，最終將不再具有奇異性，此時它的局部奇異性指數(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ù)集，再次進行同樣的局部奇異性計算，如此循環(huán)，新產(chǎn)生的c集將變得越來越光滑，趨近于一種平穩(wěn)狀態(tài)……于是最終可獲得所期望的c*集。迭代分解的示意見圖1。

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

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

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

(13)

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

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

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

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

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

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

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

設(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

計算[i][j]柵格焦元位置m個不同

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

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

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

}

}

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

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

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

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

}

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

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

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

利用C++編程語言在Windows平臺下開發(fā)了RasterDataMining 1.0軟件，實現(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) 實現(xiàn)局部奇異性常規(guī)方法和迭代方法，并支持缺失數(shù)據(jù)計算。

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

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

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

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

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

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

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

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

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

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

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

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

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

5 結(jié) 論

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

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

6 致 謝

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

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

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

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

成秋明.2006.非線性成礦預(yù)測理論:多重分形奇異性-廣義自相似性-分形譜系模型與方法[J].地球科學(xué)：中國地質(zhì)大學(xué)學(xué)報,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é)報,34(2):353-364.

成秋明.2012.覆蓋區(qū)礦產(chǎn)綜合預(yù)測思路與方法[J].地球科學(xué)：中國地質(zhì)大學(xué)學(xué)報，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é)基金項目(41272361, 40802081),中國地質(zhì)調(diào)查局工作項目(12120113089300),中央高?；究蒲袠I(yè)務(wù)費專項資金資助項目(CUG120102)

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

P628

：A

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