基于ArcGIS的Fry法自動(dòng)實(shí)現(xiàn)

倪永進(jìn),常杰*,田鵬飛,宋香鎖,梁吉坡,仵康林

(1.山東省地質(zhì)科學(xué)研究院,山東 濟(jì)南 250013;2.山東省第一地質(zhì)礦產(chǎn)勘查院,山東 濟(jì)南 250110)

0 引言

巖石有限應(yīng)變測(cè)量系利用巖石內(nèi)應(yīng)變標(biāo)志體的形態(tài)和分布來(lái)確定巖石的應(yīng)變狀態(tài),構(gòu)成了現(xiàn)代構(gòu)造地質(zhì)學(xué)研究的重要手段。反映巖石變形前和變形后幾何形態(tài)差異的實(shí)體稱為應(yīng)變標(biāo)志物[1],常用的如化石、鮞粒等。近半個(gè)世紀(jì)以來(lái),大量應(yīng)變測(cè)量方法被提出[2-3],其中一類關(guān)注應(yīng)變標(biāo)志物單體或集合體變形前后的形態(tài)或方向變化,如反向輪法、Rf-φ法、莫爾圓圖解法[3-7],另一類關(guān)注應(yīng)變標(biāo)志物的相對(duì)位置關(guān)系,如心對(duì)心法[2-3]。Fry[8-9]將心對(duì)心法發(fā)展為一種簡(jiǎn)便而優(yōu)雅的測(cè)量方法,該方法通過(guò)平移統(tǒng)計(jì)應(yīng)變標(biāo)志物中心點(diǎn),形成一個(gè)空心橢圓來(lái)代表應(yīng)變橢圓。

傳統(tǒng)的Fry法利用透明紙進(jìn)行手工制圖,操作過(guò)程簡(jiǎn)單但經(jīng)多次重復(fù),耗時(shí)長(zhǎng)、效率低,且在確定最終應(yīng)變橢圓時(shí)主觀性強(qiáng)、誤差大[10-12]。為解決這一問(wèn)題,過(guò)去30年中多種Fry法自動(dòng)實(shí)現(xiàn)方案被設(shè)計(jì)出來(lái),其中基于PASCAL[13-14]、Fortran[12,15]等語(yǔ)言由于設(shè)計(jì)年代較早,輸入、輸出以點(diǎn)的坐標(biāo)數(shù)據(jù)為主,對(duì)可視化支持較差;基于Objective-C[16]語(yǔ)言、CorelDRAW平臺(tái)[17]的僅解決了Fry圖的自動(dòng)生成,仍需人工確定最終的應(yīng)變橢圓。Waldron &Wallace提出了一種客觀擬合應(yīng)變橢圓的方法[10],即點(diǎn)密度法,該方法計(jì)算空心橢圓及與其具有相同面積的橢圓環(huán)中的點(diǎn)密度,并最小化這一比值來(lái)確定最佳擬合橢圓。

相較于CorelDRAW、AI、Photoshop等作圖軟件,ArcGIS作為一種地理信息系統(tǒng),具有良好的坐標(biāo)和空間分析能力。雖然Fry法本身并不依賴地理坐標(biāo),但ArcGIS平臺(tái)便捷的坐標(biāo)提取和空間處理功能為本方法的自動(dòng)制圖和點(diǎn)密度運(yùn)算提供了優(yōu)良的工具選擇。選擇該平臺(tái)的另一個(gè)重要原因是,ArcGIS平臺(tái)提供了強(qiáng)大的python接口,其Arcpy功能包可以經(jīng)由python語(yǔ)言調(diào)用而實(shí)現(xiàn)大量快速的重復(fù)計(jì)算。

本文利用python語(yǔ)言編制程序,調(diào)用ArcGIS平臺(tái)的空間處理功能,自動(dòng)繪制Fry圖;并通過(guò)人機(jī)交互,選擇合理的應(yīng)變橢圓參數(shù)范圍,采用點(diǎn)密度法求解最佳擬合的應(yīng)變橢圓。

1 Fry法回顧

