本發(fā)明屬于地理信息技術應用領域，具體涉及一種基于地形地質圖，自適應判定地表出露巖層產狀的方法。

背景技術：

巖層的產狀是指巖層在空間的產出狀態(tài)。除水平巖層成水平狀態(tài)產出外，一切傾斜巖層的產狀均以其走向、傾向和傾角表示，稱為巖層產狀三要素。產狀作為地質測量與成圖中常見的地質要素之一，在地質各領域有著廣泛應用，其獲取方法主要包括野外實測和間接計算兩種方法。

野外的巖層產狀測量方法，通常是使用地質羅盤、坡度儀等儀器，直接在出露的巖層層面上量測。但是，野外巖石的風化、層面和劈理、線理面區(qū)分困難等因素均會影響產狀的測量。有時受限于巖層出露情況，只能借鑒附近區(qū)域的數(shù)據(jù)，再通過室內計算方法間接獲得產狀數(shù)據(jù)。產狀間接計算中常用的方法，主要有“三點法”和“相鄰等高線法”，其原理為通過尋找共面且不共線的多點，確定平面的形態(tài)。在計算過程中，一方面，需要人工篩選合適的點和計算方法，并通過繪制輔助點和線，對產狀進行計算；另一方面，需要人為確定產狀符號位置。相關計算工作的人為化處理，導致計算效率低下，且產狀數(shù)據(jù)的質量也因人而異。

為此，本專利針對地表巖層，擬基于地形地質圖，實現(xiàn)一種地表巖層產狀的自適應計算方法。其重點在于自適應完成交點選擇、計算方法選擇和產狀點位選擇。

技術實現(xiàn)要素：

本發(fā)明所要解決的技術問題在于，克服現(xiàn)有技術存在的缺陷，提出了一種地表巖層產狀的自適應判定方法，基于地形地質圖，自適應完成交點選擇、計算方法選擇和產狀點位選擇。從而實現(xiàn)地表出露巖層產狀的自適應判定。

本專利的主要思路如下：1)提取和預處理巖層邊界線，計算巖層界線與等高線的交點；2)遍歷同一巖層界線上的鄰接交點，通過智能的交點選擇、計算方法選擇和產狀點位選擇等步驟，實現(xiàn)產狀的自適應判定；3)對產狀信息進行符號化和標注。

基于上述思路，本發(fā)明提出的一種地表出露巖層產狀的自適應判定方法，其步驟如下：

(一)交點集合獲取

步驟1.1：加載矢量格式的巖層面要素圖層和等高線線圖層，得到所有面要素集合Stra＝{s_i|i＝1,2,3,…,n}(其中：s_i為巖層面要素，n為巖層面要素的數(shù)量)和等高線的線要素集合CourLine＝{cl_i|i＝1,2,3,…,p}(其中：cl_i為等高線的線要素，p為等高線要素的數(shù)量)。

步驟1.2：讀取集合Stra中每個巖層面要素的邊界，將其預處理為單線要素，并剔除圖幅邊緣部分、第四紀巖層邊界部分、斷層界線部分以及不整合接觸界線部分，得到巖層層面邊界線集合SLine＝{sl_i|i＝1,2,3,…,c}，sl_i是巖層層面邊界線要素，c為所有巖層保留的層面出露邊界線的數(shù)量。

步驟1.3：針對任一邊界線要素sl_i遍歷等高線集合CourLine，求取線要素sl_i與所有等高線的交點集合InterP_i，對邊界線集合SLine內所有元素進行交點計算，得到集合InterP＝InterP₁∪InterP₂∪InterP₃∪…∪InterP_l，l為與等高線相交的邊界線數(shù)量。記錄交點的高程H，源巖層編號SID和源等高線編號CID，在同一集合InterP_i內，交點的SID屬性均為i。

(二)產狀自適應獲取與剔除

針對交點集合InterP_i＝{pi_j|j＝1,2,3,…,qi}，qi為交點集合InterP_i包含的交點數(shù)量，逐點進行條件判斷，獲取產狀a＝{ρ，θ，δ}，并確定產狀點位。其中：ρ為巖層走向，θ是巖層傾向，δ是巖層傾角。

由j＝1開始進行判斷，具體步驟為：

步驟2.1：讀取交點pi_j、pi_j+1和pi_j+2和pi_j+3(若下標超出qi，則使用pi_(j+1)％qi、pi_(j+2)％qi和pi_(j+3)％qi)，判別交點pi_j、pi_j+1、pi_j+2和pi_j+3是否滿足相鄰等高線法的下列選點規(guī)則；如果符合，執(zhí)行步驟2.4，如果不符合執(zhí)行步驟2.2：

a)交點pi_j和pi_j+3高程一致，交點pi_j+1和pi_j+2高程一致；

b)交點pi_j+1和pi_j+2的源等高線編號屬性CID一致；

