基于ArcGIS的Fry法自動實現

倪永進,常杰*,田鵬飛,宋香鎖,梁吉坡,仵康林

(1.山東省地質科學研究院,山東 濟南 250013;2.山東省第一地質礦產勘查院,山東 濟南 250110)

0 引言

巖石有限應變測量系利用巖石內應變標志體的形態和分布來確定巖石的應變狀態,構成了現代構造地質學研究的重要手段。反映巖石變形前和變形后幾何形態差異的實體稱為應變標志物[1],常用的如化石、鮞粒等。近半個世紀以來,大量應變測量方法被提出[2-3],其中一類關注應變標志物單體或集合體變形前后的形態或方向變化,如反向輪法、Rf-φ法、莫爾圓圖解法[3-7],另一類關注應變標志物的相對位置關系,如心對心法[2-3]。Fry[8-9]將心對心法發展為一種簡便而優雅的測量方法,該方法通過平移統計應變標志物中心點,形成一個空心橢圓來代表應變橢圓。

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

相較于CorelDRAW、AI、Photoshop等作圖軟件,ArcGIS作為一種地理信息系統,具有良好的坐標和空間分析能力。雖然Fry法本身并不依賴地理坐標,但ArcGIS平臺便捷的坐標提取和空間處理功能為本方法的自動制圖和點密度運算提供了優良的工具選擇。選擇該平臺的另一個重要原因是,ArcGIS平臺提供了強大的python接口,其Arcpy功能包可以經由python語言調用而實現大量快速的重復計算。

本文利用python語言編制程序,調用ArcGIS平臺的空間處理功能,自動繪制Fry圖;并通過人機交互,選擇合理的應變橢圓參數范圍,采用點密度法求解最佳擬合的應變橢圓。

1 Fry法回顧

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

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

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

2 原理與方法

2.1 原理

基于ArcGIS的自動Fry法實現主體分為2個步驟:

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

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

最佳橢圓的自動檢驗標準采點計數密度法[10]。該方法首先構造一個內部橢圓和與之緊鄰的外部橢圓環,二者中心點、軸比、旋轉角和面積均相同(圖2)。當擬合良好時,內部橢圓對應空心橢圓,包含最低密度的點;外部橢圓環對應點密環帶,包含最高密度的點。因此可以通過在參數空間內的系統檢驗,使外-內密度比達到最大值。

圖2 點密度法原理圖(改自文獻[10])

2.2 具體實現

(1)將照片導入arcmap平臺,并創建點要素類,以點要素標記標志物的中心點。形成原始點集pori={pn(xn,yn)|n=1,2,3……}。

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

2個橢圓的特點為軸比相同、主軸方向相同,后者面積為前者的2倍。分別計算在橢圓1內的點數Nin和在橢圓2內但不在橢圓1內的點數Nout,計算后者的占比。

C=Nout/Nin

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

值得說明的是,在包含最佳取值且排除偽影取值的前提下,取值范圍的給定不影響最終結果,只會對計算速度產生較大影響。

3 應用實例與結果

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

表1 搜索空間

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

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

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

4 討論與結論

4.1 Fry法空心橢圓及其與應變的關系

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

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

(1)應變標志點應當足夠多。

(2)應變標志點的空間分布應當是隨機,或者說是反聚類的。

(3)巖石內的應變應當是均勻且可被觀察到的(反例如晶格位錯等)。

4.2 未來的發展方向

本方法通過人機交互,實現了基于ArcGIS的Fry法自動制圖和空心橢圓提取。但仍有部分值得改進之處,如應變標志點的自動提取和最佳橢圓的快速擬合。

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

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

致謝：本文成文過程中,山東省地質科學研究院地質研究所全體同事提供了支持。