利用沙漠地表高程数据计算沙丘底部面高程的方法及系统

摘要

本发明提出了一种利用沙漠地表高程数据计算沙丘底部面高程的方法及系统，该方法包括：输入工区测点坐标与参数；计算地形趋势面；计算所述地形趋势面与沙丘谷底间的高程差函数；确定属于沙丘谷底测点的高程分界面；计算沙丘谷底趋势面以及沙丘底部面。与常规沙丘底部面的确定方法相比，该方法显著提高了效率，不受人为因素的影响，结果稳定可靠，得到的沙丘底部面就在沙丘谷底附近的下方，整体趋势平滑。该发明已经在中国西部沙漠区实际资料处理中得到广泛应用，取得了良好的效果。

说明书

技术领域

本发明属于地震勘探资料数据处理等技术领域，尤其涉及沙漠工区地震勘探静校正计算技术中的从沙漠地表高程数据中计算沙丘底部面高程的技术。

背景技术

沙漠工区地震勘探资料处理中，应用沙丘曲线法相关技术计算静校正量时需要沙丘底部面的高程数据。地震勘探测量阶段只能得到沙漠工区随沙丘起伏的地表高程，无法直接得到沙丘底部面的高程数据。沙丘底部面高程都是通过对地表高程数据的分析处理得到。现有技术是由技术人员人工选择一系列被认为是沙丘谷底的测量样点，利用内插外推的方式计算整个工区的沙丘底部面，效率较低，结果受人为因素及所选沙丘谷底的测量样点的影响。

发明内容

本发明提出一种依据沙漠工区测量得到的地表高程数据，快速准确地计算沙丘底部面高程的方法，达到提高沙丘底部面计算效率和避免人为因素对计算结果影响的目的。

根据本发明的一个方面，提供一种利用沙漠地表高程数据计算沙丘底部面高程的方法，该方法包括：

输入工区测点坐标与参数；

计算地形趋势面；

计算所述地形趋势面与沙丘谷底间的高程差函数；

确定属于沙丘谷底测点的高程分界面；

计算沙丘谷底趋势面以及沙丘底部面；

输出沙丘底部面。

进一步地，输入工区测点的平面与高程坐标(x_i,y_i,z_i)，i＝1,2,...,N；N为测点总数。

输入的参数包括：平滑半径参数R1、R2与R3，测点的最大高程大于谷底测点分界面高程的幅度t，以及高程差最大值剔除百分比α。

进一步地，以高程属性数据(x,y,z)为输入，采用逐点平面拟合的方法，计算地形趋势面(x,y,u)：

以测点(x_i,y_i)为中心，半径R1范围内的所有测点，以平面拟合方法构建u＝ax+by+c的平面方程，由该平面方程计算的该点高程即所述地形趋势面在该点的高程u_i；所有测点的地形趋势面高程构成地形趋势面；其中下标i是测点序号。

进一步地，计算所述地形趋势面与沙丘谷底间的高程差函数，包括：

构建地形趋势面与地表之间的高程差属性数据(x,y,u-z)，将其中z_i≥u_i即地表高程大于等于地形趋势面高程的测点设置成无效高程差，z_i<u_i即地表高程小于地形趋势面高程的测点设置为有效高程差；

以高程差属性数据(x,y,u-z)为输入，移动平均计算平均高程差分布函数(x,y,v)：以测点(x_i,y_i)为中心，半径R2范围内的除被设置为无效高程差的点以外的所有测点都参与计算平均值，得到该点的平均高程差v_i；所有测点的平均高程差构成平均高程差分布函数；其中下标i是测点序号。

进一步地，确定属于沙丘谷底测点的高程分界面，包括：

应用地形趋势面与沙丘谷底的平均高程差分布函数(x,y,v)，即重置u＝u-v得到新的高程属性(x,y,u)，为沙丘谷底测点的分界面。

进一步地，计算沙丘谷底趋势面，包括：

构建沙丘谷底地表测点高程属性数据(x,y,w)，当z_i-u_i<t时，w_i＝z_i，其它情况高程w_i设置成无效高程；

以高程属性数据(x,y,w)为输入，逐点平面拟合计算沙丘谷底趋势面(x,y,p)：