c)向量和的最小夾角不大于閾值ε，即判定兩向量方位近似平行；閾值ε的設定一般不大于6°，兩線夾角在此范圍內可視為近似平行。閾值的設定越小，則平行要求越嚴格，可根據(jù)需求調整閾值以滿足不同要求。

步驟2.2：讀取交點pi_j、pi_j+1和pi_j+2(若下標超出qi，則使用pi_(j+1)％qi、pi_(j+2)％qi)，判別是否滿足三點法的下列選點規(guī)則；如果符合，執(zhí)行步驟2.6；如果不符合，執(zhí)行步驟2.3:

a)交點pi_j、pi_j+1和pi_j+2中至少兩點的高程不等；

b)交點pi_j、pi_j+1和pi_j+2不共線。

步驟2.3：重新選擇判斷起點。若原線要素sl_i為環(huán)線要素，則取d＝qi為界，若原線要素sl_i非環(huán)線要素，則取d＝qi-3為界。下標j增加1之后，若j不超過d，則返回步驟2.1；若j超過d，則交點集合InterP_i上的產狀信息獲取完成，執(zhí)行步驟2.8。

步驟2.4：相鄰等高線法計算產狀?；诮稽cpi_j(x_j,y_j,h₁)、pi_j+1(x_j+1,y_j+1,h₂)、pi_j+2(x_j+2,y_j+2,h₂)和pi_j+3(x_j+3,y_j+3,h₁)使用相鄰等高線法計算產狀：

a)對于交點pi_j、pi_j+1、pi_j+2和pi_j+3，計算平面向量和分

別取和方向上的單位向量和取和向量和的方位角度ρ為巖層走向；

b)根據(jù)公式(1)計算巖層傾向θ；

c)計算交點pi_j+1和交點pi_j+2的中點，記為輔助點ap₁(x₁,y₁,h₂)，過點ap₁，沿θ方向的

直線與向量v₂交于點ap₂(x₂,y₂,h₁)，根據(jù)公式(2)計算傾角δ；

步驟2.5：相鄰等高線法產狀點位自動確定。記直線l為過輔助點ap₁且方向為θ的直線，l與交點pi_j+1與pi_j+2之間的sl_i邊界線部分相交于點ap₃，則以ap₃坐標(x₃,y₃)作為產狀標注點位坐標，存儲產狀信息a。返回步驟2.2。

步驟2.6：三點法計算產狀?；诮稽cpi_j(x₁,y₁,h₁)，pi_j+1(x₂,y₂,h₂)，pi_j+2(x₃,y₃,h₃)使用三點法計算產狀：

a)將交點pi_j、pi_j+1和pi_j+2的坐標代入公式(3)，得到巖層層面方程(公式(4)所示，

其中D為任意常數(shù))的系數(shù)A、B和C；

Ax+By+Cz+D＝0 (4)

b)根據(jù)系數(shù)A、B，確定巖層傾向線方程(公式(5)所示)，其中M為任意常數(shù)；

Bx+Ay+M＝0 (5)

c)由公式(6)、(7)分別計算、修正傾向θ；

d)根據(jù)傾向θ，由公式(8)計算巖層走向；

ρ＝θ±90 (8)

e)根據(jù)系數(shù)A、B和C，由公式(9)計算巖層傾角δ；

步驟2.7：三點法產狀標注點位自動確定。以交點pi_j+1的坐標作為產狀點位坐標，存儲產狀信息a。返回步驟2.3。

步驟2.8：遍歷巖層邊界線sl_i對應的交點集合InterP_i上的交點，得到Att_i＝{ai_j|j＝1,2,3,…,ci}，其中ci為邊界線sl_i上的產狀點數(shù)量。對邊界線集合SLine內所有巖層邊界線sl_i進行產狀計算，得到產狀集合Att＝Att₁∪Att₂∪Att₃∪…Att_s，s為有效獲取了產狀的交點集合數(shù)量。

步驟2.9：對具明顯錯誤的產狀值進行剔除。以產狀集合Att_i為單位，基于構造特征結合地質知識對錯誤的產狀值進行剔除，得到篩選后的產狀集合Att_i＝{ai_j|j＝1,2,3,…,di}，其中di為邊界線sl_i上的糾錯后保留的產狀點數(shù)量。錯誤產狀值剔除步驟如下：

a)褶皺具有“兩翼產狀較穩(wěn)定，中間部分產狀過渡”的構造特征。結合“V字型法則”，整體判斷兩翼產狀傾向。具體表現(xiàn)為：同斜褶皺兩翼產狀同向，直立向斜兩翼產狀指向外部，直立背斜兩翼產狀指向內部。剔除褶皺上不滿足上述褶皺構造特征的產狀值；

