面向構(gòu)造分析的赤平投影的計算機(jī)實現(xiàn)與應(yīng)用

毛先成

(1. 有色金屬成礦預(yù)測與地質(zhì)環(huán)境監(jiān)測教育部重點實驗室, 湖南 長沙 410006;2. 中南大學(xué)地球科學(xué)與信息物理學(xué)院, 湖南 長沙 410006)

0 引 言

傳統(tǒng)地質(zhì)構(gòu)造方面的研究有助于理解礦產(chǎn)資源的形成機(jī)制，某些礦體的產(chǎn)狀和形態(tài)與斷層之間存在直接的決定關(guān)系，斷層對產(chǎn)狀資源的分布有著重要的影響(李忠權(quán)等，2010)，礦液的匯集、滲透、分布和上涌需要節(jié)理提供空間條件(陳進(jìn)等，2018)。因此，研究礦床的地質(zhì)體數(shù)據(jù)能夠為礦產(chǎn)資源的勘探提供指導(dǎo)。

赤平投影是常用的地質(zhì)構(gòu)造分析方法，通過統(tǒng)計節(jié)理的傾向和傾角數(shù)據(jù)，可得到研究區(qū)節(jié)理的優(yōu)勢方位(陳麗芹等，2014)。目前，雖然赤平投影等地質(zhì)構(gòu)造方法的分析軟件種類繁多，但大多不是開源軟件，在處理數(shù)據(jù)時存在數(shù)據(jù)泄露的風(fēng)險，而且目前暫無系統(tǒng)介紹赤平投影在節(jié)理優(yōu)勢方位處理方面的資料。

因此，通過推導(dǎo)繪制赤平投影弧和節(jié)理優(yōu)勢方位計算模型算法，使用C++語言進(jìn)行編碼，實現(xiàn)赤平投影的可視化程序，并以膠西北多處金礦床為研究對象，驗證開發(fā)程序的正確性和實用性。

1 赤平投影原理與繪制

1.1 赤平投影

1.1.1 投影原理 赤平投影以投影球南北方向的頂點為極點，以極點為發(fā)射點，將投影球與空間中任意線狀或面狀實體相交的點或曲線投影到赤平面上。赤平投影具有簡單明了、直觀形象和綜合性強(qiáng)的特點，廣泛應(yīng)用于地球科學(xué)中(廖國華，1979)。

1.2 繪制赤平投影弧

1.2.1 計算圓交點坐標(biāo) 赤平投影弧是產(chǎn)狀數(shù)據(jù)所對應(yīng)的圓在基圓內(nèi)部的弧段，因此繪制赤平投影弧需要先計算產(chǎn)狀數(shù)據(jù)所對應(yīng)的圓與基圓的交點。計算機(jī)中的圓用半徑R和圓心坐標(biāo)(x,y)存儲，根據(jù)2個圓的半徑和圓心坐標(biāo)即可求出2個圓的交點。已知圓A和圓B相交于點C、D(圖3)，圓A的圓心坐標(biāo)為A(x1,y1)，半徑為r1；圓B的圓心坐標(biāo)為B(x2,y2)，半徑為r2。交點C和D的坐標(biāo)推導(dǎo)公式如下。

圖1 平面赤平投影原理Fig. 1 Stereographic projection principle

圖2 赤平投影圖構(gòu)成Fig. 2 Structure of stereographic projection

繪制輔助線(圖4)，其中直線AB=L，直線AB的斜率為K1，直線CD的斜率為K2，則有：

(1)

解得：

(2)

(3)

(4)

圖3 平面兩圓相交Fig. 3 Intersection of two circles on a plane

圖4 求兩圓交點示意圖Fig. 4 Calculating the points where two circles intersect

1.2.2 計算機(jī)繪制吳氏網(wǎng) 赤平投影常用吳氏網(wǎng)作為投影網(wǎng)。吳氏網(wǎng)由1個基圓和1組經(jīng)緯線構(gòu)成(圖5)，基圓圓周以北方向N為0°，順時針方向0°～360°，節(jié)理傾向方位與圓周方位對應(yīng)，圓周到圓心為0°～ 90°，節(jié)理的傾角方位與半徑方位對應(yīng)。