以测点(x_i,y_i)为中心，半径R3范围内的除被设置为无效高程的点以外的所有测点都参与平面拟合，构建p＝ax+by+c的平面方程，由该平面方程计算的该点高程，即沙丘谷底趋势面在该点的高程p_i；所有测点的沙丘谷底趋势面高程构成沙丘谷底趋势面；其中下标i是测点序号。

进一步地，计算沙丘底部面，包括：

从沙丘谷底趋势面高程(x,y,p)与地表高程(x,y,z)的差值p-z中提取一个高程调整量Δz，最终的沙丘底部面高程(x,y,q)等于沙丘谷底趋势面高程减去该调整量，即q_i＝p_i-Δz，其中下标i是测点序号。

进一步地，提取高程调整量的步骤包括：

沙丘谷底趋势面高程减去地表面高程得高程差，保留高程差正值部分即地表高程在沙丘谷底趋势面下方的部分；

对高程差正值数据从大到小排序，从排序后数据的最大值开始剔除α百分比的点，剩下点中的最大值即为高程调整量Δz。

根据本发明的另一方面，提供一种利用沙漠地表高程数据计算沙丘底部面高程的系统，该系统包括：

存储器，存储有计算机可执行指令；

处理器，所述处理器运行所述存储器中的计算机可执行指令，执行以下步骤：

输入工区测点坐标与参数；

计算地形趋势面；

计算所述地形趋势面与沙丘谷底间的高程差函数；

确定属于沙丘谷底测点的高程分界面；

计算沙丘谷底趋势面以及沙丘底部面；

输出沙丘底部面。

与常规沙丘底部面确定方法相比，该方法显著提高了效率，不受人为因素的影响，结果稳定可靠，得到的沙丘底部面就在沙丘谷底附近的下方，整体趋势平滑。该发明已经在中国西部沙漠区实际资料处理中得到广泛应用，取得了良好的效果。

附图说明

通过结合附图对本公开示例性实施方式进行更详细的描述，本公开的上述以及其它目的、特征和优势将变得更加明显，其中，在本公开示例性实施方式中，相同的参考标号通常代表相同部件。

图1示出了本发明实施例的测点地形高程图。

图2示出了本发明实施例的地形趋势面图。

图3示出了本发明实施例的地形趋势面与沙丘谷底间的平均高程差分布函数图。

图4示出了本发明实施例的属于沙丘谷底测点的分界面图。

图5示出了本发明实施例的沙丘谷底趋势面图。

图6示出了本发明实施例的沙丘底部面图。

图7示出了本发明实施例的沙丘厚度图。

图8示出了本发明实施例的方法流程图。

具体实施方式

下面将参照附图更详细地描述本公开的优选实施方式。虽然附图中显示了本公开的优选实施方式，然而应该理解，可以以各种形式实现本公开而不应被这里阐述的实施方式所限制。相反，提供这些实施方式是为了使本公开更加透彻和完整，并且能够将本公开的范围完整地传达给本领域的技术人员。

本发明提供了一种依据沙漠工区测量得到的地表高程数据，快速准确地计算沙丘底部面高程的方法，提高沙丘底部面计算效率，避免人为因素对计算结果的影响。采用逐步接近的方法计算沙丘底部面，结果在沙丘谷底附近的下方，整体趋势平滑。

本公开提出了一种利用沙漠地表高程数据计算沙丘底部面高程的方法，该方法包括：

输入工区测点坐标与参数；

计算地形趋势面；

计算所述地形趋势面与沙丘谷底间的高程差函数；

确定属于沙丘谷底测点的高程分界面；

计算沙丘谷底趋势面以及沙丘底部面；

输出沙丘底部面。

本发明方法的技术原理是，采用逐步接近的方法自动计算沙丘底部面，技术要点包括寻找沙丘谷底、计算沙丘底部趋势面、计算沙丘底部面等。

