CN104035138B - 一种全球及局部海洋扰动重力的精确快速计算方法 - Google Patents

Abstract

本发明公开了一种全球及局部海洋扰动重力的精确快速计算方法，属于大地测量技术领域。该方法对于全球海洋扰动重力计算，通过将观测数据补充至全球区域可在每一纬度圈进行FFT计算，而后在经度方向上进行求和运算，这种处理方式保证了FFT计算后的结果与原解析算法是一致的。对于局部海洋扰动重力计算，通过对核函数和观测数据补充一定量数据从而使得核函数构成了循环矩阵，进而利用FFT计算后保证了结果与原解析算法是一致的。该发明为海洋测高卫星数据计算扰动重力提供了计算公式，同时给出了精确快速方法，本发明也可以应用于物理大地测量中大部分积分问题的求解如求解大地水准面、垂线偏差、地形改正等。

Description

一种全球及局部海洋扰动重力的精确快速计算方法

技术领域

本发明涉及精确快速计算全球及局部海洋扰动重力的方法，属于大地测量技术领域。

背景技术

传统的海洋重力场确定一般是由卫星测高数据计算获得海洋区域的重力异常，随着海洋学、空间大地测量学的发展，扰动重力较重力异常表现出更多优势。首先海洋区域由于海面地形的量级较小，扰动重力的确定可以避免重力异常计算过程中的归算问题，因此扰动重力的确定要比重力异常精确，其次，从应用角度看，相同精度情况下，利用扰动重力计算重力场扰动场元的准确性要比重力异常高。基于上述原因，越来越多的学者将扰动重力作为地球重力场的基础数据。Molodensky(1960)给出了利用垂线偏差数据计算重力异常的公式，C.Hwang(1998)利用不同的推导方法也给出了利用垂线偏差数据计算重力异常的公式，然而到目前为止，没有公开文献给出利用垂线偏差数据计算扰动重力的公式。

目前，随着卫星测高技术发展，观测数据的数量和精度大幅提升，计算方法的精确性和快速性成为制约重力场应用的主要问题，因此这也是国内外学者一直关注的重点。目前反演1′分辨率的全球海洋扰动重力，若按照严格解析算法计算，则单机全球计算一遍约需280天时间。为了解决计算慢的问题，国外的Colombo、Sideris、Forsberg等学者较早研究分析了快速傅里叶变换(FFT)算法的应用，国内学者也对此开展了深入细致的研究，特别是在2维FFT算法应用上，李建成、黄谟涛等开展了详细的研究和分析。FFT算法可提高计算的速度，但是直接用于重力场反演时将产生混叠、边缘效应等问题，这些问题使得重力场反演的精度有所下降。有学者提出在实际计算中采用补零方式解决混叠、边缘效应等问题，王冰(2004)等学者指出任意长度的函数序列补零前后的傅里叶变换结果是不相同的，补零必须使得补零后的函数序列为补零前函数序列的整数倍且为2的整数次幂。目前海洋重力场计算中的一维、二维FFT算法在应用时都很难满足上述条件，因此这些不利因素使得FFT算法在实际应用时带有很大的局限性。

发明内容

本发明的目的是：为了解决海洋扰动重力的精确快速计算问题，本发明首先给出了利用垂线偏差数据计算扰动重力的公式，同时提出了一种全球及局部海洋扰动重力的精确快速计算方法，该方法在提高计算速度同时保证计算精度等同于原解析表达式。

本发明的技术方案是：

全球海洋扰动重力的精确快速计算方法，包括如下步骤：

步骤一：采用卫星测高技术获得海洋区域的大地水准面高数据H_海，采用重力场模型获得陆地区域的大地水准面高H_陆。

式中，GM表示引力常数，γ表示计算点的正常重力，a表示参考椭球长半轴，r表示地心向径，N_max表示模型最高阶数，λ、θ分别表示计算点的经度和极距， 表示重力场模型的位系数。

将H_海和H_陆组合形成全球大地水准面高数据H，H的维数是(N+1)×(M+1)。