在闡述自動(dòng)化實(shí)現(xiàn)方法前,有必要回顧Fry法的基本原理和手工作圖實(shí)現(xiàn)方法[8-9]。Fry法不考慮應(yīng)變標(biāo)志體本身的形態(tài),只衡量標(biāo)志體中心點(diǎn)相對(duì)位置[3],其基本思想是:理想狀態(tài)下,巖石截面可以簡(jiǎn)化為由半徑為r的圓緊密堆積構(gòu)成的集合體(圖1)。在該堆積狀態(tài)下,任意圓的圓心周圍都有6個(gè)距離為2r的點(diǎn)(B、C等)、6個(gè)點(diǎn)的距離為(D、E等)……構(gòu)成一系列同心圓,巖石均勻變形后,同心圓變?yōu)橥臋E圓。天然狀態(tài)下,雖然應(yīng)變標(biāo)志體半徑差異較大且堆積方式多變,但由于其在空間上的分布是隨機(jī)的,經(jīng)點(diǎn)的疊加仍可構(gòu)成近于完整的橢圓。

a—變形前緊密堆積的顆粒圓及同心圓;b—變形后緊密堆積的顆粒橢圓及同心橢圓

該方法操作過(guò)程中可能出現(xiàn)誤差的主要有2處,即標(biāo)志點(diǎn)的選擇和空心橢圓的確定。前者主要由于樣品表面不平整或亞顆粒的存在所造成,但基于大量的統(tǒng)計(jì),這種誤差的影響是微乎其微的。而空心橢圓作為Fry法最后確定應(yīng)變的直接依據(jù),通常由人為劃定,主觀性較強(qiáng)。大量工作討論過(guò)客觀自動(dòng)識(shí)別空心橢圓的方法[10-12],但也指出了這種自動(dòng)識(shí)別的關(guān)鍵挑戰(zhàn)--偽影[18],即一個(gè)與真實(shí)空心橢圓存在交角的、范圍較小的應(yīng)變標(biāo)志點(diǎn)空缺區(qū)。

2 原理與方法

2.1 原理

基于ArcGIS的自動(dòng)Fry法實(shí)現(xiàn)主體分為2個(gè)步驟:

(1)自動(dòng)生成Fry圖。本文中,將照片導(dǎo)入ArcGIS平臺(tái)后,創(chuàng)建原始點(diǎn)要素類,采用目視解譯,人工繪制標(biāo)志體中心點(diǎn)位置。之后以arcpy工具包讀取原始點(diǎn)要素類中各點(diǎn)的坐標(biāo),并按照Fry法,自動(dòng)繪制Fry圖。

(2)人機(jī)交互擬合最佳橢圓。由于其尺寸的局限性和顯著的角度偏差,偽影的存在通常不會(huì)給人工確定空心橢圓造成困擾,只會(huì)影響缺乏直觀感覺(jué)的自動(dòng)算法;與之相反,不同觀察者對(duì)空白橢圓軸比和旋轉(zhuǎn)角的判斷卻可能存在較大差異,這是自動(dòng)算法可以避免的。因此本文采用交互方法,先由人工依據(jù)第(1)步自動(dòng)生成的Fry圖解給出應(yīng)變橢圓可能的參數(shù)取值范圍,再由程序在該范圍內(nèi)進(jìn)行檢驗(yàn),確定最佳參數(shù)并生成橢圓。

最佳橢圓的自動(dòng)檢驗(yàn)標(biāo)準(zhǔn)采點(diǎn)計(jì)數(shù)密度法[10]。該方法首先構(gòu)造一個(gè)內(nèi)部橢圓和與之緊鄰的外部橢圓環(huán),二者中心點(diǎn)、軸比、旋轉(zhuǎn)角和面積均相同(圖2)。當(dāng)擬合良好時(shí),內(nèi)部橢圓對(duì)應(yīng)空心橢圓,包含最低密度的點(diǎn);外部橢圓環(huán)對(duì)應(yīng)點(diǎn)密環(huán)帶,包含最高密度的點(diǎn)。因此可以通過(guò)在參數(shù)空間內(nèi)的系統(tǒng)檢驗(yàn),使外-內(nèi)密度比達(dá)到最大值。

圖2 點(diǎn)密度法原理圖(改自文獻(xiàn)[10])

2.2 具體實(shí)現(xiàn)

(1)將照片導(dǎo)入arcmap平臺(tái),并創(chuàng)建點(diǎn)要素類,以點(diǎn)要素標(biāo)記標(biāo)志物的中心點(diǎn)。形成原始點(diǎn)集pori={pn(xn,yn)|n=1,2,3……}。

(3)經(jīng)觀察步驟(2)生成的圖像,輸入空白橢圓半長(zhǎng)軸a、半短軸b、旋轉(zhuǎn)角k可能的取值范圍rangea、rangeb、rangek。以半長(zhǎng)軸(a)、半短軸(b),生成2個(gè)橢圓,并旋轉(zhuǎn)角度(k):