首先，计算地形趋势面。所有地表测点坐标数据为(x_i,y_i,z_i)，i＝1,2,...,N，N为总点数，用向量表示成(x,y,z)。用输入的地表测点数据构建地表测点的高程属性数据(x,y,z)，所有测点的高程都是有效的。以高程属性数据(x,y,z)为输入，采用逐点较大区域范围内求平均(即移动平均)或平面拟合的方法，计算地形趋势面(x,y,u)。考虑到移动平均在工区边界存在边界效应，采用逐点平面拟合方法更适合取得沙漠工区沙丘地形的趋势。就具体测点而言，以该点为中心给定区域范围内的所有测点都参与平面拟合，构建一个u＝ax+by+c的平面方程，由该平面方程计算的该点高程即地形趋势面在该点的高程。所有测点的地形趋势面的高程构成地形趋势面。给定区域范围可以是圆、矩形等，这里采用给定较大半径的圆，该较大的半径至少要大于沙丘起伏横向跨度的2倍。

计算地形趋势面与沙丘谷底间的平均高程差分布函数。构建地形趋势面与地表之间的高程差属性数据(x,y,u-z)，将其中z_i≥u_i即地表高程大于等于地形趋势面高程的测点设置成无效高程差，z_i<u_i即地表高程小于地形趋势面的测点设置为有效高程差，式中下标i是测点序号。以高程差属性数据(x,y,u-z)为输入，采用移动平均的方法，计算平均高程差分布函数(x,y,v)。就具体测点而言(即使其高程差被设置成无效)，以该点为中心给定区域范围内的除被设置为无效高程差的点以外的所有测点都参与算术平均，得到该点的平均高程差。所有测点的平均高程差构成平均高程差分布函数。实际上是计算从地形趋势面到沙丘谷底还要下移的初步量。给定区域范围可以是圆、矩形等，这里采用给定相对较小的半径的圆，该相对较小的半径一般应不大于上述“较大半径”的二分之一。

确定属于沙丘谷底测点的分界面。对地形趋势面高程(x,y,u)，应用地形趋势面与沙丘谷底的平均高程差分布函数(x,y,v)，即重置u＝u-v得新的高程属性(x,y,u)即沙丘谷底测点的分界面。高程在沙丘谷底测点分界面(x,y,u)以下的所有测点被看成是属于沙丘谷底的点。

计算沙丘谷底趋势面。构建新的地表测点高程属性数据(x,y,w)，其中w依据u与地表高程z的关系而定。当z_i-u_i<t时，w_i＝z_i，其它情况高程w_i设置成无效高程，式中下标i是测点序号。其中t为给定的门槛值(比如2)，表示参与沙丘谷底趋势面计算的测点的最大高程大于谷底测点分界面高程的幅度，保证有效数据均来自沙丘谷底及附近的点。以高程属性数据(x,y,w)为输入，采用逐点较大区域范围内求平均(即移动平均)或平面拟合的方法，计算沙丘谷底趋势面(x,y,p)。考虑到移动平均在工区边界存在边界效应，采用逐点平面拟合方法更适合取得沙漠工区沙丘谷底的趋势。就具体测点而言(即使其高程w被设置成无效)，以该点为中心给定区域范围内的除被设置为无效高程的点以外的所有测点都参与平面拟合，构建一个p＝ax+by+c的平面方程，由该平面方程计算该点的高程，即沙丘谷底趋势面在该点的高程。所有测点的沙丘谷底趋势面的高程构成沙丘谷底趋势面。给定区域范围可以是圆、矩形等，这里采用给定较大半径的圆。

计算沙丘底部面。从沙丘谷底趋势面高程(x,y,p)与地表高程(x,y,z)的差值p-z中提取一个高程调整量Δz，最终的沙丘底部面高程(x,y,q)等于沙丘谷底趋势面高程减去该调整量，即q_i＝p_i-Δz，式中下标i是测点序号。提取该调整量的步骤包括如下2步：第1步，沙丘谷底趋势面高程减去地表面高程得高程差，保留高程差正值部分即地表高程在沙丘谷底趋势面下方的部分；第2步，对高程差正值数据从大到小排序，从排序后数据的最大值开始剔除总测点数一定比例(比如α＝0.2表示剔除比例0.2％)的点，剩下点中的最大值即为高程调整量Δz。剔除一定比例的最大值，在工区存在沟槽时有必要。如果工区不存在沟槽，剔除比例可以设定为0即不剔除任何最大值。如此提取高程调整量的目的是保证除了仅有极少数测点的地表高程在沙丘底部面下方外，其它所有测点的沙丘底部面就在沙丘谷底附近的下方，而且整体趋势平滑。