圖5 吳氏投影網(wǎng)Fig. 5 Wulff net

首先繪制基圓,其半徑為20 cm，在計算機(jī)制圖時，為方便尺度換算，半徑設(shè)為200 px。

繪制經(jīng)線，吳氏網(wǎng)中每條經(jīng)線均為圓心坐標(biāo)在X軸的半徑不同的圓在基圓內(nèi)部的弧段。以這些圓的圓心至C點(或D點)的距離為半徑，繪制一系列通過C、D點的圓弧，得到吳氏網(wǎng)一側(cè)的經(jīng)線，同理繪制另一側(cè)經(jīng)線(馮明權(quán)，2009)。設(shè)基圓的半徑為R，產(chǎn)狀數(shù)據(jù)投影到基圓的圓弧的圓心為F(圖6)，其坐標(biāo)為(x,y)，則:

圖6 投影網(wǎng)經(jīng)線繪制Fig. 6 Meridian drawing of the projection net

(5)

圓弧的半徑為r，長度為：

(6)

式(5)、(6)中，β為產(chǎn)狀平面的傾角，(°)。

將基圓與圓F的半徑和圓心坐標(biāo)代入公式(3)和(4)，求得交點C、D的坐標(biāo)，在基圓內(nèi)部的圓弧部分即為組成吳氏網(wǎng)的南北經(jīng)向部分。赤平投影弧的繪制僅用到吳氏網(wǎng)經(jīng)線的繪制方法，不再介紹計算機(jī)繪制吳氏網(wǎng)緯線的方法。

1.2.3 繪制任意產(chǎn)狀的赤平投影弧 在手工繪制某節(jié)理的赤平投影弧時，需要固定吳氏網(wǎng)，根據(jù)節(jié)理的走向旋轉(zhuǎn)相應(yīng)的角度，再根據(jù)傾角選取相應(yīng)南北走向的經(jīng)線作為該節(jié)理的赤平投影弧。

在計算機(jī)中，為避免手工繪制之類的操作，可采用以下方法，即在繪制吳氏網(wǎng)的南北向經(jīng)線時，每條經(jīng)線都代表傾向相同的節(jié)理，其圓心均在X軸上，只是每個節(jié)理的傾角不同，所以圓心的位置不同。

如圖7所示，以某節(jié)理赤平投影弧所在圓的圓心O2至基圓圓心O1的距離L12為半徑，以O(shè)1為圓心繪制圓C1，在圓C1上傾角處處相同，根據(jù)傾向可推導(dǎo)出該赤平投影弧所在圓的圓心O2的坐標(biāo)，推導(dǎo)公式如下。

圖7 任意產(chǎn)狀赤平投影弧Fig. 7 Stereographic projection arc of arbitrary occurrence

由公式(6)得到圓C1的半徑，即L12的長度，圓C1的標(biāo)準(zhǔn)方程為：

(7)

赤平投影弧所在圓的圓心在圓C1上，其圓心坐標(biāo)的推導(dǎo)公式為:

(8)

式(8)中,β為節(jié)理的傾向，(°)。

1.3 優(yōu)勢方位計算模型

1.3.1 繪制等密圖 計算節(jié)理組產(chǎn)狀數(shù)據(jù)的優(yōu)勢方位通常采用等密圖的方式。等密圖是屬性特征為密度的等值線圖，后者的繪制通常在規(guī)則格網(wǎng)或三角網(wǎng)中對離散點進(jìn)行插值處理，然后通過追尋等值線，最終繪制成圖(龔恒等，2015)。等值線圖的繪制方法已經(jīng)很成熟。