2個(gè)橢圓的特點(diǎn)為軸比相同、主軸方向相同,后者面積為前者的2倍。分別計(jì)算在橢圓1內(nèi)的點(diǎn)數(shù)Nin和在橢圓2內(nèi)但不在橢圓1內(nèi)的點(diǎn)數(shù)Nout,計(jì)算后者的占比。

C=Nout/Nin

在rangea、rangeb、rangek內(nèi)組合取值,直至遍歷整個(gè)空間。由于內(nèi)部橢圓1和外部圓環(huán)的面積相等,因此面積項(xiàng)被約掉,二者范圍內(nèi)點(diǎn)的個(gè)數(shù)比即為密度比。選取c取最大值時(shí)內(nèi)部橢圓(橢圓1)的參數(shù),即為空白橢圓參數(shù)。

值得說(shuō)明的是,在包含最佳取值且排除偽影取值的前提下,取值范圍的給定不影響最終結(jié)果,只會(huì)對(duì)計(jì)算速度產(chǎn)生較大影響。

3 應(yīng)用實(shí)例與結(jié)果

選擇2個(gè)實(shí)例作為應(yīng)用實(shí)例進(jìn)行實(shí)驗(yàn)。實(shí)例1照片(圖3a)來(lái)自文獻(xiàn)[10],為加拿大英屬哥倫比亞地區(qū)上Kaza組砂巖,碳酸鹽膠結(jié),分選良好。在照片中采集到了142個(gè)點(diǎn),并生成了20022個(gè)點(diǎn)構(gòu)成的Fry圖。為尋找最合適的橢圓,搜索空間如表1,在a-b-k空間中共搜索了5408個(gè)點(diǎn),最終得到的橢圓a=1.7,b=0.9,c=1.15(圖3c),橢圓內(nèi)包含1個(gè)點(diǎn),橢圓外環(huán)包含點(diǎn)79個(gè),c值(橢圓外環(huán)點(diǎn)數(shù)占比)為98.73%。

表1 搜索空間

a—砂巖鏡下照片及中心點(diǎn);b—礫巖露頭照片及中心點(diǎn);c—砂巖的Fry圖解和最佳擬合的空心橢圓;d—礫巖的Fry圖和最佳擬合的空心橢圓

實(shí)例2照片(圖3b)來(lái)自文獻(xiàn)[12],為瑞典Gaddede地區(qū)礫巖,顆粒支撐,礫石長(zhǎng)軸方向近水平。該照片中采集到了68個(gè)點(diǎn),并生成了由4556個(gè)點(diǎn)構(gòu)成的Fry圖。注意目視該圖中空心橢圓近水平(0°),但在左上—右下方向(130°附近)存在一個(gè)偽影區(qū)。為了尋找最合適的橢圓,在a-b-k空間中搜索了13104個(gè)點(diǎn)(表1),最終得到的橢圓a=10.2,b=1.8,c=0(圖3d),橢圓內(nèi)包含2個(gè)點(diǎn),橢圓外環(huán)包含點(diǎn)61個(gè),c值為96.72%。

表1中a、b、k分別為擬合橢圓的半長(zhǎng)軸、半短軸和旋轉(zhuǎn)角度。由于Fry法對(duì)絕對(duì)長(zhǎng)度不敏感,僅關(guān)注a與b的比值,因此長(zhǎng)度單位可以是任意的。k采弧度值,以指向右為0,順時(shí)針旋轉(zhuǎn)為負(fù),逆時(shí)針旋轉(zhuǎn)為正。

4 討論與結(jié)論

4.1 Fry法空心橢圓及其與應(yīng)變的關(guān)系

通常認(rèn)為Fry圖空心橢圓以外應(yīng)當(dāng)有一個(gè)點(diǎn)密環(huán)帶,表征了距離應(yīng)變標(biāo)志點(diǎn)最近的質(zhì)點(diǎn)在變形后的分布。但從本文2個(gè)實(shí)例對(duì)比來(lái)看,實(shí)例1中的點(diǎn)密環(huán)帶相對(duì)清晰,c值可達(dá)到98.73%;實(shí)例2中點(diǎn)密環(huán)帶表現(xiàn)并不明顯,c值也低于實(shí)例1。該環(huán)帶的顯著性可能受2種因素制約,即應(yīng)變均勻程度和顆粒接觸關(guān)系。早期觀點(diǎn)認(rèn)為空心橢圓代表應(yīng)變標(biāo)志體的平均應(yīng)變,點(diǎn)密環(huán)帶代表全巖總應(yīng)變,只有空心橢圓而無(wú)點(diǎn)密環(huán)帶的Fry圖表示標(biāo)志體初始分布不均勻或變形不均勻[15,19]。另一種觀點(diǎn)則認(rèn)為該環(huán)帶的顯著性可能與顆粒的接觸關(guān)系有關(guān),即其在顆粒支撐環(huán)境下的顯著性要高于基質(zhì)支撐環(huán)境[17-18,20]。

