重磁勘探实验,质子磁力仪,matlab日变校正

质子磁力仪实验
实验一：质子磁力仪操作及数据采集实验
试验目的：了解质子磁力仪的使用以及野外数据采集
试验步骤：
1、布置测线：
这次试验的场地是教四楼前面的小树林，测线方向是东西向，测线间距为4m共有5条，测点自西向东增加，每条测线26个测点,间距为2m。

下图为测网布置的草图：
2、组装仪器，同步各个仪器的时间。

3、找一个没有外部干扰的地方布置日变观测站。

4、建立测量基站。

5、测量数据。

6、处理数据,并做解释。

实验结果：
1、日变曲线：
由日变曲线可看出，日变站的磁场变化值大概在200nT左右，在60000秒---60400秒间磁场出现了较大的扰动，并且持续了大约两分钟，由此，对各个测点的日变校正是非常重要的，而且在处理日变之前，必须先进行曲线平滑处理以去除日变意外的干扰（本次试验对平滑处理不作要求）。

以下是个测线未经过日变校正的的剖面图。

2、各测线剖面图：
图1 测线1
由上图可知，磁场在2m附近出现一个极小值点大概在44000nT左右，并且在4m处磁场快速回升；在10m和17m处有两个极大值。

图2 测线2
如上图所示，2号测线所观测到的磁场值变化幅度大，在12m到20m处磁场变化非常剧烈，有300nT左右的变化，而且可以判断此处的磁场梯度比较大；20m 到30m处曲线相对平缓，之后出现小的波动。

图3 测线3
3号测线的场值大小总体上是由西向东减小的，主要是在1、2号测点附近磁场下降的比较快，此处梯度比较大，向东磁场变化比较平缓，波动较少。

图4 测线4
4号测线的波动比3号测线大，在20m和35m处出现极小值点，而且这两个极小值点附近的磁场变化梯度都比较大。

与3号测线一样，4号测线左端是高磁场值，右端是低场值。

图5 测线5
5号测线在15m至25m处出现一个极小值点，在该处向两侧的磁场变化很大，约在200nT左右，而且梯度也很大。

图6 剖面平面图
由剖面平面图大致可以判定五条测线的基本情况：南边的场值为正异常，北边为负异常，中间异常比较低。

2、测区的平面等值线图：
图6 平面等值线图
上图为测区未经过日变校正的平面等值线图，其中X轴为北向，Y轴为测线方向，由西向东。

由图可以发现一条磁场值为46900nT的等值线把测区分成南北两侧，在南侧17m和32m都出现了一个极大值区域，而且由等值线的疏密程度可以判定这个区域附近的场值梯度很大；在46900nT等值线的北端出现两个极小值区域与南端的极大值区域几乎是对称的。

由此可以大胆的判定，这个对称异常的之间必
定存在某种磁性体，使得测区附近的磁场出现如此对称的高低异常。

并且由异常的分布状况看，这个异常体应该不会超过1号测线和5号测线(16m),因为对称磁异常的出是因为磁性体的磁化造成磁体边缘场值与地磁场叠加产生，从而产生对称分布，此处磁化方向为X负方向，磁性体被磁化之后，自身产生磁力线，这些磁力线在磁体北侧是背离磁化方向的，也就是X轴负方向，故在北边出现了负异常；而在磁性体的南侧，磁性体的磁力线方向是向北方的，顺着地球磁力线方向即X轴正向，所以该处为正异常。

在测区的最西侧也有两个对称分布的异常，而且情况与上述完全相反，由于测区的范围比较窄只能通过梯度来判断异常的存在，但是不能确定异常体的分布状况在如果是真正进行磁测工作，应该增加测线数量并且延长测线长度，进一步对测区进行调查。

试验二：日变校正程序编制实验
实验目的：掌握日变校正过程，并实现日变校正自动处理
思路：
因为这次试验是个小范围的测区，所以使用日变观测站作为基站，也就是说把日变观测站的正常场值To作为测区的正常场值。

设Td为日变观测站的观测值，Tr为当天的日变值，To为该区的背景值；Tc为观测点的观测值，△T为异常值。

则有如下式子：
Td=Tr+To
Tc=To+△T+Tr
△T=Tc- Td
那么，我们要求日变校正之后的异常值△T，只需要用Tc- Td就行了。

