扇形分区的地形校正方法

1、目 录1、 基本原理部分2 1、近区理论2 2、中区理论32、 地形改正值表的制作4 1、制表理论4 2、数值表及用法53、 量板的使用方法64、 编程部分7 1、近区编程及图形8 2、中区编程及图形10五、误差分析与结论12地形校正量板法的计算机实现一、基本原理部分：在山区或丘陵区施工时，各测点周围的地形不同，地形对各点重力勘探值的影响也各异，这往往使重力布格异常曲线产生假异常或曲线明显畸变，以致影响异常的分辨和解释。为了提高数据的精度，通常需要进行改正，借此来消除测点所在平面上、下多余或缺失部分物质的影响，使测点处于上面无物质、下面完全被物质充填的平面。根据万有引力定律，质量元dm对A点产

2、生引力的垂直分量为： （1）式中：G-万有引力常数； -dm至A点的距离； 与轴间的夹角。为了达到以上校正的目的，通过以下来两种方法进行改进。1 近区地改020米。分为05米，510米和1020米三个环，每环分为8个方位； 近区扇形域划分 扇形柱体地形校正计算图通常近区第一环的地改值是以扇形锥公式来计算的，即： 式中：n测点周围所等分的扇形锥数目； R近区地改半径； i扇形锥上下侧面中心线的夹角。在相同地形条件下，如果把测点近区的地形抽象为一个漏斗且测点处于其底端，或抽象为一个圆锥且测点处于锥顶，那么测点周围几个扇形锥的角i均为近区的最大地形倾角，即。由上式计算并求和，即可得到测点近区的地改值

3、为： 由上式推导不难看出，当测点近区的地形倾角时，其近区地改值。为了更确切地反映近区的地形影响，在理论计算时将近区划分为05m、510m、1020m 3个环，采用扇形域地形改正计算公式推导如下：令式中： 第环内半径； 第环外半径； 扇形柱体相对于测点A的平均高差； n测点周围整个圆环分成的等分数，这里n=8。 所以可得：A点的总地改值 2 中区地改20700米，分为2050米、50100米、100200米、200300米、300500米、500700米六环。其中，前三环分为8个扇形柱，后三环分为16个扇形柱；用于改正的公式为扇形公式： 3远区地改值变化平缓，将其作为区域重力场背景看待，在异常划

4、分时予以消除，不会对局部重力异常明显影响。 2、 地形改正值表的制作：A、原理：利用的公式为： 制作的思路：为加快计算速度取=1g/cm3和给出不同的高度（h）时，按上式可将其制成g(T)-h表格，计算地形影响值时只需查表格即可。当工区地表岩石平均密度为时，只需将结果乘以即可。地形改正值表半径为0-20米，角度为1-14度时，的关系RiR(i+1)2*pi/n(度)数值转成弧度高程（h）2*6.67*pi/ndg(T)*106(mGal)02010.017453293-3500.1164134612.2618014602020.034906585-3300.2328269224.5155605

5、5602030.052359878-3100.3492403836.75972526402040.06981317-2900.4656538448.99231725602050.087266463-2700.58206730611.2107749702060.104719755-2500.69848076713.4117219302070.122173048-2300.81489422815.5906153502080.13962634-2100.93130768917.7411962502090.157079633-1901.0477211519.85459685020100.1745329

6、25-1701.16413461121.91782857020110.191986218-1501.28054807223.91108573020120.20943951-1301.39696153325.80262776020130.226892803-1101.51337499427.53828113020140.244346095-901.62978845629.01765889半径为20-200米，角度为1-14度时，的关系205010.017453293-3500.1164134613.145209098205020.034906585-3300.2328269226.2488696

7、42205030.052359878-3100.3492403839.303107931205040.06981317-2900.46565384412.29794489205050.087266463-2700.58206730615.2205496205060.104719755-2500.69848076718.05415512205070.122173048-2300.81489422820.77644682205080.13962634-2100.93130768923.35710338205090.157079633-1901.0477211525.753933742050100.

8、174532925-1701.16413461127.906607212050110.191986218-1501.28054807229.726100632050120.20943951-1301.39696153331.076221792050130.226892803-1101.51337499431.739921652050140.244346095-901.62978845631.355733425010010.017453293-3500.1164134614.6039074685010020.034906585-3300.2328269229.068048065010030.05

9、2359878-3100.34924038313.367671555010040.06981317-2900.46565384417.472052325010050.087266463-2700.58206730621.342663865010060.104719755-2500.69848076724.930689135010070.122173048-2300.81489422828.17360045010080.13962634-2100.93130768930.990440955010090.157079633-1901.0477211533.2753319450100100.1745

10、32925-1701.16413461134.8886745250100110.191986218-1501.28054807235.6457449350100120.20943951-1301.39696153335.3035761650100130.226892803-1101.51337499433.5511789750100140.244346095-901.62978845630.0208026710020010.017453293-3500.1164134617.08871914210020020.034906585-3300.23282692213.723555781002003

11、0.052359878-3100.34924038319.8412592910020040.06981317-2900.46565384425.3683250210020050.087266463-2700.58206730630.2195658610020060.104719755-2500.69848076734.296822510020070.122173048-2300.81489422837.4881585510020080.13962634-2100.93130768939.6681820910020090.157079633-1901.0477211540.70063603100