步骤二：根据H按照下两式得到全球N×M格网的垂线偏差南北分量ξ_ij和东西分量η_ij，i＝1,2,…N，j＝1,2,…M。M＝360°/Δλ，表示南北方向格网边长，Δλ表示东西方向格网边长。

${\mathrm{\η}}_{ij}=-\frac{H_{i,j+1}-H_{i,j}}{{\mathrm{ds}}_{\mathrm{\λ}}}---\left(3\right)$

式中，ds_λ分别表示南北方向和东西方向的格网边长。

步骤三：k＝1；

步骤四：按照下式计算第1至第N个纬度圈的核函数序列IV_ξ,ij、IV_η,ij，i＝1,2,…N，j＝1,2,…M。

IV_ξ,ij＝K'·cosα (4)

IV_η,ij＝K'·sinα (5)

式中各符号含义如下：

$K^{\′}=\frac{-cos\frac{\mathrm{\ψ}}{2}}{2{\sin }^2\frac{\mathrm{\ψ}}{2}}-\frac{cos\frac{\mathrm{\ψ}}{2}(1+2\sin \frac{\mathrm{\ψ}}{2})}{2\sin \frac{\mathrm{\ψ}}{2}(1+sin\frac{\mathrm{\ψ}}{2})}+3sin\mathrm{\ψ}$

式中， 表示被计算纬度圈中心点的纬度，表示第i个纬度圈中心点的纬度。

该步骤是本发明的核心之一，核函数K'是确定扰动重力的关键部分，由本发明推导获得。

步骤五：将每个纬度圈的IV_ξ,ij、IV_η,ij序列进行FFT变换，将ξ_ij、η_ij与需要计算的第i个纬度圈的相乘后进行FFT变换，而后按照下式计算得到第k个纬度圈的扰动重力δg_kj，j＝1,2,…M。

式中，FFT表示快速傅里叶变换，IFFT表示逆FFT计算，表示南北方向格网边长。

现有技术中计算第k个纬度圈的δg_kj采用的解析公式为：

${\mathrm{\δg}}_{kj}=\frac{\mathrm{\γ}}{4{\mathrm{\πR}}^2}\underset{\mathrm{\σ}}{\∫\∫}({\mathrm{IV}}_{\mathrm{\ξ}}\·\mathrm{\ξ}cos\mathrm{\α}+{\mathrm{IV}}_{\mathrm{\η}}\·\mathrm{\η}sin\mathrm{\α})d\mathrm{\σ}---\left(7\right)$

如果采用式(7)按照传统算法进行第k个纬度圈的扰动重力δg_kj，会产生混叠、边缘效应、计算效率低的问题，而使用本发明中提出的式(6)进行δg_kj计算，则能避免上述问题。主要原因在于步骤一中大地水准面高数据H覆盖全球，使得IV_ξ和IV_η组成了循环矩阵，而IV_ξ,ij、IV_η,ij正是循环矩阵的第一列数值，而根据循环矩阵的特性可知，式(6)与式(7)的数值计算是严格等价的。因此，步骤四中进行δg_kj的计算也是是本发明的核心之一。

步骤六：判断是否满足k<N，如果满足，则k＝k+1，进入步骤四；否则，全球第1至第N个纬度圈的扰动重力δg_kj计算全部完成，进入步骤七；

步骤七：在全球第1至第N个纬度圈的扰动重力δg_kj中，去除陆地扰动重力数据，得到全球海洋区域扰动重力Δg_global。

该步骤需要说明的是，使用全球区域大地水准面高数据计算扰动重力不仅保证了步骤四过程的严密性，而且避免了现有算法在海洋、陆地交界区域存在边缘效应的问题，从而进一步提高了近海区域海洋扰动重力计算的准确性。

局部扰动重力的精确快速方法，包括如下步骤：

步骤一：由卫星测高技术获得局部海洋区域的大地水准面高数据H，维数是(N+1)×(M+1)，按照下式计算得到局部海洋区域的垂线偏差南北分量ξ_ij和东西分量η_ij，i＝1,2,…N，j＝1,2,…M。