首先通過公式(9)將產(chǎn)狀數(shù)據(jù)投影到施密特網(wǎng)。施密特網(wǎng)與吳氏網(wǎng)的結(jié)構(gòu)相似(陳勝同等,1996)，不同之處在于施密特網(wǎng)從圓心到圓周對應(yīng)的是0°～90°，每條產(chǎn)狀數(shù)據(jù)在圖上均為一個離散點。為保證施密特網(wǎng)圓周附近的離散點能夠按照疊加對折原理處理，需將圓周邊緣的離散點對折到圓的另一端(余淳梅等，2008)。

(9)

將施密特網(wǎng)的外接矩形劃分為若干等分的規(guī)則格網(wǎng)，統(tǒng)計落入每個2×2單元大小格網(wǎng)內(nèi)切圓中的離散點個數(shù)，將離散點個數(shù)除以總散點數(shù)的結(jié)果作為該格網(wǎng)中心點位置的密度值；順序追蹤圖中密度值相同的格網(wǎng)中心點，連接所有等密度線，得到等密圖。

1.3.2 計算優(yōu)勢方位 追蹤得到的每條等密度線均包括一些格網(wǎng)中心點，通過對等密度線內(nèi)的格網(wǎng)中心點進(jìn)行處理，得到局部區(qū)域的優(yōu)勢中心。從若干條封閉型等值線中選取達(dá)到等值線最高值一半的等值線(李國慶等，2008)，統(tǒng)計落入這些等值線的格網(wǎng)點，計算坐標(biāo)的加權(quán)平均值，作為局部區(qū)域的優(yōu)勢中心。

設(shè)某條等值線所圍的格網(wǎng)點坐標(biāo)為(xi,yi),Hi為該格網(wǎng)點的密度值，則

(10)

求得的(x,y)即為施密特網(wǎng)下優(yōu)勢方位的直角坐標(biāo)，代入式(11)得到節(jié)理的優(yōu)勢方位(β,θ)。

(11)

2 赤平投影程序設(shè)計與實現(xiàn)

2.1 程序設(shè)計

2.1.1 交互設(shè)計 所設(shè)計程序在Microsoft Visual Studio 2012(vs2012)平臺進(jìn)行開發(fā)，在Windows桌面平臺運行。程序界面交互設(shè)計選擇微軟基礎(chǔ)類庫(MFC)框架，它不僅是一個以C++類的形式封裝的Windows API類庫，還是一個桌面程序開發(fā)框架。在vs2012開發(fā)平臺中創(chuàng)建一個MFC項目不需要考慮消息傳遞等底層邏輯，僅需考慮程序的實現(xiàn)邏輯,程序邏輯的實現(xiàn)使用C++語言，能夠直接繼承MFC函數(shù)庫中的類和方法。數(shù)據(jù)存儲使用CSV文本文件，它具有輕量級、易安裝的特點，非常適用于小型程序。程序的交互界面設(shè)計見圖8。

圖8 程序交互界面Fig. 8 Program interface

2.1.2 運行流程 (1) 讀取某研究區(qū)的節(jié)理數(shù)據(jù)，然后根據(jù)各節(jié)理的傾向和傾角數(shù)據(jù)繪制對應(yīng)的赤平投影弧。

(2) 對數(shù)據(jù)進(jìn)行處理，將產(chǎn)狀數(shù)據(jù)轉(zhuǎn)化為施式網(wǎng)直角坐標(biāo)系下的數(shù)據(jù)，并通過規(guī)則格網(wǎng)進(jìn)行統(tǒng)計。

(3) 根據(jù)設(shè)置的等值線等級，在規(guī)則格網(wǎng)中進(jìn)行等值線追蹤，選擇等值線中最高級別一半等級所包含的區(qū)域，計算優(yōu)勢方位。

(4) 繪制優(yōu)勢方位所對應(yīng)的赤平投影弧，并將可視化結(jié)果輸出。

程序運行流程見圖9。

圖9 程序運行流程圖Fig. 9 Program flow chart

2.2 關(guān)鍵過程實現(xiàn)

赤平投影的繪制過程主要由赤平投影弧繪制和優(yōu)勢方位計算組成。

