礦產(chǎn)資源儲(chǔ)量估算邊界線垂直縱投影角的計(jì)算方法、Excel程序編制及應(yīng)用研究

李光明，褚敬

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

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

1 方法原理

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

2 公式推導(dǎo)

2.1 求交線l的產(chǎn)狀

2.1.1 根據(jù)兩相交平面各自產(chǎn)狀，求交線l的產(chǎn)狀

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

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

(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)

產(chǎn)狀為γi∠αi的平面，其法向量ni在3個(gè)坐標(biāo)軸上的分量分別是：

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的產(chǎn)狀(圖1)：

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

(10)

(11)

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

圖2 側(cè)伏角為δ的側(cè)伏線l

若求礦界π2和礦體交線的產(chǎn)狀，由于礦界為直立平面，傾角為90°，假設(shè)其走向?yàn)棣龋瑒t:

γ=θ

(12)

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

(13)

2.1.2 根據(jù)礦體產(chǎn)狀和側(cè)伏角，求礦體邊界線(側(cè)伏線)l的產(chǎn)狀

設(shè)礦體產(chǎn)狀為γ1∠α1，側(cè)伏角為δ(圖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)式中，側(cè)伏方向向左取“+”號(hào)，向右取“-”號(hào)。(15)式結(jié)果若<0則加360°，>360°則減360°。

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

在圖1中，礦體π1產(chǎn)狀為γ1∠α1，相交面π2產(chǎn)狀為γ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 計(jì)算示例

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

4 Excel程序及應(yīng)用

4.1 Excel程序

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

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

軟件包括1個(gè)調(diào)用命令、1個(gè)主程序和3個(gè)函數(shù)：

調(diào)用命令CmdCount_Click()

主程序ProjectionAngle()

度轉(zhuǎn)換為弧度函數(shù)dtor()

弧度轉(zhuǎn)換為度函數(shù)rtod()

度轉(zhuǎn)換為度分秒函數(shù)dtodms()

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

表1 投影角計(jì)算示例

圖3 程序界面

4.1.1 調(diào)用命令

Private Sub CmdCount_Click()

Call ProjectionAngle

End Sub

4.1.2 主程序

Public nPI As Double '全局變量

Sub ProjectionAngle() '投影角計(jì)算

nPI = WorksheetFunction.Pi

'檢查數(shù)據(jù)的合規(guī)性

i=1 '序號(hào)

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

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

Select Case nType

Case 1 '一般情況

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

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

'輸入0或未輸入值，系統(tǒng)都認(rèn)為是Empty，結(jié)果是0。下面語句的作用：輸入的0為值，不判斷為“空”

'邏輯結(jié)果為1或0。4個(gè)參數(shù)單元格全部為真值時(shí)，結(jié)果為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 & "行類型選擇有誤或參數(shù)不全！", 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 & "行兩平面產(chǎn)狀相同！", vbOKOnly)

End

End If

Case 2 '礦界

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

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 & "行類型選擇有誤或參數(shù)不全！", vbOKOnly)

End

End If

Case 3 '側(cè)伏角

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 & "行類型選擇有誤或參數(shù)不全！", vbOKOnly)

End

End If

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

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

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

End

End If

End Select

i = i + 1

Loop

'計(jì)算

i = 1

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

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

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

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

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

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) '交線傾伏角(度)

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

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°，礦界方位應(yīng)調(diào)轉(zhuǎn)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-側(cè)伏角

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

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

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

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 度轉(zhuǎn)換為弧度函數(shù)

Function dtor(degree) '度轉(zhuǎn)換為弧度

dtor = degree * nPI / 180

End Function

4.1.4 弧度轉(zhuǎn)換為度函數(shù)

Function rtod(radium) '弧度轉(zhuǎn)換為度

rtod = radium / nPI * 180

End Function

4.1.5 度轉(zhuǎn)換為度分秒函數(shù)

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

Degrees = Int(degree) '度

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

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

Degrees = LTrim(Str(Degrees))

minutes = LTrim(Str(minutes))

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

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

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

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 程序應(yīng)用

在圖3中輸入有關(guān)參數(shù)。

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

一般類型需輸入礦體編號(hào)、傾向、傾角和交面傾向、傾角；

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

側(cè)伏線類型需輸入礦體編號(hào)、傾向、傾角和側(cè)伏角。向左側(cè)伏的側(cè)伏角輸入正值，向右側(cè)伏的側(cè)伏角輸入負(fù)值。

名稱起標(biāo)識(shí)作用，可以輸入，也可以不輸入。

圖5 計(jì)算結(jié)果

需要程序的讀者可聯(lián)系作者索取。