${\mathrm{\&eta;}}_{ij}=-\frac{H_{i,j+1}-H_{i,j}}{{\mathrm{ds}}_{\mathrm{\&lambda;}}}---\left(9\right)$

步骤二：k＝1；

步骤三：计算第1至第N个纬度圈的核函数序列IV_ξ,ij、IV_η,ij，i＝1,2,…N，j＝1,2,…M。

IV_ξ,ij＝K'·cosα (10)

IV_η,ij＝K'·sinα (11)

式中符号含义见上，在局部区域，将核函数序列IV_ξ,ij、IV_η,ij进行改造，每一个序列增加M-2个元素，然后变成新的序列具体如下。

${\mathrm{IV}}_{\mathrm{\&xi;},i,j}^{new}={\mathrm{IV}}_{\mathrm{\&xi;},i,0},{\mathrm{IV}}_{\mathrm{\&xi;},i,1},{\mathrm{IV}}_{\mathrm{\&xi;},i,2},{\mathrm{IV}}_{\mathrm{\&xi;},i,3}...\mathrm{..}{\mathrm{IV}}_{\mathrm{\&xi;},i,M-1},{\mathrm{IV}}_{\mathrm{\&xi;},i,M-2},{\mathrm{IV}}_{\mathrm{\&xi;},i,M-3}\mathrm{.......}{\mathrm{IV}}_{\mathrm{\&xi;},i,1}---\left(12\right)$

${\mathrm{IV}}_{\mathrm{\&eta;},i,j}^{new}={\mathrm{IV}}_{\mathrm{\&eta;},i,0},{\mathrm{IV}}_{\mathrm{\&eta;},i,1},{\mathrm{IV}}_{\mathrm{\&eta;},i,2},{\mathrm{IV}}_{\mathrm{\&eta;},i,3}...\mathrm{..}{\mathrm{IV}}_{\mathrm{\&eta;},i,M-1},{\mathrm{IV}}_{\mathrm{\&eta;},i,M-2},{\mathrm{IV}}_{\mathrm{\&eta;},i,M-3}\mathrm{.......}{\mathrm{IV}}_{\mathrm{\&eta;},i,1}---\left(13\right)$

步骤四：将每一个纬度圈的垂线偏差南北分量ξ_ij和东西分量η_ij增加M-2个零值，形成新的序列

步骤五：将每个纬度圈的序列进行FFT变换，将与需要计算的第i个纬度圈的相乘后进行FFT变换，而后按照下式计算得到第k个纬度圈的扰动重力

该步骤是本发明的核心所在，在局部区域，式(7)中的核函数IV_ξ和IV_η所组成的矩阵并不是循环矩阵，因此局部区域在使用FFT算法时是不准确的，会产生混叠、边缘效应等问题。通过在核函数和观测数据序列中增加元素，使之满足循环矩阵的特性，从而很好地解决了上述存在的问题。该发明并不要求观测数据序列长度是2的整数次幂，对于任意长度序列都可实现精确快速计算。

步骤六：判断是否满足k<N，如果满足，则k＝k+1，进入步骤三，否则，则第1至第N个纬度圈的扰动重力均计算完毕，进入步骤七；

步骤七：得到局部海洋区域的扰动重力对于每一个纬度圈取前M个数据得到所需要的局部海洋扰动重力δg_local。

本发明的有益效果是：本发明深入研究了海洋扰动重力确定理论，给出了由垂线偏差计算扰动重力的公式。对于全球海洋扰动重力计算，通过将观测数据补充至全球区域可在每一纬度圈进行FFT计算，而后在经度方向上进行求和运算，这种处理方式保证了FFT计算后的结果与原解析算法是一致的。对于局部海洋扰动重力计算，通过对核函数和观测数据补充一定量数据从而使得核函数构成了循环矩阵，进而利用FFT计算后保证了结果与原解析算法是一致的。通过本发明，海洋扰动重力的计算速度相比较解析方法提高约20倍，且避免了海洋重力场反演中存在的混叠、边缘效应问题，计算精度等同于原解析算法，更加重要的是对观测数据的序列长度没有硬性要求，使得该发明的应用更加灵活。该发明为海洋测高卫星数据计算扰动重力提供了计算公式，同时给出了精确快速方法，本发明也可以应用于物理大地测量中大部分积分问题的求解如求解大地水准面、垂线偏差、地形改正等。