2.2.1 赤平投影弧 僅存儲1條赤平投影弧，包含該弧段的起點和終點位置、寬度和顏色等屬性、計算交點及繪制方法等。具體結(jié)構(gòu)如下。

Class sgp_pjt_arc{

int m_r; //弧段所在圓的半徑

intm_width; //弧段的寬度

Color m_rgb; //弧段的顏色

double m_dip, m_dip_dir; //傾角和傾向

Point m_start_p, m_end_p; //起點終點位置

…

void draw_arc(); //繪制弧段

sgp_pjt_arc( double dip, double dip_dir );

//類的構(gòu)造函數(shù)

Pair cal_intersection_point(

double r1, double r2, Point p1, Point p2 ); //求兩個圓的交點

int cal_arc_radius(double R, double dip_dir ); //計算圓弧所在圓的半徑

…

}

2.2.2 優(yōu)勢方位計算模型 包括節(jié)理數(shù)據(jù)的轉(zhuǎn)換、極點統(tǒng)計與計算、等值線追蹤和優(yōu)勢方位計算等過程。具體結(jié)構(gòu)如下。

Classpri_orientiation{

Line m_max_countor; //級別最高等值線

Pair m_pri_ori; //優(yōu)勢方位

Array m_points; //存儲轉(zhuǎn)化后直角坐標(biāo)數(shù)據(jù)

Array m_dip_arr,m_dip_dir_arr; //存儲傾向傾角數(shù)據(jù)

Array < Pair m_grid; //規(guī)則格網(wǎng)，包括每個格網(wǎng)點的密度

Array > m_countors; //存儲追蹤到的等值線，包括等值//線的等級

…

Array > trace_countor(

Array < pair > grid ); //等值線追蹤

voidstatistics_density( Array points ); //極點密度統(tǒng)計

Array coor_transf(

Array dip, Array dip_dir ); //節(jié)理數(shù)據(jù)轉(zhuǎn)換

Pair cal_ pri_orientiation ( Line countor ); //計算優(yōu)勢方位

…

}

3 赤平投影應(yīng)用實例

3.1 研究區(qū)地質(zhì)背景

膠西北金礦礦集區(qū)礦產(chǎn)資源豐富，已探明金礦資源儲量占膠東地區(qū)總儲量的90%以上，中型以上金礦床達(dá)28處，膠東地區(qū)的特大型金礦均產(chǎn)于該地區(qū)(宋明春等，2011)。膠東金礦床具有區(qū)域集中、規(guī)模大、富集強(qiáng)度高和成礦期短的特點，區(qū)內(nèi)包括三山島斷裂帶、焦家斷裂帶、招平斷裂帶3條主要含礦斷裂帶(圖10)。

三山島斷裂帶位于膠西北最西部三山島—倉上一線，陸地延伸約10 km，斷裂的兩側(cè)延入渤海。該斷裂帶嚴(yán)格控制著一系列金礦床的產(chǎn)出，如三山島、倉上等金礦(肖風(fēng)利等，2018)。斷裂的總體走向為NE-NNE，局部為NE-NEE。

圖10 膠西北地質(zhì)簡圖(據(jù)Yang et al., 2016; Mao et al., 2019修改)1-第四系；2-上太古界膠東群；3-早白堊世郭家?guī)X花崗閃長巖；4-晚侏羅世欒家河二長花崗巖；5-晚侏羅世玲瓏花崗巖；6-區(qū)域斷裂Fig. 10 Simplified map of northwestern Jiaodong Peninsula(modified from Yang et al., 2016; Mao et al., 2019)

焦家斷裂帶位于平里店—黃山館一帶，總長度為27 km，寬80～500 m，區(qū)內(nèi)發(fā)育新城、紅布、寺莊等金礦床(蘇旭亮等，2016；王金利等，2020)。

招平斷裂帶北起龍口顏家莊，南至平度、洪山一帶，斷裂全長200 km，寬150～200 m，區(qū)內(nèi)發(fā)育大尹格莊、后倉、夏甸等金礦床，斷裂總體走向為NE-NNE，從北到南整體呈S形展布。

