一种基于卫星测高重力数据的局部重力异常提取方法与流程
本发明涉及网函数插值方法，具体说是一种基于卫星测高重力数据的局部重力异常提取方法。
背景技术：
一般获得的重力数据为布格重力异常数据，它是当代重力研究中获取最为简单且研究最为广泛的重力异常数据。布格重力异常是由区域重力异常和局部重力异常叠加组合而成。其中，由分布范围较大且影响较深的地质因素(例如：区域地质构造)所引起的重力异常称为区域重力异常，由地质因素(例如：矿产，油气)所引起的小范围的重力异常称为局部重力异常。由此可见，局部重力异常的提取不仅是资源勘探的重要研究内容，同时也有助于地质构造的研究。在此基础下，如何有效地提取出局部重力异常也就成为重要问题。
目前用于提取局部重力异常的方法有很多，其中应用最为广泛的方法就是圆周法与趋势分析法。两种方法的基本原理基本一致，皆是通过计算点位的区域重力异常值，再用重力异常值减去区域重力异常值得到局部重力异常值。这两种方法的实验原理较为简易但是在实验过程中会有较多局限性。例如，会出现边缘缺失或是边缘拟合精度低的现象，由此可见，传统提取方法很难有效地提取出局部重力异常的边缘信息，而网函数插值则很好地利用边缘数据进行数据插值，为局部重力异常提取提供了新的研究思路。本发明就此提出了网函数插值在局部重力异常上提取应用。
技术实现要素：
(一)要解决的问题
实现网函数插值法在局部重力异常提取上的应用，充分考虑网函数插值法的使用以及布格重力异常的特性。通过网函数插值法计算得到区域重力异常，获取局部重力异常。
(二)技术方案
本发明包含以下步骤：
(1)利用重力数据g和sobel算子得到重力梯度数据g′，进而得到p个重力异常区域。
(2)在各重力异常区域中，提取重力异常区域边界线，得到各重力异常区域的最小外接矩形abcd和与之对应的旋转矩阵r(θ)。将各最小外接矩形abcd经过旋转变换得到变换矩形g1g2g3g4。
(3)利用步骤(2)得到的变换矩形g1g2g3g4，通过改进的曲面拟合法计算各变换矩形g1g2g3g4内的重力异常值gk。
(4)将步骤(3)得到的变换矩形g1g2g3g4进行网函数插值，得到插值g。利用旋转矩阵的逆矩阵r(θ)′，将变换矩形g1g2g3g4旋转回原最小外接矩形abcd，并从插值g中提取出重力异常区域的区域异常值gk。
(5)利用步骤(1)中的原始重力数据g与步骤(4)得到的区域重力异常值gk，计算各异常区域的局部重力异常值δgk。
进一步，所述步骤(1)中，重力梯度数据g′阈值设置为1.35mgal/km，大于重力梯度数据阈值的区域为重力异常区域，计算重力梯度数据g′公式为：
式中，gx为重力数据g的横向梯度，gy为重力数据g的纵向梯度。
进一步，所述步骤(2)中，得到各重力异常区域的最小外接矩形abcd过程如下：
a.提取重力异常区域的边界，计算边界每两相邻点的方位角θi′，并对
b.根据步骤a中两相邻点，做该异常区域的外接矩形aibicidi(i＝1，2，3...)，记录外接矩形aibicidi(i＝1，2，3...)的面积si(i＝1，2，3...)，其中面积最小的外接矩形为最小外接矩形abcd，其对应的旋转矩阵为r(θ)。
c.利用最小外接矩形abcd与相应的旋转矩阵r(θ)计算得到变换矩形g1g2g3g4，其中变换矩形中横纵坐标的最小值为x0，y0，最大值为x1，y1，变换矩形的四个角点坐标为g1(x0，y0)，g2(x1，y0)，g3(x1，y1)，g4(x0，y1)。
计算变换矩形g1g2g3g4角点坐标公式为：
式中(x1，y1)，(x2，y2)即是重力异常区域两相邻边界点坐标，θi′，θi为相邻边界点的方位角与倾斜角，r(θ)是旋转矩阵，(x，y)是重力异常区域内点的坐标，(x，y)为点(x，y)经旋转矩阵处理后的坐标。
进一步，所述步骤(3)中，在各重力异常区域中，最小外接矩形abcd内的点经过旋转矩阵r(θ)处理后分辨率可能发生改变，需要对旋转后的点进行重力异常数据插值。改进的曲面拟合法具体步骤如下：
a.当倾斜角θ＝0时，gk＝g，其中g为重力数据。
b.当倾斜角θ≠0时，设置一个a行a列的移动窗口w，使得在移动窗口w内有6个已知重力数据g。根据移动窗口w内已有的重力数据g，利用二阶曲面函数计算拟合点的拟合值g(xj，yk)t。将移动窗口w按照从左到右，从上到下的顺序移动，分别计算各拟合点的拟合值且记录各点的拟合次数n。
二阶曲面函数公式为：
g(x，y)＝a0+a1x+a2y+a3x2+a4y2+a5xy
式中g(x，y)即是每次移动窗口w拟合的二阶曲面函数，计算拟合点的拟合值g(xj，yk)，其公式为：
式中拟合点坐标为(xj，yk)，其第t次拟合结果为g(xj，yk)t(1≤t≤n)，n为点(xj，yk)的迭代曲面拟合的次数。g(xj，yk)即是点(xj，yk)经过迭代曲面拟合得到拟合结果，拟合结果g(xj，yk)即作为重力异常值gk。
进一步，所述步骤(4)中，将变换矩形g1g2g3g4内的重力异常值进行网函数插值，得到插值g，其公式为：
式中qn(n＝1，2，3，4)是过插值点q(x，y)与变换矩形g1g2g3g4平行的两条直线在变换矩形g1g2g3g4边界上截得的四个点的重力异常值，gn(n＝1，2，3，4)是变换矩形g1g2g3g4四个角点的重力异常值：an(n＝1，2，3，4，5)是这两条直线将变换矩形g1g2g3g4分成四个小矩形的面积，其中：
a1＝a5＝(x-x0)(y-y0)，a2＝(x1-x)(y-y0)，
a3＝(x1-x)(y1-y)，a4＝(x-x0)(y1-y)，
a＝a1+a2+a3+a4是变换矩形g1g2g3g4的总面积。得到了插值g后，将矩形区域内的数据旋转回各重力异常区域的最小外接矩形内。从插值g中提取出各最小外接矩形abcd格网点上值作为该点(x，y)的区域重力异常值gk。其公式为：
[xy]＝[xy]×r(θ)′
进一步，所述步骤(5)中其公式为：
δgk＝g-gk
其中δgk是各重力异常区域的原始重力数据g与区域重力异常值gk的差值。
(三)有益效果
本发明的优点体现在：
本发明采用网函数插值法来提取局部重力异常数据，与传统的圆周法、趋势分析法相比，具有边框可变、无边缘丢失，异常间水平干扰小，计算简单等优点，具有较强的可操作性与实用性，对提取局部重力异常信息具有重要意义。
附图说明
图1为本发明实施的步骤流程图
图2为坐标系旋转示意图
图3为网函数插值示意图
具体实施方式
为使本发明的目的、内容、和优点更加清楚，下面结合附图和实施例，对本发明的具体实施方式作进一步详细描述：
参照图1，本发明的具体实施步骤为：
(1)利用重力数据g和sobel算子得到重力梯度数据g′，进而得到p个重力异常区域。
其中，重力梯度数据g′阈值设置为1.35mgal/km，大于重力梯度数据阈值的区域为重力异常区域，计算重力梯度数据g′公式为：
式中，gx为重力数据g的横向梯度，gy为重力数据g的纵向梯度。
(2)在各重力异常区域中，提取重力异常区域边界线，得到各重力异常区域的最小外接矩形abcd和与之对应的旋转矩阵r(θ)，参照图2，将各最小外接矩形abcd经过旋转变换得到变换矩形g1g2g3g4。
其中，得到各重力异常区域的最小外接矩形abcd过程如下：
a.提取重力异常区域的边界，计算边界每两相邻点的方位角θi′，并对
b.根据步骤a中两相邻点，做该异常区域的外接矩形aibicidi(i＝1，2，3...)，记录外接矩形aibicidi(i＝1，2，3...)的面积si(i＝1，2，3...)，其中面积最小的外接矩形为最小外接矩形abcd，其对应的旋转矩阵为r(θ)。
c.利用最小外接矩形abcd与相应的旋转矩阵r(θ)计算得到变换矩形g1g2g3g4，其中变换矩形中横纵坐标的最小值为x0，y0，最大值为x1，y1，变换矩形的四个角点坐标为g1(x0，y0)，g2(x1，y0)，g3(x1，y1)，g4(x0，y1)。
计算变换矩形g1g2g3g4角点坐标公式为：
式中(x1，y1)，(x2，y2)即是重力异常区域两相邻边界点坐标，θi′，θi为相邻边界点的方位角与倾斜角，r(θ)是旋转矩阵，(x，y)是重力异常区域内点的坐标，(x，y)为点(x，y)经旋转矩阵处理后的坐标。
(3)利用步骤(2)得到的变换矩形g1g2g3g4，通过改进的曲面拟合法计算各变换矩形g1g2g3g4内的重力异常值gk。
其中，在各重力异常区域中，最小外接矩形abcd内的点经过旋转矩阵r(θ)处理后分辨率可能发生改变，需要对旋转后的点进行重力异常数据插值。改进的曲面拟合法具体步骤如下：
a.当倾斜角θ＝0时，gk＝g，其中g为重力数据。
b.当倾斜角θ≠0时，设置一个a行a列的移动窗口w，使得在移动窗口w内有6个已知重力数据g。根据移动窗口w内已有的重力数据g，利用二阶曲面函数计算拟合点的拟合值g(xj，yk)t。将移动窗口w按照从左到右，从上到下的顺序移动，分别计算各拟合点的拟合值且记录各点的拟合次数n。
二阶曲面函数公式为：
g(x，y)＝a0+a1x+a2y+a3x2+a4y2+a5xy
式中g(x，y)即是每次移动窗口w拟合的二阶曲面函数，计算拟合点的拟合值g(xj，yk)，其公式为：
式中拟合点坐标为(xj，yk)，其第t次拟合结果为g(xj，yk)t(1≤t≤n)，n为点(xj，yk)的迭代曲面拟合的次数。g(xj，yk)即是点(xj，yk)经过迭代曲面拟合得到拟合结果，拟合结果g(xj，yk)即作为重力异常值gk。
(4)参照图3，将步骤(3)得到的变换矩形g1g2g3g4进行网函数插值，得到插值g。利用旋转矩阵的逆矩阵r(θ)′，将变换矩形g1g2g3g4旋转回原最小外接矩形abcd，并从插值g中提取出重力异常区域的区域异常值gk。
其中，将变换矩形g1g2g3g4内的重力异常值进行网函数插值，得到插值g，其公式为：
式中qn(n＝1，2，3，4)是过插值点q(x，y)与变换矩形g1g2g3g4平行的两条直线在变换矩形g1g2g3g4边界上截得的四个点的重力异常值，gn(n＝1，2，3，4)是变换矩形g1g2g3g4四个角点的重力异常值：an(n＝1，2，3，4，5)是这两条直线将变换矩形g1g2g3g4分成四个小矩形的面积，其中：
a1＝a5＝(x-x0)(y-y0)，a2＝(x1-x)(y-y0)，
a3＝(x1-x)(y1-y)，a4＝(x-x0)(y1-y)，
a＝a1+a2+a3+a4是变换矩形g1g2g3g4的总面积。得到了插值g后，将矩形区域内的数据旋转回各重力异常区域的最小外接矩形内。从插值g中提取出各最小外接矩形abcd格网点上值作为该点(x，y)的区域重力异常值gk。其公式为：
[xy]＝[xy]×r(θ)′
(5)利用步骤(1)中的原始重力数据g与步骤(4)得到的区域重力异常值gk，计算各异常区域的局部重力异常值δgk。
其中：其公式为：
δgk＝g-gk
其中δgk是各重力异常区域的原始重力数据g与区域重力异常值gk的差值。
- 还没有人留言评论。精彩留言会获得点赞！