附图说明

图1是本实施例中得到的全球海洋扰动重力图；

图2是本实施例中得到的局部海洋扰动重力图。

具体实施例

实施例一：该实施例是对全球海洋扰动重力进行计算。

算例数据包括2.5′分辨率的全球海洋区域格网大地水准面高数据，该数据由DNSC08模型得到。陆地部分采用由EGM2008模型生成2.5′分辨率的大地水准面高数据，经度范围是0°-360°，纬度范围是-90°-90°。使用的常参数如下表。

地心引力常数GM	3986004.418×108m3s-2
		长半轴a	6378137.0m
正常重力γ	9.8m/s2
		地球平均半径R	6371000.0m
π	3.141592653589793

包括如下步骤：

步骤一：将海域大地水准面高数据和陆地大地水准面高数据组合形成2.5′分辨率的全球格网大地水准面高数据，格网数是4321×8641。

步骤二：按照式(2)、式(3)形成全球2.5′分辨率的垂线偏差南北分量ξ和东西分量η，格网数分别是4320×8640和4320×8640。

步骤三：k＝1；

步骤四：计算得到第k个纬度圈格网中心点处的纬度值，Δλ＝j·0.0208333°,j＝0,1…8639。将每个纬度圈的IV_ξ,ij、IV_η,ij序列进行FFT变换，将ξ_ij、η_ij与第i个纬度圈的相乘后进行FFT变换，而后按照式(6)计算得到第1个纬度圈的扰动重力δg_1j。

步骤五：计算得到第2个纬度圈格网中心点处的纬度值,Δλ＝j·0.0208333°,j＝0,1…8639。重复步骤三过程得到第二个纬度圈的扰动重力δg_2j。

步骤六：判断是否满足k<4320，如果满足，则k＝k+1，进入步骤四；否则，全球第1至第4320个纬度圈的扰动重力δg_kj计算全部完成，进入步骤七。

步骤七：得到全球范围的扰动重力，经度范围是0°-360°，纬度范围是-90°-90°，去除陆地区域扰动重力数据后得到海洋区域扰动重力数据。

本实施例中得到的扰动重力数据如下表。

本实施例中得到的海洋扰动重力见图1。

实施例二：该实施例是对局部海洋扰动重力进行计算。

算例数据包括2.5′分辨率的局部海洋区域格网大地水准面高数据，格网数是241×241，该数据由DNSC08模型得到，数据范围是东经110°-东经120°，纬度范围是北纬10°-北纬20°。使用的常参数如下表。

地心引力常数GM	3986004.418×108m3s-2
		长半轴a	6378137.0m
正常重力γ	9.8m/s2
		地球平均半径R	6371000.0m
π	3.141592653589793

包括如下步骤：

步骤一：按照式(8)、式(9)由大地水准面高数据形成局部区域2.5′分辨率的垂线偏差南北分量ξ和东西分量η，格网数分别是240×240和240×240。

步骤二：按照式(12)、式(13)对每一个纬度圈的IV_ξ,ij、IV_η,ij增加238个元素，对每一个纬度圈的ξ_ij、η_ij增加238个零值。

步骤三：k＝1；

步骤四：计算得到第1个纬度圈格网中心点处的纬度值Δλ＝j·0.0208333°,j＝0,1…239。将每个纬度圈的序列进行FFT变换，将与第i个纬度圈的相乘后进行FFT变换，而后按照式(14)计算得到第1个纬度圈的扰动重力

步骤五：计算得到第2个纬度圈格网中心点处的纬度值,Δλ＝j·0.0208333°,j＝0,1…239。重复步骤三过程得到第二个纬度圈的扰动重力

