礦產資源儲量估算邊界線垂直縱投影角的計算方法、Excel程序編制及應用研究

李光明，褚敬

(1.山東省地質科學研究院，山東 濟南 250013；2.山東魯能軟件技術有限公司，山東 濟南 250001)

地質技術人員在制作資源儲量估算垂直縱投影圖時，往往需要進行礦體邊界、斷層、邊坡、受護體移動面(線狀壓覆)、側伏線、礦界等投影。為此，通常要制作多個剖面圖，確定它們和礦體的交點位置，再將交點位置投影到縱投影圖上并進行連接，得到交線的投影線。這個過程費時費力，效率低下，增加了成本。該文通過計算的方法直接得到投影角，通過地面交點位置直接在垂直縱投影圖上畫出投影線，從而解決上述問題，大大簡化或省略了圖件制作過程。該文還提供了Excel計算程序，以提高計算效率。

1 方法原理

利用矢量和空間解析幾何方法，根據礦體和其他相交面(線)的產狀，建立平面法向量，推導和礦體交線的方向向量。再根據交線方向向量在3個坐標軸上的分量，求出交線的產狀(方位角和傾伏角)。最后根據交線產狀和投影面方位(礦體走向)，求出交線在垂直縱投影面上的投影角。

2 公式推導

2.1 求交線l的產狀

2.1.1 根據兩相交平面各自產狀，求交線l的產狀

已知兩平面π1,π2法向量分別為n1,n2;則π1,π2交線l的方向向量n應同時垂直于n1,n2，所以n=n1×n2。

設i,j,k分別為x,y,z軸上的單位矢量;Px1,Py1,Pz1為n1在3個坐標軸上的分量;Px2,Py2,Pz2為n2在3個坐標軸上的分量;Px,Py,Pz為n在3個坐標軸上的分量，則

(Py1×Pz2-Py2×Pz1)i+(Px2×Pz1-

Px1×Pz2)j+(Px1×Py2-Px2×Py1)k

即:

Px=Py1×Pz2-Py2×Pz1

(1)

Py=Px2×Pz1-Px1×Pz2

(2)

Pz=Px1×Py2-Px2×Py1

(3)

產狀為γi∠αi的平面，其法向量ni在3個坐標軸上的分量分別是：

Pxi=sinαi×cosγi

Pyi=sinαi×sinγi

Pzi=cosαi

則

Px1=sinα1×cosγ1

(4)

Py1=sinα1×sinγ1

(5)

Pz1=cosα1

(6)

Px2=sinα2×cosγ2

(7)

Py2=sinα2×sinγ2

(8)

Pz2=cosα2

(9)

交線l的產狀(圖1)：

傾伏向γ=tan-1(Py/Px)

(10)

(11)

圖1 兩平面交線l及垂直縱投影角φ

圖2 側伏角為δ的側伏線l

若求礦界π2和礦體交線的產狀，由于礦界為直立平面，傾角為90°，假設其走向為θ，則:

γ=θ

(12)

β=tan-1(tanα1×cos(γ1-θ))

(13)

2.1.2 根據礦體產狀和側伏角，求礦體邊界線(側伏線)l的產狀

設礦體產狀為γ1∠α1，側伏角為δ(圖2)，求礦體邊界線l的傾伏向γ(方位角)和傾伏角β。

l'=h/sinα1

l=l'/sinδ=h/(sinα1·sinδ)

β=sin-1(h/l)=sin-1(sinα1·sinδ)

ω=cos-1(tanα1/tanβ)

(14)

γ=γ1±ω=γ1±cos-1(tanα1/tanβ)

(15)

(15)式中，側伏方向向左取“+”號，向右取“-”號。(15)式結果若<0則加360°，>360°則減360°。

2.2 求交線的垂直縱投影角φ

在圖1中，礦體π1產狀為γ1∠α1，相交面π2產狀為γ2∠α2，礦體垂直縱投影面為π3，交線l的傾伏角為β，l與礦體傾向的夾角為ω，與π3的夾角為ω′，則:

tanβ=h/s

h=s·tanβ

cosω'=s'/s

s'=s·cosω'=s·cos(90-ω)=s·sinω

tanφ=h/s'=s·tanβ/(s·sinω)=tanβ/sinω

φ=tan-1(tanβ/sinω)

(16)

3 計算示例