而且只要Tc和Td时间上严格同步，我们就可可以直接作差，但是Td和Tc的观测时间并不是一一同步的，所以，不同步的值就需要做线性插值进行计算，如下图：
t1为前一个日变观测值为T1的记录时间，t2为后一个日变观测值T2的记录时间（本次试验t1-t2=30s），我们要求的是观测点t时刻的日变值T，也就是上述公式中的Td。

下面通过软件实现日变校正的计算，我们使用MATLA做矩阵运算进行日变校正。

首先把所有的测点数据包括记录时间是和场值，保存在单独的表格中；把日变记录数据也保存在单独的表格中（包括时间和场值）。

然后我们把这两个数
据分别倒入matlab中，生成两个i行两列的矩阵，其中测点的数据保存在C矩阵中，日变数据保存在D矩阵中，矩阵的第一列保存的是时间（秒作单位），第二列保存的是场值，矩阵的行数i由测点或日变的记录数量决定。

下面是具体的操作步骤。

步骤：
1、导出数据，由于本次试验是手工记录，需要把数据输入Excel表格中，按照
测线C1、C2、C3、C4、C5进行保存,记录的数据包括点号、测线号、时间、测量值，然后我们把点号测线号换成直角坐标，第一条测线Y为0第五条测线为16，第一个测点X为0，最后一个点为50。

为了一次性对所有数据进行处理我们可以把所有测点的数据整合起来，如下图。

以同样的方式把日变观测数据记录在D表格中。

2、首先要把各个测点的时间由时分秒换算成秒的格式：如左图所示，建立一个
文件名为T.txt的记事本,然后把时间这一列全部复制上去，如左图。

然后导入MATLAB中生成i行3列的矩阵T，i有数据的个数决定，如右图：
T(i,1)表示矩阵的第i行第1列，也就是第一个时间的小时数比如16个小时，那么T(i,2)就表示分，T(i,3)表示秒。

根据公式：秒=时*3600+分*60+秒，进行计算，新建一个m文件用来保存转换函数，方程如下：
function zhuanhuan(T)%建立转换时间的函数
for i=1:length(T)%获取矩阵列数，确定循环次数
t(i,1)=T(i,1)*3600+T(i,2)*60+T(i,3);%计算第i行时间保存于t矩阵
end
方程编辑完成后保存如图：
然后在command window中调用我们建立好的方程，输入t矩阵以显示结果，结果如图：
然后把转换后的时间数据导出，并替换原来的时间数据，这就完成了时间的转换了，下面开始进行日变校正。

用上述方式把测量值数据和日变观测数据导入MATLAB中建立以时间t为第一列，场值T为第二列的矩阵如图为测量值的矩阵C.txt日变值保存为D.txt 如下图 C.txt。

那么我们就有了两个矩阵，下面介绍各个函数变量的所表示的意义：C(i,1)表示第i个测点的记录时间，D(j,1)、D(j+1,1)分别表示两个相邻日变记录的时间，R(i,2)表示第i个测点的日变校正结果。

左图所示是MATLAB中生
成的C矩阵。

下面建立方程如上图：
function jiaozheng(C,D)%建立函数
for i=1:length(C);%获取测点数据的列数一确定循环次数
for j=1: length(D)-1;%获取日变数据列数以确定循环次数
if D(j,1)<=C(i,1) & C(i,1)<D(j+1,1)%限定i测点时间在两个相邻日变时间之间
R(i,2)=C(i,2)-((D(j+1,2)-D(j,2))*(C(i,1)-D(j,1))/(D(j+1,1)-D(j,1)));%时间差之比*日变差值得到日变改正值
R(i,1)=C(i,1);%为矩阵R附上测点对应的时间。

end
end
end
处理结果如下图所示：
我们可以进一步计算，求出日变化曲线数据，结果如图所示：
值△T：
测线1
测线1的异常值主要是正异常，而且在15m和35m处异常值相对较高。

测线2
测线2的异常是正异常，而且表现出较大的波动，异常值达到400nT的高异常。

测线3
测线3的异常值也是正异常，最高异常达到800nT,各个测点总体上维持300nT左右的正异常。

测线4
测线4的异常呈现处正负异常交替出现的现象，整个测线的异常值波动较大，而且负异常出现在中间。

测线5
测线5呈现两边正异常，中间负异常的现象，而且正负异常的落差很大。

平面图
经过日变校正之后的的平面图和原始数据呈现的样式差不多，都表现为对称异常，正负异
常都很大，梯度较大。

剖面平面图
由剖面平面图大致可以判定五条测线的基本情况：南边的场值为正异常，北边为负异常，中间异常比较低。