步骤六：判断是否满足k<240，如果满足，则k＝k+1，进入步骤四，否则，则第1至第240个纬度圈的扰动重力均计算完毕，进入步骤七；

步骤七：得到新的局部范围的扰动重力数据，将每一个纬度圈的扰动重力数据取前240个元素得到所求海洋区域的扰动重力数据，其经度范围是110°-120°，纬度范围是10°-20°。

本实施例中得到的扰动重力数据如下表。

序号	经度(°)	纬度(°)	扰动重力(mGal)
				1	110.020833	19.979167	-22.607613
2	110.062500	19.979167	-23.597222
				3	110.104167	19.979167	-24.207804
4	110.145833	19.979167	-23.696602
				5	110.187500	19.979167	-22.722102
6	110.229167	19.979167	-22.487323
				7	110.270833	19.979167	-23.147792
8	110.312500	19.979167	-23.727742
				9	110.354167	19.979167	-23.462396
10	110.395833	19.979167	-22.496458
				11	110.437500	19.979167	-21.222454
12	110.479167	19.979167	-19.712754
				13	110.520833	19.979167	-18.100934
……	……	……	……

本实施例中得到的海洋扰动重力见图2。

Claims (2)

1.全球海洋扰动重力的计算方法，其特征在于包括如下步骤：

步骤一：采用卫星测高技术获得海洋区域的大地水准面高数据H_海，采用重力场模型获得陆地区域的大地水准面高H_陆：

式中，GM表示引力常数，γ表示计算点的正常重力，a表示参考椭球长半轴，r表示地心向径，N_max表示模型最高阶数，λ、θ分别表示计算点的经度和极距， 表示重力场模型的位系数；

将H_海和H_陆组合形成全球大地水准面高数据H，H的维数是(N+1)×(M+1)；

步骤二：根据H按照下两式得到全球N×M格网的垂线偏差南北分量ξ_ij和东西分量η_ij，i＝1,2,…N，j＝1,2,…M，M＝360°/Δλ，表示南北方向格网边长，Δλ表示东西方向格网边长；

${\mathrm{\&eta;}}_{ij}=-\frac{H_{i,j+1}-H_{i,j}}{{\mathrm{ds}}_{\mathrm{\&lambda;}}}---\left(3\right)$

式中，ds_λ分别表示南北方向和东西方向的格网边长；

步骤三：k＝1；

步骤四：按照下式计算第1至第N个纬度圈的核函数序列IV_ξ,ij、IV_η,ij，i＝1,2,…N，j＝1,2,…M；

IV_ξ,ij＝K'·cosα (4)

IV_η,ij＝K'·sinα (5)

式中各符号含义如下：

$K^{\&prime;}=\frac{-cos\frac{\mathrm{\&psi;}}{2}}{2{\sin }^2\frac{\mathrm{\&psi;}}{2}}-\frac{cos\frac{\mathrm{\&psi;}}{2}(1+2\sin \frac{\mathrm{\&psi;}}{2})}{2\sin \frac{\mathrm{\&psi;}}{2}(1+\sin \frac{\mathrm{\&psi;}}{2})}+3\sin \mathrm{\&psi;}$

式中， 表示被计算纬度圈中心点的纬度，表示第i个纬度圈中心点的纬度；

步骤五：将每个纬度圈的IV_ξ,ij、IV_η,ij序列进行FFT变换，将ξ_ij、η_ij与需要计算的第i个纬度圈的相乘后进行FFT变换，而后按照下式计算得到第k个纬度圈的扰动重力δg_kj，j＝1,2,…M；

式中，FFT表示快速傅里叶变换，IFFT表示逆FFT计算，表示南北方向格网边长；

步骤六：判断是否满足k<N，如果满足，则k＝k+1，进入步骤四；否则，全球第1至第N个纬度圈的扰动重力δg_kj计算全部完成，进入步骤七；

步骤七：在全球第1至第N个纬度圈的扰动重力δg_kj中，去除陆地扰动重力数据，得到全球海洋区域扰动重力Δg_global。

2.局部海洋扰动重力的计算方法，其特征在于包括如下步骤：