b)若巖層厚度均勻，則兩側產狀基本為同向。對于巖層形狀簡單，兩側界線基本平行的巖層，首先可根據(jù)周圍構造或結合“V字型法則”大致確定一側界線上的產狀；繼而推測另一側界線上的產狀。剔除簡單巖層上不滿足其產狀特征的產狀值；

c)形態(tài)較平直的界線部分，若界線上產狀傾向不穩(wěn)定，統(tǒng)計界線上產狀傾向在不同方向區(qū)間上的頻率，保留高頻并連續(xù)的產狀區(qū)間數(shù)據(jù)。區(qū)間大小根據(jù)實際情況而定，在產狀個數(shù)較少的情況下，可取較大的區(qū)間。取一平直界線上的的16個產狀數(shù)據(jù)為例(如下表所示)，以10°為間隔對產狀傾向做雷達圖(如圖2所示)。其產狀數(shù)據(jù)傾向連續(xù)且較多分布于90°到130°之間，則保留90°到130°之間的產狀數(shù)據(jù)，剔除其他低頻、離散的數(shù)據(jù)值

步驟3、產狀數(shù)據(jù)符號化

步驟3.1：對于產狀點ai_j，過其坐標(axi_j,ayi_j)點，沿其傾向θ繪制單位長度的線，垂直于傾向θ繪制單位長度的線，完成產狀符號的繪制。

步驟3.2：將傾向θ由極坐標系轉為大地坐標系，對產狀要素傾角δ進行標注。

本發(fā)明方法針對地表出露巖層，基于地形地質圖，自適應完成交點選擇、計算方法選擇和產狀點位選擇的產狀判定，實現(xiàn)了地表出露巖層產狀的自適應判定，實用性強，適用性廣。

附圖說明

圖1本發(fā)明方法流程圖。

圖2示例產狀傾向雷達圖

圖3實施例實驗數(shù)據(jù)。

圖4巖層邊界線和等高線交點圖。

圖5示例交點圖。

圖6示例三點法自適應選點結果圖。

圖7示例相鄰等高線法自適應選點結果圖。

圖8實施例全圖三點法自適應選點結果圖。

圖9實施例全圖相鄰等高線法自適應選點結果圖。

圖10實施例全圖產狀符號化和標注結果圖。

圖11示例實測產狀數(shù)據(jù)與計算產狀數(shù)據(jù)對比圖。

具體實施方式

下面結合實施例和附圖，對本發(fā)明方法作進一步詳細說明。

實施例：選用某地1:50000地形地質圖(如圖2所示)數(shù)據(jù)，下面結合附圖，基于shp格式的地質體面圖層和等高線圖層，描述本方法方法實施過程并進一步說明本發(fā)明的效果。

(一)交點集合獲取

步驟1.1：加載巖層面要素圖層和等高線圖層，得到所有面要素集合Stra＝{s_i|i＝1,2,3,…,29}，等高線要素集合CourLine＝{cl_i|i＝1,2,3,…,21}。

步驟1.2：讀取集合Stra中巖層面要素的邊界線分別進行預處理，去除第四紀地層邊界線及地層與第四紀接觸的邊界線部分，并去除圖幅邊界的部分，得到邊界線集合SLine＝{sl_i|i＝1,2,3,…,47}。

步驟1.3：對邊界線集合SLine和等高線集合CourLine求取交點，邊界線要素sli上的交點存入一個集合InterP_i，得到交點集合InterP₁＝{(400875.27,3271721.62),(400957.22,3271801.49),(401056.83,3271892.77)}，InterP₂＝{(401235.55,3272047.50),(401334.05,3272128.06),(401434.51,3272198.33),(401551.99,3272129.90),(401497.99,3272039.65),(401398.55,3272039.69),(401376.85,3271795.00),(401606.74,3271790.62)}，…，InterP₄₄＝{(400144.21,3274334.45),(400102.31,3274263.53),(400051.63,3274186.33),(400050.79,3273975.12),(400131.35,3274013.82),(400196.16,3274044.96),(400278.58,3274088.92),(400354.36,3274135.09)}，交點記錄等高線高程H，源巖層編號SID和源等高線編號CID，在同一集合InterP_i內，交點的SID屬性一致為i。所有交點集合InterP＝InterP₁∪InterP₂∪InterP₃∪…∪InterP₄₄，包括交點398個(如圖3所示)。

(二)產狀自適應計算與剔除

針對交點集合自適應選點、計算及點位確定，在此設定向量平行判斷的角度閾值ε為6°。以巖層面要素s₇的邊界線為例，其邊界包括巖層邊界線sl₁₁、sl₁₂、sl₁₃和sl₁₄，與等高線分別存在16個、3個、0個和9個交點，則邊界線sl₁₁、sl₁₂和sl₁₄分別對應交點集合InterP₉＝{p9_j|j＝1,2,3,…,16}、InterP₁₀＝{p10_j|j＝1,2,3}和InterP₁₁＝{p11_j|j＝1,2,3,…,9}(如圖4所示)。以交點集合InterP₉為例，交點坐標值如下表所示。由交點p9₁開始進行判斷和計算；

