矿产资源储量估算边界线垂直纵投影角的计算方法、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 计算结果

需要程序的读者可联系作者索取。