12、200100.174532925-1701.16413461140.44521944100200110.191986218-1501.28054807238.77187592100200120.20943951-1301.39696153335.58757527100200130.226892803-1101.51337499430.88260096100200140.244346095-901.62978845624.8039436 用法：读出每个扇形柱的平均高程，将每一个扇形柱的高程减去测点高ho得到高差，然后根据高差以及方位角的度数从地改表的相应环中查出它的校正值。例如0-20米的环中高

13、程h=-90（米），方位角为14度，对应的改正值为2.91*10-（10）再将各环各区校正值相加便得到一个测点的总的地形校正值。三、中区地形改正量板的使用方法 地形改正区域划分图 量板使用示意图附表1三、编程部分近区地改020米。分为05米，510米和1020米三个环，每环分为8个方位，其具体编程为：clear allclose allt=-10*rand(40,4);e=10*rand(40,4);x=t,e;m,n=size(x);G=6.67*10-6;%单位为mGalrow=2.67*103;R1=5;R2=10;R3=20;N=n;M=m/4;%计算的点数for k=1:1:M fo

14、r i=1:1:4 for j=1:1:8 B(i,j,k)=abs(x(i+4*(k-1),j);%数据转换成每个点所包含的各方位各环地形值 end endendfor k=1:1:M for j=1:1:N D(k,1,j)=(2*pi*G*row/N)*R1*(1-R1/sqrt(R12+B(1,j,k)2);%0-5（锥形公式） D(k,2,j)=(2*pi*G*row/N)*R2*(1-R2/sqrt(R22+B(2,j,k)2);% 0-10（锥形公式） D(k,3,j)=(2*pi*G*row/N)*(sqrt(R12+B(3,j,k)2)-R1-sqrt(R22+B(3,j,k

15、)2)+R2);%5-10（扇形公式） D(k,4,j)=(2*pi*G*row/N)*(sqrt(R22+B(4,j,k)2)-R2-sqrt(R32+B(4,j,k)2)+R3);%10-20（扇形公式） end gg(k)=(sum(D(k,1,:)+sum(D(k,3,:)+sum(D(k,2,:)/2+sum(D(k,4,:)%各点总地改值J=var(gg)%均方差plot(1:length(gg),gg,'r:*')xlabel('点号')ylabel('g （ 1 0 -6 m / s2 ）')grid onend附表2中区地改20

16、700米，分为2050米、50100米、100200米、200300米、300500米、500700米六环，前三环分为8个扇形柱，后三环分为16个扇形柱，其具体编程如下：3.1 前三环分为8个扇形clear allclose allt=-10*rand(40,4);e=10*rand(40,4);x=t,e;m,n=size(x);G=6.67*10-6;%单位为mGalrow=2.67*103;R1=20;R2=50;R3=100;R4=200N=n;M=m/3;%计算的点数for k=1:1:M for i=1:1:3 for j=1:1:8 B(i,j,k)=abs(x(i+3*(k-1

17、),j);%数据转换成每个点所包含的各方位各环地形值 end endendfor k=1:1:M for j=1:1:N D(k,1,j)=(2*pi*G*row/N)*(sqrt(R12+B(1,j,k)2)-R1-sqrt(R22+B(1,j,k)2)+R2;% 20-50（扇形公式） D(k,2,j)=(2*pi*G*row/N)*(sqrt(R22+B(2,j,k)2)-R2-sqrt(R32+B(2,j,k)2)+R3);%50-100（扇形公式） D(k,3,j)=(2*pi*G*row/N)*(sqrt(R32+B(3,j,k)2)-R3-sqrt(R42+B(3,j,k)2)+

18、R4);%100-200（扇形公式） End gg(k)=sum(D(k,1,:)+sum(D(k,2,:)+sum(D(k,3,:)%各点总地改值J=var(gg)%均方差plot(1:length(gg),gg,'r:*')xlabel('点号')ylabel('g （ 1 0 -6 m / s2 ）')grid onend3.2、后三环分为16个扇形clear allclose allt=-10*rand(60,8);e=10*rand(60,8);x=t,e;m,n=size(x);G=6.67*10-6;%单位为mGalrow=2.67

19、*103;R1=200;R2=300;R3=500;R4=700；N=n;M=m/3;%计算的点数for k=1:1:M for i=1:1:3 for j=1:1:16 B(i,j,k)=abs(x(i+3*(k-1),j);%数据转换成每个点所包含的各方位各环地形值 end endendfor k=1:1:M for j=1:1:N D(k,1,j)=(2*pi*G*row/N)*(sqrt(R12+B(1,j,k)2)-R1-sqrt(R22+B(1,j,k)2)+R2);%200-300（扇形公式） D(k,2,j)=(2*pi*G*row/N)*(sqrt(R22+B(2,j,k)2

20、)-R3-sqrt(R22+B(2,j,k)2)+R3);%300-500（扇形公式） D(k,3,j)=(2*pi*G*row/N)*(sqrt(R32+B(3,j,k)2)-R3-sqrt(R42+B(4,j,k)2)+R4);%500-700（扇形公式） end gg(k)=sum(D(k,1,:)+sum(D(k,3,:)+sum(D(k,2,:)+sum(D(k,4,:)%各点总地改值J=var(gg)%均方差plot(1:length(gg),gg,'r:*')xlabel('点号')ylabel('g （ 1 0 -6 m / s2 ）')grid onend4、 误差与结果分析1、地形校正的计算分为近区、中区、远区分别计算（1）近区