表1中列出了礦體和斷層、邊坡、側伏線及礦界等情形的4個計算示例。表中8～13項計算公式見(4)～(9)，14～16項計算公式見(1)～(3)，17～21項計算公式見(14)～(15)，22項計算公式見(16)。投影角計算結果見第22項。

4 Excel程序及應用

4.1 Excel程序

為方便投影角計算，利用Excel編寫了計算程序,程序界面見圖3。

界面中將礦體垂直縱投影分為3種類型：一般類型、礦界類型和側伏線類型。礦界類型為計算礦界鉛垂面和礦體交線投影角的情形；側伏線類型為計算礦體側伏線投影角的情形；一般類型為其他情形，如計算礦體邊界、斷層、邊坡、受護體移動面(線狀壓覆)和礦體交線投影角的情形等。

軟件包括1個調用命令、1個主程序和3個函數：

調用命令CmdCount_Click()

主程序ProjectionAngle()

度轉換為弧度函數dtor()

弧度轉換為度函數rtod()

度轉換為度分秒函數dtodms()

各程序(函數)代碼如下：

表1 投影角計算示例

圖3 程序界面

4.1.1 調用命令

Private Sub CmdCount_Click()

Call ProjectionAngle

End Sub

4.1.2 主程序

Public nPI As Double '全局變量

Sub ProjectionAngle() '投影角計算

nPI = WorksheetFunction.Pi

'檢查數據的合規性

i=1 '序號

Do While Cells(i + 3, "C") <> Empty '礦體編號不能為空

nType = Cells(i + 3, "B") '類型：1-一般、2-礦界、3-側伏角

Select Case nType

Case 1 '一般情況

Cells(i + 3, "F") = Empty '無需輸入礦體側伏角

Cells(i + 3, "H") = Empty '無需輸入交面走向

'輸入0或未輸入值，系統都認為是Empty，結果是0。下面語句的作用：輸入的0為值，不判斷為“空”

'邏輯結果為1或0。4個參數單元格全部為真值時，結果為1+1+1+1=4。

lCond = (Cells(i + 3, "D") >= 0) * (Cells(i + 3, "D") <> "") + _

(Cells(i + 3, "E") >= 0) * (Cells(i + 3, "E") <> "") + _

(Cells(i + 3, "I") >= 0) * (Cells(i + 3, "I") <> "") + _

(Cells(i + 3, "J") >= 0) * (Cells(i + 3, "J") <> "")

If lCond< 4 Then

Range(Cells(i + 3, "A"), Cells(i + 3, "M")).Select

msg = MsgBox("第" &i + 3 & "行類型選擇有誤或參數不全！", vbOKOnly)

End

End If

If Cells(i + 3, "D") = Cells(i + 3, "I") And _

Cells(i + 3, "E") = Cells(i + 3, "J") Then

Range(Cells(i + 3, "A"), Cells(i + 3, "M")).Select

msg = MsgBox("第" &i + 3 & "行兩平面產狀相同！", vbOKOnly)

End

End If

Case 2 '礦界

Cells(i + 3, "F") = Empty '無需輸入礦體側伏角

Cells(i + 3, "I") = Empty '無需輸入交面傾向

Cells(i + 3, "J") = Empty '無需輸入交面傾角

lCond = (Cells(i + 3, "D") >= 0) * (Cells(i + 3, "D") <> "") + _

(Cells(i + 3, "E") >= 0) * (Cells(i + 3, "E") <> "") + _

(Cells(i + 3, "H") >= 0) * (Cells(i + 3, "H") <> "")

If lCond< 3 Then

Range(Cells(i + 3, "A"), Cells(i + 3, "M")).Select

msg = MsgBox("第" &i + 3 & "行類型選擇有誤或參數不全！", vbOKOnly)

End

End If

Case 3 '側伏角

Cells(i + 3, "H") = Empty '無需輸入交面走向

Cells(i + 3, "I") = Empty '無需輸入交面傾向

Cells(i + 3, "J") = Empty '無需輸入交面傾角

lCond = (Cells(i + 3, "D") >= 0) * (Cells(i + 3, "D") <> "") + _

(Cells(i + 3, "E") >= 0) * (Cells(i + 3, "E") <> "") + _

(Abs(Cells(i + 3, "F")) >= 0) * (Cells(i + 3, "F") <> "")

