一维电测深反演实验

煤田电法勘探数字煤田电法勘探数字 解释基础解释基础 实验报告实验报告 姓姓 名名： 龙龙 斌斌 学学 号：号： 05102344 班班 级：级： 地物地物 10-2 班班 指导老师：指导老师：姜志海、刘志新姜志海、刘志新 2013 年 11 月 中国矿业大学 实验一实验一 一维电测深反演实验一维电测深反演实验 一、一、实验目的实验目的 通过对课堂学习的一维电测深反演的基础理论知识，对一维电测深反演的方 法（拟合核函数法、直接拟合视电阻率法、DZ 法）进一步的加强学习并练习， 通过对程序的编制更深层次的理解一维电测深反演的基本思想与基本步骤。 二、二、实验内容实验内容 利用软件编制一维电测深反演程序，并对最终的成果图进行理解与解释，同 时对于数据采集、处理及解释的程序进一步的理解与体会。 三、三、实验原理实验原理 1. 拟合核函数法拟合核函数法 基本思想： 确定初始模型：（层参数初值） ； 0 P 由实测曲线视电阻率转换函数曲线； ( ) s r 数字滤法 T 由初值计算理论曲线； 0 P T 按照一定法则进行曲线对比，判断与拟合程度是否达到要求。根据拟合差 TT 确定是否进行层参数调查； 调整，以改善与的拟合程度； 0 P TT 重复步骤，直到与曲线达到拟合要求为止。 TT 2. 视电阻率转化函数视电阻率转化函数 基本思想与拟合核函数法相似。 3. DZDZ 法自动反演法自动反演 基本思想： 根据实测 曲线给出初值s 0 p 由 计算理论 曲线，计算 DZ 曲线（只计算 DZ 点坐标） 0 p s 比较 与 ，修正 DZ 曲线s s 由 DZ 曲线 计算新的模型参数 1 p 把 当作 ，重复 步骤，直到 与 之间拟合差 1 p 0 ps s 满足精度要求。 三、实验成果图三、实验成果图 1、模型参数设置： 2、运行曲线 3、反演结果 五、实验体会五、实验体会 通过对本次实验的学习，对一维电测深反演有了更深刻的，对 VB 的编写也有了一定 的提升，但是由于各方面原因，在实验编写的过程当中也遇到了很多不可避免的问题，比 如说，对反演运算过程的一些认识和程序的编写，最后通过请教其他同学才得以完成，总 的来说，这次实验的收获还是不小。 附加程序：附加程序： 一维电测深反演 VB 程序 Public zbx As Integer Public n As Integer Public xps() As Double Public nj1 As Integer Public nj2 As Integer Public ppsz() As Double Public ppsf()As Double Form1 部分 Dim x11As Double Dim x22 As Double Dim y11As Double Dim y22 As Double Private Sub Command1_Click() Dim nn As Integer Dim hh()As Double Dim pp()As Double Dim zpp() As Double Dim zhh() As Double nn = Val(Text1.Text) ReDim hh(nn - 1)As Double ReDim pp(nn)As Double ReDim zpp(nn)As Double ReDim zhh(nn - 1)As Double ss1 = Text2.Text ss1 = ss1 + ss2 = Text3.Text ss2 = ss2 + For i = 1 To nn - 1 ss1 = LTrim(ss1) hh(i) = Val(Left(ss1, InStr(ss1, ) - 1) ss1 = Mid(ss1, InStr(ss1, ) + 1) Next i For i = 1 To nn ss2 = LTrim(ss2) pp(i) = Val(Left(ss2, InStr(ss2, ) - 1) ss2 = Mid(ss2, InStr(ss2, ) + 1) Next i znxs = Val(Text4.Text) ddcshu = Val(Text5.Text) yyy = 9 Text6.Text = Dim ssAs String Dim ccf(9)As Double Dim ccz(18) As Double fname = 滤波系数施伦贝尔反演滤波系数 & .txt OpenApp.Path & fname For InputAs #1 For i = 9 To 1 Step -1 Input #1, ccf(i) ccf(i) = ccf(i) * 10 (-4) Next i For i = 0 To 18 Input #1, ccz(i) ccz(i) = ccz(i) * 10 (-4) Next i Close #1 Dim ztt() As Double ReDim ztt(nj1 To nj2)As Double For j = nj1 To nj2 For i = -9 To 18 If (j - i) = 0 Then ls1 = ppsz(j - i) Else ls1 = ppsf(Abs(j - i) End If If i = 0 Then ls2 = ccz(i) Else ls2 = ccf(Abs(i) End If ztt(j) = ztt(j) + ls1 * ls2 Next i Next j For i = nj1 To nj2 - 1 Picture1.Line (js(Exp(i * Log(10) / 6), js(ztt(i)-(js(Exp(i + 1) * Log(10) / 6), js(ztt(i + 1), RGB(0, 0, 10) Next i If ddcshu = 0 Then Dim htt() As Double ReDim htt(nj1 To nj2)As Double For i = nj1 To nj2 nmta = Exp(-1 * i * Log(10) / 6) ttt = pp(nn) For t = nn - 1 To 1 Step -1 vvv = pp(t) * (1 - Exp(-2 * nmta * hh(t) / (1 + Exp(-2 * nmta * hh(t) ttt = (vvv + ttt) / (1 + (vvv * ttt / pp(t) 2) Next t htt(i) = ttt Next i For i = nj1 To nj2 - 1 Picture1.Line (js(Exp(i * Log(10) / 6), js(htt(i)-(js(Exp(i + 1) * Log(10) / 6), js(htt(i + 1), RGB(0, 0, 255) Next i ElseIf ddcshu = 1 Then For ii = 1 To ddcshu For ee = 1 To 2 Dim llt()As Double ReDim llt(1 To nn, nj1 To nj2)As Double For i = nj1 To nj2 For j = nn To 1 Step -1 nmta = Exp(-1 * i * Log(10) / 6) ttt = pp(nn) For t = j - 1 To 1 Step -1 vvv = pp(t) * (1 - Exp(-2 * nmta * hh(t) / (1 + Exp(-2 * nmta * hh(t) ttt = (vvv + ttt) / (1 + (vvv * ttt / pp(t) 2) Next t llt(nn - j + 1, i) = ttt Next j Next i For i = nj1 To nj2 - 1 Picture1.Line (js(Exp(i * Log(10) / 6), js(llt(1, i)-(js(Exp(i + 1) * Log(10) / 6), js(llt(1, i + 1), RGB(0, 0, 255) Next i nnn = 2 * nn - 1 mmm = nj2 - nj1 + 1 Dim ppp()As Double Dim jjj()As Double Dim aaa()As Double Dim ggg() As Double Dim iii()As Double Dim bbb()As Double ReDim ppp(nnn)As Double ReDim jjj(mmm, nnn)As Double ReDim aaa(nnn, nnn)As Double ReDim ggg(mmm) As Double ReDim iii(nnn, nnn)As Double ReDim bbb(nnn)As Double For i = 1 To nnn For j = 1 To nnn If i = j Then If i Mod 2 = 1 Then iii(i, j) = 1 / (pp(i 2 + 1) 2 Else iii(i, j) = 1 / (hh(i 2) 2 End If Else iii(i, j) = 0 End If Next j Next i k = nj1 For i = 1 To mmm ggg(i) = Log(ztt(k) / llt(1, k) k = k + 1 Next i Dim ptt() As Double Dim ptp()As Double ReDim ptt(nn - 1)As Double ReDim ptp(nnn)As Double For t = nj1 To nj2 nmta = Exp(-1 * t * Log(10) / 6) For i = 1 To nn - 1 af = Exp(nmta * hh(i) bt = (1 - Exp(-2 * nmta * hh(i) / (1 + Exp(-2 * nmta * hh(i) ptt(i) = (1 - bt 2) / (1 + bt * llt(i + 1, t) / pp(i) 2 ptp(i) = (bt * (1 + 2 * bt * (llt(i + 1, t) / pp(i) + (llt(i + 1, t) / pp(i) 2) / (1 + bt * llt(i + 1, t) / pp(i) ptp(2 * i) = pp(i) * (1 - (llt(i + 1, t) / pp(i) 2) / (1 + bt * llt(i + 1, t) / pp(i) 2 * 4 * nmta / (af + 1 / af) 2 Next i ptp(i) = 1 For i = 1 To nn - 1 x = 1 For j = 1 To i - 1 x = x * ptt(j) Next j ptp(i) = x * ptp(i) ptp(2 * i) = x * ptp(2 * i) Next i x = 1 For j = 1 To i - 1 x = x * ptt(j) Next j ptp(i) = x * ptp(i) For k = 1 To nnn jjj(t - nj1 + 1, k) = (1 / llt(1, t) * ptp(k) Next k Next t For i = 1 To nnn For j = 1 To nnn For t = 1 To mmm aaa(i, j) = aaa(i, j) + jjj(t, i) * jjj(t, j) Next t Next j Next i For i = 1 To nnn For j = 1 To nnn aaa(i, j) = aaa(i, j) + znxs * iii(i, j) Next j Next i For i = 1 To nnn For j = 1 To mmm bbb(i) = bbb(i) + jjj(j, i) * ggg(j) Next j Next i ddcs = 100 Dim mAs Double m = 0 Dim y() As Double ReDim y(nnn)As Double For h = 1 To ddcs For i = 1 To nnn For j = 1 To i - 1 m = m + aaa(i, j) * ppp(j) Next j For j = i + 1 To nnn m = m + aaa(i, j) * ppp(j) Next j y(i) = (1 / aaa(i, i) * (-m + bbb(i) m = 0 Next i For i = 1 To nnn ppp(i) = y(i) Next i Next h If ee = 1 Then wucha1 = 0 For i = nj1 To nj2 wucha1 = wucha1 + (ztt(i) - llt(1, i) 2 Next i For i = 1 To nn - 1 zpp(i) = pp(i) zhh(i) = hh(i) Next i zpp(i) = pp(i) End If If ee = 2 Then wucha2 = 0 For i = nj1 To nj2 wucha2 = wucha2 + (ztt(i) - llt(1, i) 2 Next i End If If ee = 1 Then For i = 1 To nn - 1 pp(i) = pp(i) + ppp(i) hh(i) = hh(i) + ppp(2 * i) Next i pp(i) = pp(i) + ppp(i) End If Next ee If wucha1 wucha2 Then znxs = znxs / yyy Else znxs = znxs * yyy For i = 1 To nn - 1 pp(i) = zpp(i) hh(i) = zhh(i) Next i pp(i) = zpp(i) End If Text8.Text = Text9.Text = For i = 1 To nn - 1 Text8.Text = Text8.Text & hh(i) & Text9.Text = Text9.Text & pp(i) & Next i Text9.Text = Text9.Text & pp(i) ss = ss & wucha1 & vbCrLf Next ii Text6.Text = ss End If End Sub Private Sub Command2_Click() Picture1.Cls Call hkd End Sub Private Sub Command3_Click() Dim filename As String CommonDialog1.Filter = BMP 文件(*.bmp)|*.bmp|JPG 文件(*.jpg)|*.jpg CommonDialog1.FilterIndex = 1 CommonDialog1.ShowSave filename = CommonDialog1.filename If Right(CommonDialog1.filename, 4) = .bmp Or Right(CommonDialog1.filename, 4) = .jpg Then SavePicture Picture1.Image, filename End If End Sub Private Sub Command4_Click() End End Sub Private Sub Command5_Click() CommonDialog1.Filter = BMP 文件(*.bmp)|*.bmp|JPG 文件(*.jpg)|*.jpg CommonDialog1.FilterIndex = 1 CommonDialog1.ShowOpen If Right(CommonDialog1.filename, 4) = .bmp Or Right(CommonDialog1.filename, 4) = .jpg Then Picture1.Picture = LoadPicture(CommonDialog1.filename) End If Call pictiaozheng1 End Sub Private Sub Command6_Click() zbx = 3 MsgBox 请单击定位直线起始点 End Sub Private Sub Command7_Click() zbx = 2 End Sub Private Sub Command8_Click() CommonDialog1.Filter = txt 文件(*.txt)|*.txt CommonDialog1.FilterIndex = 1 CommonDialog1.ShowOpen If Right(CommonDialog1.filename, 4) = .txt Then Open CommonDialog1.filename For InputAs #1 Erase xab Erase xps Erase ppsz Erase ppsf n = 0 Do While Not EOF(1) n = n + 1 ReDim Preserve xab(n) As Double ReDim Preserve xps(n) As Double Input #1, xab(n) Input #1, xps(n) Loop Close #1 Picture1.DrawStyle = 4 For j = 1 To n - 1 Picture1.Line (js(xab(j), js(xps(j)-(js(xab(j + 1), js(xps(j + 1), RGB(0, 20, 40) Next j IfAbs(Log(xab(1) / (Log(10) / 6) - Int(Log(xab(1) / (Log(10) / 6) 0.001 Then 确定实测电阻率曲线的重采样起始点序号 nj1 = CInt(Log(xab(1) / (Log(10) / 6) Else nj1 = Int(Log(xab(1) / (Log(10) / 6) + 1 End If IfAbs(Log(xab(n) / (Log(10) / 6) - Int(Log(xab(n) / (Log(10) / 6) 0.001 Then 确定实测电阻率曲线的重采样终点序号 nj2 = CInt(Log(xab(n) / (Log(10) / 6) Else nj2 = Int(Log(xab(n) / (Log(10) / 6) End If ReDim Preserve ppsz(nj2 + 9)As Double ReDim Preserve ppsf(18 - nj1)As Double For i = nj1 To nj2 ppsz(i) = chazhi(xab(), xps(), n, i) Next i For i = 18 - nj1 To 1 Step -1 ppsf(i) = ppsz(nj1) Next i For i = 0 To nj1 - 1 ppsz(i) = ppsz(nj1) Next i For i = nj2 + 1 To nj2 + 9 ppsz(i) = ppsz(nj2) + Exp(i * Log(10) / 6) - Exp(nj2 * Log(10) / 6) Next i For i = 0 To nj2 + 9 - 1 Picture1.Line (js(Exp(i * Log(10) / 6), js(ppsz(i)-(js(Exp(i + 1) * Log(10) / 6), js(ppsz(i + 1), RGB(0, 0, 255) Next i End If End Sub Private Sub Form_Load() Call hkd End Sub Public Sub pictiaozheng1() VScroll1.Min = 0 HScroll1.Min = 0 | VScroll1.Max = (Picture1.Height - Picture2.Height) / 400 | HScroll1.Max = (Picture1.Width - Picture2.Width) / 400 | If HScroll1.Max 0 Then | HScroll1.Enabled = False | Else | HScroll1.Enabled = True | End If | If VScroll1.Max = 1 And x 10 And x 100 And x 1000 Then js = Log(x / 1000) / Log(10) * 5 + 15 End If End Function Public Function chazhi(ByRef xab() As Double, ByRef xps() As Double, ByVal nn As Integer, ByVal ii As Integer) As Double xxx = Exp(ii * Log(10) / 6) For i = 1 To nn - 1 If xab(i) = xxx Then chazhi = (xps(i + 1) - xps(i) / (xab(i + 1) - xab(i) * (xxx - xab(i) + xps(i) Exit For End If Next i If ii = nj2 Then i = nn - 1 chazhi = (xps(i + 1) - xps(i) / (xab(i + 1) - xab(i) * (xxx - xab(i) + xps(i) End If End Function 实验二、一维电测深正演实验实验二、一维电测深正演实验 一、实验目的一、实验目的 利用所讲的一维电测深的理论基础，应用计算机软件实现一维电测深曲线。 给定地电模型并合理地假定每层介质的层参数（电阻率值和层厚度），已知地 质体的形状、埋深及其与围岩的物性参数，求取该种地球物理场的剖面曲线或 分布。利用不同的测深装置，其中包括对称四极装置、二极装置等，实现不同 装置的电测深曲线。 二、实验内容二、实验内容 给定地电模型并合理地假定每层介质的层参数（电阻率值和层厚度），利用 不同的测深装置，其中包括对称四极装置、二极装置等，实现不同装置的电测 深曲线。 三、实验步骤三、实验步骤 一维电测深正演程序 Public zbx As Integer Form1 的代码： Dim x11As Double Dim x22 As Double Dim y11As Double Dim y22 As Double Private Sub Command1_Click() Picture1.DrawStyle = 0 Dim n As Integer Dim hh()As Double Dim pp()As Double Dim ccz() As Double Dim ccf()As Double Dim zzxs As Integer Dim pps(24)As Double n = Val(Text1.Text) ReDim hh(n - 1)As Double ReDim pp(n)As Double ss1 = Text2.Text ss1 = ss1 + ss2 = Text3.Text ss2 = ss2 + For i = 1 To n - 1 ss1 = LTrim(ss1) hh(i) = Val(Left(ss1, InStr(ss1, ) - 1) ss1 = Mid(ss1, InStr(ss1, ) + 1) Next i For i = 1 To n ss2 = LTrim(ss2) pp(i) = Val(Left(ss2, InStr(ss2, ) - 1) ss2 = Mid(ss2, InStr(ss2, ) + 1) Next i rr = Val(Text4.Text) gg = Val(Text5.Text) bb = Val(Text6.Text) zzxs = Combo1.ListIndex Select Case zzxs Case -1, 0 xf = 10 xz = 17 xs = 10 (-4) fname = 滤波系数对称四极梯度装置 & .txt s = 1.4678 / 1000 Case 1 xf = 6 xz = 7 xs = 1 fname = 滤波系数对称四极温纳装置 & .txt s = 0.0895 Case 2 xf = 25 xz = 10 xs = 1 fname = 滤波系数二极装置 & .txt s = -0.046339794 Case 3 xf = 1 xz = 19 xs = 1 fname = 滤波系数垂直偶极装置 & .txt s = -0.0959 Case 4 xf = 1 xz = 21 xs = 1 fname = 滤波系数径向偶极装置 & .txt s = -0.0863 End Select ReDim ccf(xf) As Double ReDim ccz(xz) As Double OpenApp.Path & fname For InputAs #1 For i = xf To 1 Step -1 Input #1, ccf(i) ccf(i) = ccf(i) * xs Next i For i = 0 To xz Input #1, ccz(i) ccz(i) = ccz(i) * xs Next i Close #1 For j = 0 To 24 ps = 0 For i = -xf To xz nmta = Exp(-1 * (j - i) * Log(10) / 6 - s) ttt = pp(n) For t = n - 1 To 1 Step -1 vvv = pp(t) * (1 - Exp(-2 * nmta * hh(t) / (1 + Exp(-2 * nmta * hh(t) ttt = (vvv + ttt) / (1 + (vvv * ttt / pp(t) 2) Next t If i 0 Then ps = ps + ttt * ccf(Abs(i) Else ps = ps + ttt * ccz(i) End If Next i pps(j) = ps Next j For j = 0 To 23 Picture1.Line (js(Exp(j * Log(10) / 6), js(pps(j)-(js(Exp(j + 1) * Log(10) / 6), js(pps(j + 1), RGB(rr, gg, bb) Next j OpenApp.Path & 成果文件 & Text7.Text & .txt For Output As #1 For i = 0 To 24 Print #1, pps(i) Next i Close #1 End Sub Private Sub Command2_Click() Picture1.Cls Call hkd End Sub Private Sub Command3_Click() Dim filename As String CommonDialog1.Filter = BMP 文件(*.bmp)|*.bmp|JPG 文件(*.jpg)|*.jpg CommonD