上述所述平面拟合，是指针对一组成对的数据(x₁,y₁,z₁,x₂,y₂,z₂,......,x_M,y_M,z_M)，M为参与平面拟合的点数，找出一个z＝ax+by+c的平面方程，使达到最小。针对本发明技术，测点的(x,y)就是该成对数据的(x,y)，测点的各种高程值就是该成对数据中的z。

上述所述以向量表示的测点的各种高程属性数据，如(x,y,u)、(x,y,w)等等。(x,y,u)代表测点数据(x_i,y_i,u_i)，i＝1,2,...,N；(x,y,w)代表测点数据(x_i,y_i,w_i)，i＝1,2,...,N，其中N为总点数，其余类推。

根据本发明的另一实施方式，提供一种利用沙漠地表高程数据计算沙丘底部面高程的系统，该系统包括：

存储器，存储有计算机可执行指令；

处理器，所述处理器运行所述存储器中的计算机可执行指令，执行以下步骤：

输入工区测点坐标与参数；

计算地形趋势面；

计算所述地形趋势面与沙丘谷底间的高程差函数；

确定属于沙丘谷底测点的高程分界面；

计算沙丘谷底趋势面以及沙丘底部面；

输出沙丘底部面。

为便于理解本发明实施例的方案及其效果，以下给出一个具体应用示例。本领域技术人员应理解，该示例仅为了便于理解本发明，其任何具体细节并非意在以任何方式限制本发明。

如图8所示，图8示出了本发明实施例的方法流程图。

首先，输入工区测点坐标与平滑半径等参数。

沙漠工区所有测点的平面与高程坐标(x_i,y_i,z_i)，i＝1,2,...,N，3个平滑半径参数R1、R2与R3，参与沙丘谷底趋势面计算的测点的最大高程大于谷底测点分界面高程的幅度t，以及用于从沙丘谷底趋势面计算沙丘底部面所用的高程调整量时高程差最大值剔除百分比α。

接下来，计算地形趋势面。

构建地表测点的高程属性数据(x,y,z)，所有测点的高程都是有效的。以高程属性数据(x,y,z)为输入，逐点平面拟合计算地形趋势面(x,y,u)。就具体测点而言，以该点为中心半径R1范围内的所有测点都参与平面拟合，构建u＝ax+by+c的平面方程，由该平面方程计算的该点高程即地形趋势面在该点的高程。

接下来，计算地形趋势面与沙丘谷底间的平均高程差分布函数。

构建地形趋势面与地表之间的高程差属性数据(x,y,u-z)，将其中z_i≥u_i即地表高程大于等于地形趋势面高程的测点设置成无效高程差，z_i<u_i即地表高程小于地形趋势面高程的测点设置为有效高程差，式中下标i是测点序号。

以高程差属性数据(x,y,u-z)为输入，移动平均计算平均高程差分布函数(x,y,v)。就具体测点而言(即使其高程差被设置成无效)，以该点为中心半径R2范围内的除被设置为无效高程差的点以外的所有测点都参与计算平均值，得到该点的平均高程差。所有测点的平均高程差构成平均高程差分布函数。

接下来，确定属于沙丘谷底测点的分界面。

对地形趋势面高程属性数据(x,y,u)，应用地形趋势面与沙丘谷底间的平均高程差分布函数(x,y,v)，即重置u＝u-v得新的高程属性数据(x,y,u)即沙丘谷底测点的分界面。

接下来，计算沙丘谷底趋势面。

构建沙丘谷底地表测点高程属性数据(x,y,w)，其中w依据u与地表高程z的关系而定。当z_i-u_i<t时，w_i＝z_i，其它情况高程w_i设置成无效高程，式中下标i是测点序号。