步驟2.1：由交點p9₁開始，先對交點p9₁、p9₂、p9₃和p9₄進行相鄰等高線法的選點規(guī)則判斷，并不符合相鄰等高線法的選點規(guī)則；則對交點p9₁、p9₂和p9₃進行三點法的選點規(guī)則判斷，并不符合三點法的選點規(guī)則，則繼續(xù)判斷交點p9₂。

步驟2.2：由交點p9₂開始，對交點p9₂、p9₃、p9₄和p9₅進行相鄰等高線法的選點規(guī)則判斷，并不符合相鄰等高線法的選點規(guī)則；則對交點p9₂、p9₃和p9₄進行三點法的選點規(guī)則判斷，三點滿足如下條件：

1)三點的高程分別為250、250和300，高程不全相等；

2)三點不共線；

因此，三點滿足三點法選點規(guī)則，進行三點法計算。

步驟2.3：將交點p9₂、p9₆和p9₇代入公式計算，得到所在平面方程的系數(shù)：

計算得到產狀結果：

步驟2.4：產狀符號位置取交點p9₃的坐標，繼續(xù)判斷交點p9₃。

步驟2.5：由交點p9₃開始，對交點p9₃、p9₄、p9₅和p9₆進行相鄰等高線法的選點規(guī)則判斷，并不符合相鄰等高線法的選點規(guī)則；則對交點p9₃、p9₄和p9₅進行三點法的選點規(guī)則判斷，三點滿足三點法選點規(guī)則，進行三點法計算，得到產狀結果：

產狀符號位置取交點p9₄的坐標，繼續(xù)判斷交點p9₄。

由交點p9₄開始，對交點p9₄、p9₅、p9₆和p9₇進行相鄰等高線法的選點規(guī)則判斷，并不符合；則對交點p9₄、p9₅和p9₆進行三點法的選點規(guī)則判斷，三點滿足三點法選點規(guī)則，進行三點法計算，得到產狀結果：

產狀符號位置取交點p9₅的坐標，繼續(xù)判斷交點p9₅。

步驟2.6：由交點p9₅開始，對交點p9₅、p9₆、p9₇和p9₈進行相鄰等高線法的選點規(guī)則判斷，四點滿足：

1)交點p9₅和交點p9₈高程一致為350m；交點p9₆和交點p9₇高程一致為400m；

2)交點p9₆和交點p9₇的源等高線編號屬性CID一致；

3)交點p9₅指向交點p9₈的方位角為214.75°，交點p9₆指向交點p9₇的方位角為217.06°，兩方位角相差2.31°，小于閾值ε。

因此，四點滿足相鄰等高線法，利用相鄰等高線法進行產狀計算。計算得到產狀結果：

步驟2.7：計算得到產狀符號位置坐標(400450.92,3272412.92)，繼續(xù)判斷交點p9₅、p9₆和p9₇是否滿足三點法選點規(guī)則，直至遍歷完成。

步驟2.8：交點集InterP₉產狀計算完成，得到邊界線sl₁₁的產狀集合Att₆＝{a6_j|j＝1,2,3,…,15}。其中，使用三點法計算得到的有13個，其三點組合形式如圖5所示；使用相鄰等高線法計算得到的有2個，其四點組合形式如圖6所示。

步驟2.9：對所有InterP_i計算產狀，除去因不符合條件而沒有計算產狀值的交點集，最后得到Att＝Att₁∪Att₂∪…∪Att₂₆?？紤]位置重復且產狀值一致的產狀點為同一個，則共包含產狀點151個。其中，使用三點法計算得到的有132個，其三點組合形式如圖7所示；使用相鄰等高線法計算得到的有19個，其四點組合形式如圖8所示。通過對錯誤的產狀計算值進行剔除，得到114個產狀計算值(如圖9所示)。

(三)產狀數(shù)據(jù)符號化

加載實測產狀數(shù)據(jù)與計算產狀數(shù)據(jù)進行對比(如圖10所示)，實測數(shù)據(jù)包含32個實測產狀值。取實測值300m以內的計算值，或與實測值在同一巖層界線上的計算值與實測值作比較，對比結果如下表：

基于上表計算產狀值與實測產狀值的誤差值，將傾向、傾角劃分為“正確”、“尚可”和“錯誤”三個等級，統(tǒng)計其數(shù)量。

通過32個產狀實測值與其附近的計算值的對比，存在傾向或傾角計算結果誤差較大的產狀點有5個，傾向傾角完全計算錯誤的點位0個，其余數(shù)據(jù)誤差在可接受范圍內，總體正確率為84.4％。