這一現(xiàn)象也引發(fā)了Fry法局限性的討論。盡管該方法被用在多種巖性的應(yīng)變測(cè)量中,甚至是用以評(píng)估礦化空間的分布[21],但應(yīng)當(dāng)認(rèn)識(shí)到,只有當(dāng)符合如下假設(shè)時(shí)空心橢圓才可被視為應(yīng)變橢圓:

(1)應(yīng)變標(biāo)志點(diǎn)應(yīng)當(dāng)足夠多。

(2)應(yīng)變標(biāo)志點(diǎn)的空間分布應(yīng)當(dāng)是隨機(jī),或者說(shuō)是反聚類的。

(3)巖石內(nèi)的應(yīng)變應(yīng)當(dāng)是均勻且可被觀察到的(反例如晶格位錯(cuò)等)。

4.2 未來(lái)的發(fā)展方向

本方法通過(guò)人機(jī)交互,實(shí)現(xiàn)了基于ArcGIS的Fry法自動(dòng)制圖和空心橢圓提取。但仍有部分值得改進(jìn)之處,如應(yīng)變標(biāo)志點(diǎn)的自動(dòng)提取和最佳橢圓的快速擬合。

ArcGIS平臺(tái)強(qiáng)大的空間處理能力可以快速的識(shí)別照片(柵格數(shù)據(jù))中的顆粒(面),并提取應(yīng)變標(biāo)志體(面)的中心點(diǎn)(點(diǎn))。但在實(shí)際實(shí)驗(yàn)中,照片的識(shí)別普遍被干擾,且干擾源不可控,即不同照片的干擾因素差別巨大。典型的,對(duì)于露頭或手標(biāo)本尺度的照片,干擾因素可以包括露頭或標(biāo)本表面不平整造成的照片折線、應(yīng)變標(biāo)志物與基質(zhì)的顏色對(duì)比過(guò)小或邊界不清晰;對(duì)于顯微照片,干擾因素可以報(bào)考亞顆粒的存在、雙晶或多晶的邊界以及不均勻消光等。近幾年在計(jì)算機(jī)視覺(jué)和大數(shù)據(jù)地質(zhì)學(xué)領(lǐng)域,已有大量工作嘗試解決圖像識(shí)別問(wèn)題[22-25],但距離本工作的實(shí)際應(yīng)用尚有效率和準(zhǔn)確率問(wèn)題需要進(jìn)一步解決:在幾十至百余個(gè)點(diǎn)的數(shù)量級(jí)上,人工修改點(diǎn)與新建點(diǎn)的速度相差無(wú)幾。

本方法計(jì)算過(guò)程主要耗時(shí)之處在于遍歷可能的長(zhǎng)軸—短軸—旋轉(zhuǎn)角組合,由于反復(fù)調(diào)用arcpy的點(diǎn)—面生成方法、空間連接函數(shù)和讀取要素字段值,造成該步驟耗時(shí)過(guò)長(zhǎng)。在本文提供的方法中,該步驟相當(dāng)于生成一個(gè)由a-b-k和擬合度構(gòu)成的四維空間,在該空間內(nèi),前三者形成一個(gè)四維曲面,通過(guò)遍歷該曲面上的節(jié)點(diǎn),尋找到擬合度軸上的最大值點(diǎn)。而在當(dāng)前的機(jī)器學(xué)習(xí)領(lǐng)域,該步驟似乎可以通過(guò)梯度下降解決,即在該曲面上隨機(jī)選擇一點(diǎn),計(jì)算該點(diǎn)的梯度,并向梯度方向移動(dòng)一定距離,該距離與梯度的大小呈正比,通過(guò)該過(guò)程可避免計(jì)算全部曲面節(jié)點(diǎn),進(jìn)而大幅度提高計(jì)算效率。

致謝：本文成文過(guò)程中,山東省地質(zhì)科學(xué)研究院地質(zhì)研究所全體同事提供了支持。