步骤一：由卫星测高技术获得局部海洋区域的大地水准面高数据H，维数是(N+1)×(M+1)，按照下式计算得到局部海洋区域的垂线偏差南北分量ξ_ij和东西分量η_ij，i＝1,2,…N，j＝1,2,…M，M＝360°/Δλ，表示南北方向格网边长，Δλ表示东西方向格网边长；

${\mathrm{\&eta;}}_{ij}=-\frac{H_{i,j+1}-H_{i,j}}{{\mathrm{ds}}_{\mathrm{\&lambda;}}}---\left(8\right)$

式中，ds_λ分别表示南北方向和东西方向的格网边长；

步骤二：k＝1；

步骤三：计算第1至第N个纬度圈的核函数序列IV_ξ,ij、IV_η,ij，i＝1,2,…N，j＝1,2,…M；

IV_ξ,ij＝K'·cosα (9)

IV_η,ij＝K'·sinα (10)

式中各符号含义如下：

$K^{\&prime;}=\frac{-cos\frac{\mathrm{\&psi;}}{2}}{2{\sin }^2\frac{\mathrm{\&psi;}}{2}}-\frac{cos\frac{\mathrm{\&psi;}}{2}(1+2\sin \frac{\mathrm{\&psi;}}{2})}{2\sin \frac{\mathrm{\&psi;}}{2}(1+\sin \frac{\mathrm{\&psi;}}{2})}+3\sin \mathrm{\&psi;}$

式中， 表示被计算纬度圈中心点的纬度，表示第i个纬度圈中心点的纬度；

在局部区域，将核函数序列IV_ξ,ij、IV_η,ij进行改造，每一个序列增加M-2个元素，然后变成新的序列具体如下：

${\mathrm{IV}}_{\mathrm{\&xi;},i,j}^{new}={\mathrm{IV}}_{\mathrm{\&xi;},i,0},{\mathrm{IV}}_{\mathrm{\&xi;},i,1},{\mathrm{IV}}_{\mathrm{\&xi;},i,2},{\mathrm{IV}}_{\mathrm{\&xi;},i,3}...\mathrm{..}{\mathrm{IV}}_{\mathrm{\&xi;},i,M-1},{\mathrm{IV}}_{\mathrm{\&xi;},i,M-2},{\mathrm{IV}}_{\mathrm{\&xi;},i,M-3}.......{\mathrm{IV}}_{\mathrm{\&xi;},i,1}---\left(11\right)$

${\mathrm{IV}}_{\mathrm{\&eta;},i,j}^{new}={\mathrm{IV}}_{\mathrm{\&eta;},i,0},{\mathrm{IV}}_{\mathrm{\&eta;},i,1},{\mathrm{IV}}_{\mathrm{\&eta;},i,2},{\mathrm{IV}}_{\mathrm{\&eta;},i,3}...\mathrm{..}{\mathrm{IV}}_{\mathrm{\&eta;},i,M-1},{\mathrm{IV}}_{\mathrm{\&eta;},i,M-2},{\mathrm{IV}}_{\mathrm{\&eta;},i,M-3}.......{\mathrm{IV}}_{\mathrm{\&eta;},i,1}---\left(12\right)$

步骤四：将每一个纬度圈的垂线偏差南北分量ξ_ij和东西分量η_ij增加M-2个零值，形成新的序列i＝1,2,…N，j＝1,2,…2M-2；

步骤五：将每个纬度圈的序列进行FFT变换，将与需要计算的第i个纬度圈的相乘后进行FFT变换，而后按照下式计算得到第k个纬度圈的扰动重力j＝1,2,…2M-2；

式中，FFT表示快速傅里叶变换，IFFT表示逆FFT计算，表示南北方向格网边长；

步骤六：判断是否满足k<N，如果满足，则k＝k+1，进入步骤三，否则，则第1至第N个纬度圈的扰动重力均计算完毕，进入步骤七；

步骤七：得到局部海洋区域的扰动重力对于每一个纬度圈取前M个数据得到所需要的局部海洋扰动重力δg_local。