3.2 節(jié)理優(yōu)勢方位處理與分析

研究范圍包括三山島、寺莊和夏甸3個礦區(qū)，處理對象為3個礦區(qū)成礦期的節(jié)理數(shù)據(jù)。

(1) 繪制礦區(qū)節(jié)理數(shù)據(jù)的赤平投影弧(圖11)。圖中從左到右依次為三山島、寺莊和夏甸不含優(yōu)勢方位的赤平投影圖。

(2) 統(tǒng)計各礦區(qū)的節(jié)理數(shù)據(jù)，繪制各礦區(qū)的極點等密圖(圖12)。通過優(yōu)勢方位計算模型，得到各礦區(qū)優(yōu)勢方位的數(shù)據(jù)(表1)。

表1 各礦區(qū)節(jié)理優(yōu)勢方位

(3) 根據(jù)得到的優(yōu)勢方位，繪制各礦區(qū)節(jié)理數(shù)據(jù)的赤平投影圖(圖13)。

圖13顯示：三山島礦床位于三山島斷裂帶，走向為NE-NNE(圖13a)，處理得到的結(jié)果符合實際區(qū)域地質(zhì)構(gòu)造；寺莊礦床位于焦家主斷裂帶的南部，走向為NNE-NE(圖13b)，與焦家主斷裂的走向相同；夏甸礦床位于招平斷裂帶，結(jié)果顯示礦床斷裂走向為NE-NNE(圖13c)，與招平斷裂帶的走向一致。

Dips是一個輕量級、易操作的地質(zhì)繪圖軟件，常用于繪制玫瑰圖、極點圖和等密圖等。使用Dips軟件處理3個礦區(qū)的產(chǎn)狀數(shù)據(jù)，利用其等密圖功能計算3個礦區(qū)的優(yōu)勢方位，并與自編軟件處理得到的優(yōu)勢方位進(jìn)行比較(表2)。

由表2可知，自編軟件計算的優(yōu)勢方位與Dips軟件處理得到的結(jié)果誤差在±2°以內(nèi)，結(jié)果的精度符合實際需求，且各礦床的優(yōu)勢方位符合預(yù)期，與該礦床多數(shù)產(chǎn)狀數(shù)據(jù)的方位相似，所編寫程序繪制的赤平投影圖能夠滿足實際需求。

圖11 不含優(yōu)勢方位的赤平投影圖Fig. 11 Stereographic projection without dominant orientation (a) Sanshandao mining area; (b) Sizhuang mining area; (c) Xiadian mining area

圖12 各礦區(qū)節(jié)理等密圖Fig. 12 Isodense map of joints in each mining area(a) Sanshandao mining area; (b) Sizhuang mining area; (c) Xiadian mining area

圖13 各礦區(qū)節(jié)理數(shù)據(jù)赤平投影圖Fig. 13 Stereographic projection of joints in each mining area(a) Sanshandao mining area; (b) Sizhuang mining area; (c) Xiadian mining area

表2 各礦區(qū)節(jié)理優(yōu)勢方位Dips軟件處理結(jié)果

4 結(jié) 論

(1) 詳細(xì)推導(dǎo)了赤平投影弧的繪制和優(yōu)勢方位計算模型的算法，闡述了程序設(shè)計和實現(xiàn)的過程，提供了系統(tǒng)的赤平投影在節(jié)理優(yōu)勢方位處理方面的資料。

(2) 以膠西北3個金礦床為研究對象，使用自編程序繪制每個礦區(qū)的赤平投影圖，展示各礦區(qū)的優(yōu)勢方位，結(jié)果符合各礦區(qū)的區(qū)域地質(zhì)構(gòu)造特征。與Dips軟件處理結(jié)果相比，礦床優(yōu)勢方位誤差精度均在±2°以內(nèi)，驗證了所編程序的正確性和實用性。