毛先成

(1. 有色金屬成礦預(yù)測與地質(zhì)環(huán)境監(jiān)測教育部重點(diǎn)實(shí)驗(yà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ù)能夠?yàn)榈V產(chǎn)資源的勘探提供指導(dǎo)。

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

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

1 赤平投影原理與繪制

1.1 赤平投影

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

1.2 繪制赤平投影弧

1.2.1 計(jì)算圓交點(diǎn)坐標(biāo) 赤平投影弧是產(chǎn)狀數(shù)據(jù)所對應(yīng)的圓在基圓內(nèi)部的弧段，因此繪制赤平投影弧需要先計(jì)算產(chǎn)狀數(shù)據(jù)所對應(yīng)的圓與基圓的交點(diǎn)。計(jì)算機(jī)中的圓用半徑R和圓心坐標(biāo)(x,y)存儲，根據(jù)2個(gè)圓的半徑和圓心坐標(biāo)即可求出2個(gè)圓的交點(diǎn)。已知圓A和圓B相交于點(diǎn)C、D(圖3)，圓A的圓心坐標(biāo)為A(x1,y1)，半徑為r1；圓B的圓心坐標(biāo)為B(x2,y2)，半徑為r2。交點(diǎn)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 求兩圓交點(diǎn)示意圖Fig. 4 Calculating the points where two circles intersect

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

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

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

繪制經(jīng)線，吳氏網(wǎng)中每條經(jīng)線均為圓心坐標(biāo)在X軸的半徑不同的圓在基圓內(nèi)部的弧段。以這些圓的圓心至C點(diǎn)(或D點(diǎn))的距離為半徑，繪制一系列通過C、D點(diǎn)的圓弧，得到吳氏網(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)，求得交點(diǎn)C、D的坐標(biāo)，在基圓內(nèi)部的圓弧部分即為組成吳氏網(wǎng)的南北經(jīng)向部分。赤平投影弧的繪制僅用到吳氏網(wǎng)經(jīng)線的繪制方法，不再介紹計(jì)算機(jī)繪制吳氏網(wǎng)緯線的方法。

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

在計(jì)算機(jī)中，為避免手工繪制之類的操作，可采用以下方法，即在繪制吳氏網(wǎng)的南北向經(jīng)線時(shí)，每條經(jīng)線都代表傾向相同的節(jié)理，其圓心均在X軸上，只是每個(gè)節(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)勢方位計(jì)算模型

1.3.1 繪制等密圖 計(jì)算節(jié)理組產(chǎn)狀數(shù)據(jù)的優(yōu)勢方位通常采用等密圖的方式。等密圖是屬性特征為密度的等值線圖，后者的繪制通常在規(guī)則格網(wǎng)或三角網(wǎng)中對離散點(diǎn)進(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ù)在圖上均為一個(gè)離散點(diǎn)。為保證施密特網(wǎng)圓周附近的離散點(diǎn)能夠按照疊加對折原理處理，需將圓周邊緣的離散點(diǎn)對折到圓的另一端(余淳梅等，2008)。

(9)

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

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

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

(10)

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

(11)

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

2.1 程序設(shè)計(jì)

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

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

2.1.2 運(yùn)行流程 (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)計(jì)。

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

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

程序運(yùn)行流程見圖9。

圖9 程序運(yùn)行流程圖Fig. 9 Program flow chart

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

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

2.2.1 赤平投影弧 僅存儲1條赤平投影弧，包含該弧段的起點(diǎn)和終點(diǎn)位置、寬度和顏色等屬性、計(jì)算交點(diǎn)及繪制方法等。具體結(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; //起點(diǎn)終點(diǎn)位置

…

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 ); //求兩個(gè)圓的交點(diǎn)

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

…

}

2.2.2 優(yōu)勢方位計(jì)算模型 包括節(jié)理數(shù)據(jù)的轉(zhuǎn)換、極點(diǎn)統(tǒng)計(jì)與計(jì)算、等值線追蹤和優(yōu)勢方位計(jì)算等過程。具體結(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)，包括每個(gè)格網(wǎng)點(diǎn)的密度

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

…

Array > trace_countor(

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

voidstatistics_density( Array points ); //極點(diǎn)密度統(tǒng)計(jì)

Array coor_transf(

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

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

…

}

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

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

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

三山島斷裂帶位于膠西北最西部三山島—倉上一線，陸地延伸約10 km，斷裂的兩側(cè)延入渤海。該斷裂帶嚴(yán)格控制著一系列金礦床的產(chǎn)出，如三山島、倉上等金礦(肖風(fēng)利等，2018)。斷裂的總體走向?yàn)镹E-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ā)育大尹格莊、后倉、夏甸等金礦床，斷裂總體走向?yàn)镹E-NNE，從北到南整體呈S形展布。

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

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

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

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

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

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

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

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

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

圖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)勢方位計(jì)算模型的算法，闡述了程序設(shè)計(jì)和實(shí)現(xiàn)的過程，提供了系統(tǒng)的赤平投影在節(jié)理優(yōu)勢方位處理方面的資料。

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