以高程属性数据(x,y,w)为输入，逐点平面拟合计算沙丘谷底趋势面(x,y,p)。就具体测点而言(即使其高程w被设置成无效)，以该点为中心半径R3范围内的除被设置为无效高程的点以外的所有测点都参与平面拟合，构建p＝ax+by+c的平面方程，由该平面方程计算的该点高程，即沙丘谷底趋势面在该点的高程。

接下来，计算沙丘底部面。

从沙丘谷底趋势面高程(x,y,p)与地表高程(x,y,z)的差值p-z中提取一个高程调整量Δz，最终的沙丘底部面高程(x,y,q)等于沙丘谷底趋势面高程减去该调整量，即q_i＝p_i-Δz，式中下标i是测点序号。

提取该高程调整量的步骤包括如下2步：第1步，沙丘谷底趋势面高程减去地表面高程得高程差，保留高程差正值部分即地表高程在沙丘谷底趋势面下方的部分；第2步，对高程差正值数据从大到小排序，从排序后数据的最大值开始剔除α百分比的点，剩下点中的最大值即为高程调整量Δz。

最后，输出沙丘底部面数据。

实施例

本实施例为新疆沙漠三维工区的实际资料。图1是测点地表高程。输入参数R1＝6000m、R2＝3000m、R3＝6000m、t＝2、α＝0。

图2是执行“计算地形趋势面”步骤得到的地形趋势面。以图1所示的测点高程属性数据为输入，逐点平面拟合计算地形趋势面。就具体测点而言，以该点为中心半径R1＝6000m范围内的所有测点都参与平面拟合，构建一个平面方程，计算该点的高程即地形趋势面在该点的高程。

图3是执行“计算地形趋势面与沙丘谷底间的平均高程差分布函数”步骤得到的地形趋势面与沙丘谷底间的平均高程差分布函数。以地形趋势面高程减去地表高程的高程差中的正值部分,即地表高程小于地形趋势面高程的测点的高程差数据为输入，采用移动平均的方法计算所有测点的平均高程差。就具体测点而言(即使其地表高程不小于地形趋势面高程)，以该点为中心半径R2＝3000m范围内所有地表高程小于地形趋势面高程的测点都参与计算平均值，得到该点的平均高程差。

图4是执行“确定属于沙丘谷底测点的分界面”步骤得到的属于沙丘谷底测点的分界面，是图2的地形趋势面减去图3的地形趋势面与沙丘谷底间的平均高程差后的结果。

图5是执行“计算沙丘谷底趋势面”步骤得到的沙丘谷底趋势面。以高程不高于沙丘谷底测点分界面t米(t＝2)的所有测点的高程数据为输入，逐点平面拟合计算沙丘谷底趋势面。就具体测点而言(即使其高程大于等于沙丘谷底测点分界面加上2米)，以该点为中心半径R3＝6000m范围内所有高程不高于沙丘谷底测点分界面2米的测点都参与平面拟合，构建一个平面方程，由该平面方程计算的该点高程即沙丘谷底趋势面在该点的高程。

图6是执行“计算沙丘底部面”步骤得到的沙丘底部面。从图5所示沙丘谷底趋势面高程减去图1所示地表高程的高程差中提取一个高程调整量。由于该工区不存在明显的沟槽，在提取高程差时并没有在排序后的高程差数据中剔除一定比例的最大值(α＝0)，而是直接取高程差的最大值为高程调整量。沙丘底部面是图5所示的沙丘谷底趋势面减去高程调整量后的结果。

图7是地表面到沙丘底部面的高程差，即沙丘厚度图。

本发明提供了一种依据沙漠工区测量得到的地表高程数据，快速准确地计算沙丘底部面高程的方法，提高沙丘底部面计算效率，避免人为因素对计算结果的影响。采用逐步接近的方法计算沙丘底部面，结果在沙丘谷底附近的下方，整体趋势平滑。

以上已经描述了本公开的各实施例，上述说明是示例性的，并非穷尽性的，并且也不限于所披露的各实施例。在不偏离所说明的各实施例的范围和精神的情况下，对于本技术领域的普通技术人员来说许多修改和变更都是显而易见的。本文中所用术语的选择，旨在最好地解释各实施例的原理、实际应用或对市场中的技术的改进，或者使本技术领域的其它普通技术人员能理解本文披露的各实施例。