If lCond< 3 Then

Range(Cells(i + 3, "A"), Cells(i + 3, "M")).Select

msg = MsgBox("第" &i + 3 & "行類型選擇有誤或參數不全！", vbOKOnly)

End

End If

If Not (Abs(Cells(i + 3, "F")) > 0 And Abs(Cells(i + 3, "F")) < 90) Then '側伏角0～±90°

Range(Cells(i + 3, "A"), Cells(i + 3, "M")).Select

msg = MsgBox("第" &i + 3 & "行側伏角應介于0～±90°之間！", vbOKOnly)

End

End If

End Select

i = i + 1

Loop

'計算

i = 1

Do While Cells(i + 3, "C") <> Empty '礦體編號

nType = Cells(i + 3, "B") '類型：1-一般、2-礦界(鉛直面)、3-側伏角

nR1 = dtor(Cells(i + 3, "D")) '礦體傾向(弧度)

nA1 = dtor(Cells(i + 3, "E")) '礦體傾角(弧度)

nPitch = dtor(Cells(i + 3, "F")) '礦體側伏角(弧度)

nStrike = dtor(Cells(i + 3, "H")) '礦界方位角(弧度)

nR2 = dtor(Cells(i + 3, "I")) '交面傾向(弧度)

nA2 = dtor(Cells(i + 3, "J")) '交面傾角(弧度)

If nType = 1 Or nType = 2 Then '類型：1-一般、2-礦界

nPx1 = Round(Sin(nA1) * Cos(nR1), 4) 'Px1

nPy1 = Round(Sin(nA1) * Sin(nR1), 4) 'Py1

nPz1 = Round(Cos(nA1), 4) 'Pz1

If nType = 1 Then '一般情況

nPx2 = Round(Sin(nA2) * Cos(nR2), 4) 'Px2

nPy2 = Round(Sin(nA2) * Sin(nR2), 4) 'Py2

nPz2 = Round(Cos(nA2), 4) 'Pz2

nPx = Round(nPy1 * nPz2 - nPy2 * nPz1, 4) 'Px

nPy = Round(nPx2 * nPz1 - nPx1 * nPz2, 4) 'Py

nPz = Round(nPx1 * nPy2 - nPx2 * nPy1, 4) 'Pz

If (Cells(i + 3, "D") - Cells(i + 3, "I")) Mod 180 = 0 Then '傾向相互平行

If Cells(i + 3, "D") = Cells(i + 3, "I") Then

nTrend = Cells(i + 3, "D") + 90 '交線方位角(度)

Else

nTrend = (Cells(i + 3, "D") + Cells(i + 3, "I")) / 2 '交線方位角(度)

End If

If nTrend> 90 And nTrend<= 270 Then

nTrend = (nTrend + 180) Mod 360 '使方位角保持偏北或正東方向

End If

nObliquity = 0 '交線傾伏角

nProjectionAngle = 0 '投影角

Else

nTrend = Round(Atn(nPy / nPx), 4) '交線方位角(弧度)

nObliquity = Round(Atn(Abs(nPz) / Sqr(nPx * nPx + nPy * nPy)), 4) '交線傾伏角(弧度)

If (Cells(i + 3, "D") - Cells(i + 3, "I")) Mod 90 = 0 Then '傾向相互垂直

nProjectionAngle = Cells(i + 3, "J") '投影角(度)

Else

nProjectionAngle = rtod(Round(Atn(Tan(nObliquity) / Abs(Sin(Abs(nR1 - nTrend)))), 4)) '投影角(度)

End If

Select Case True

Case nPx< 0 And nPy> 0 '第Ⅱ象限(原nTrend<0)

nTrend = nPI + nTrend

Case nPx< 0 And nPy< 0 '第Ⅲ象限(原nTrend>0)

nTrend = nPI + nTrend

Case nPx> 0 And nPy< 0 '第Ⅳ象限(原nTrend<0)

nTrend = nPI * 2 + nTrend

End Select

nTrend = rtod(nTrend) '交線方位角(度)

nObliquity = rtod(nObliquity) '交線傾伏角(度)

'調整交線方位角，使其和礦體、交面傾向總體方位一致

nRmax = WorksheetFunction.Max(Cells(i + 3, "D"), Cells(i + 3, "I")) '較大傾向(度)

nRmin = WorksheetFunction.Min(Cells(i + 3, "D"), Cells(i + 3, "I")) '較小傾向(度)

If nRmax - nRmin< 180 And (nTrendnRmax) Then

nTrend = (nTrend + 180) Mod 360

Else

If nRmax - nRmin> 180 And (nTrend>nRmin Or nTrend

nTrend = 180 - nTrend

End If

End If

End If

Else '礦界

'如果礦界方位和礦體傾向夾角>90°，礦界方位應調轉180°

If (Abs(nStrike - nR1) >nPI / 2) And _

(Abs(nStrike - nR1)

If nStrike>= nPI Then

nStrike = nStrike - nPI

Else

nStrike = nStrike + nPI

End If

End If

nTrend = nStrike '交線方位角(弧度)

nObliquity = Round(Atn(Tan(nA1) * Cos(nR1 - nStrike)), 4) '交線傾伏角(弧度)

If nR1 = nTrend Then '礦界走向同礦體傾向

nProjectionAngle = 90 '投影角(度)

Else

nProjectionAngle = rtod(Round(Atn(Tan(nObliquity) / Abs(Sin(Abs(nR1 - nTrend)))), 4)) '投影角(度)

End If

nTrend = rtod(nTrend) '交線方位角(度)

nObliquity = rtod(nObliquity) '交線傾伏角(度)

End If

Else '3-側伏角

nObliquity = Round(WorksheetFunction.Asin(Sin(nA1) * Sin(Abs(nPitch))), 4)'側伏線傾伏角(弧度)

nSign = IIf(nPitch< 0, -1, 1)

nTrend = Round(nR1 + nSign * Atn(Tan(nA1) / (Sin(nA1) * Tan(Abs(nPitch)))), 4) '側伏線方位角(弧度)

nTrend = IIf(nTrend< 0, 2 * nPI + nTrend, _

IIf(nTrend> 2 * nPI, nTrend - 2 * nPI, nTrend)) '<0則加2π；>2π則減2π

nProjectionAngle = rtod(Round(Atn(Tan(nObliquity) / Abs(Sin(Abs(nR1 - nTrend)))), 4)) '投影角(度)

nTrend = rtod(nTrend) '交線方位角(度)

nObliquity = rtod(nObliquity) '交線傾伏角(度)

End If

Cells(i + 3, "K") = dtodms(nTrend) '交線方位角(度分秒)

Cells(i + 3, "L") = dtodms(nObliquity) '交線傾伏角(度分秒)

Cells(i + 3, "M") = dtodms(nProjectionAngle) '投影角(度分秒)

i = i + 1

Loop

End Sub

4.1.3 度轉換為弧度函數

Function dtor(degree) '度轉換為弧度

dtor = degree * nPI / 180

End Function

4.1.4 弧度轉換為度函數

Function rtod(radium) '弧度轉換為度

rtod = radium / nPI * 180

End Function

4.1.5 度轉換為度分秒函數

Function dtodms(degree) '度轉換為度分秒

Degrees = Int(degree) '度

minutes = Int((degree - Degrees) * 60) '分

sections = (degree - Degrees - minutes / 60) * 3600 '秒(包括小數)

Degrees = LTrim(Str(Degrees))

minutes = LTrim(Str(minutes))

section1 = Int(sections) '秒(整數)，數值

section2 = WorksheetFunction.Text(sections - section1, ".000") '秒的小數部分,如".012"

section1 = LTrim(Str(section1)) '秒(整數)，字符

If Len(minutes) = 1 Then minutes = "0" + minutes

If Len(section1) = 1 Then section1 = "0" + section1

dtodms = Degrees + "°" + minutes + "'" + section1 + section2 + Chr(34)

End Function

圖4 類型選擇組合框

4.2 程序應用

在圖3中輸入有關參數。

類型通過單擊組合框右側按鈕進行選擇(圖4)。

一般類型需輸入礦體編號、傾向、傾角和交面傾向、傾角；

礦界類型需輸入礦體編號、傾向、傾角和礦界方位(走向)；

側伏線類型需輸入礦體編號、傾向、傾角和側伏角。向左側伏的側伏角輸入正值，向右側伏的側伏角輸入負值。

名稱起標識作用，可以輸入，也可以不輸入。

圖5 計算結果

需要程序的讀者可聯系作者索取。