基于VC++的GPS高程曲面拟合方法的实现

摘要：GPS定位测量技术可以得到高精度的平面位置和大地高差,但在实际应用中,地面点的高程常采用的是正常高高程系统,因此,需将GPS大地高转换成正常高。利用VC++6.0实现GPS高程多项式曲面拟合程序的实现。

关键词：GPS高程，高程拟合，VC++实现

一 前言

GPS以其全天候、高精度、自动化、高效率、多功能的显著特点已经广泛的应用于大地测量、工程测量、航空摄影测量等,给测绘领域带来一场深刻的革命。大量的实践数据表明，GPS测量的平面坐标精度是可靠的，能达到工程测量的要求，其平面相对定位的精度能够达到(0.1~1)10-6D mm甚至更高(D是基线长度)。对于G P S高程,为了满足工程的需要,须将GPS中的大地高转换为工程中直接应用的正常高,才能真正实现GPS的三维测量的目的。目前,多项式曲面拟合GPS水准高程,以其计算简单、容易编程实现,已经广泛应用于工程实践。

二GPS高程多项式曲面拟合原理

当测区GPS点布设成一定区域面时,用几何水准联测部分GPS点的正常高,用数值拟合的方法求出测区似大地水准面,计算出未联测几何水准GPS点的高程异常,即:根据测区已知点的平面坐标x,y(或大地坐标B,L)和ζ值,用数值拟合原理拟合出测区似大地水准面,再内插出待求点的ζ值,从而求出待求点的正常高。

则其高程异常一般表示为ξ=f（x，y），根据测区的实际情况，f（x，y）对应有平面拟合、二次曲面拟合、三次曲面拟合3种常用数学模型。

一般来说，根据测区实际情况，确定适宜的高程异常拟合数学模型，并根据已知联测控制点的高程异常，用最小二乘拟合法确定拟合数学模型的系数。

假设某GPS测区内有i个高程联测点，异常拟合数学模型采用多项式曲面拟合：

f(x,y)= …….（1）

式中-------多项式系数，i=0,1,…,5；（x,y）-------GPS高程控制点的平面坐标。

根据n个GPS点的可以确定多项式(1)中的系数。当n≥3且n

通过最小二乘法可以求解出多项式的系数：

x=- （4）

其中P为权阵，它可以根据水准高程和GPS所测得的大地高程的精度来确定。计算出多项式系数即拟合数学模型系数后，即可以计算出任意拟合点的高程异常。

三 多项式曲面拟合程序实现

本次多项式曲面拟合程序使用VC++6.0编程，VC++6.0功能强大，编程效率高，是科学计算编程的首选编程语言。

1. 程序界面如图1所示

图1

2. 程序的编写

定义已知控制点类CKnowdate，并定义数组储存数据，

class CKnowdate : public CObject

{

public:

double x,y,h,H;

CString namepoint;

public:

CKnowdate();

virtual ~CKnowdate();

};

typedef CTypedPtrArray KNOWDATE;

定义待求数据类CDATEGPS，并定义数组储存数据。

class CDATEGPS : public CObject

{

public:

BOOL ReaddategpsCASS(CString str);//读取数据

BOOL ReadDateGps(CString StrLine);

double x,y,h,v,Hv;

CString pro;

CString Pointname ;

CDATEGPS();

virtual ~CDATEGPS();

};

typedef CTypedPtrArray DATEGPS;

定义矩阵类CMatrix，用于矩阵计算。

class CMatrix

{

CMatrix();

CMatrix(int nRows, int nCols); // 指定行列构造函数

BOOLInit(int nRows, int nCols); // 初始化矩阵

BOOL SetElement(int nRow, int nCol, double value); // 设置指定元素的值

double GetElement(int nRow, int nCol) const; // 获取指定元素的值

voidSetData(double value[]); // 设置矩阵的值

int GetNumColumns() const; // 获取矩阵的列数

int GetNumRows() const; // 获取矩阵的行数

// 数学操作

CMatrix& operator=(const CMatrix& other);

BOOL operator==(const CMatrix& other) const;

BOOL operator!=(const CMatrix& other) const;

CMatrix operator+(const CMatrix& other) const;

CMatrix operator-(const CMatrix& other) const;

CMatrix operator*(double value) const;

CMatrix operator*(const CMatrix& other) const;

// 复矩阵乘法

BOOL CMul(const CMatrix& AR, const CMatrix& AI, const CMatrix& BR, const CMatrix& BI, CMatrix& CR, CMatrix& CI) const;

// 矩阵的转置

CMatrix Transpose() const;

// 对称正定矩阵的求逆，返回值：BOOL型，求逆是否成功

// //////////////////////////////////////////////////////////////////////

BOOL CMatrix::InvertSsgj()

{

int i, j ,k, m;

double w, g, *pTmp;

pTmp = new double[m_nNumColumns];

// 逐列处理

for (k=0; k

{

w=m_pData[0];

if (w == 0.0)

{

delete[] pTmp;

return FALSE;

}

m=m_nNumColumns-k-1;

for (i=1; i

{

g=m_pData[i*m_nNumColumns];

pTmp[i]=g/w;

if (i

pTmp[i]=-pTmp[i];

for (j=1; j

m_pData[(i-1)*m_nNumColumns+j-1]=m_pData[i*m_nNumColumns+j]+g*pTmp[j];

}

m_pData[m_nNumColumns*m_nNumColumns-1]=1.0/w;

for (i=1; i

m_pData[(m_nNumColumns-1)*m_nNumColumns+i-1]=pTmp[i];

}

// 行列调整

for (i=0; i

for (j=i+1; j

m_pData[i*m_nNumColumns+j]=m_pData[j*m_nNumColumns+i];

delete[] pTmp;

return TRUE;

}

定义计算类CCAL，用于管理操作数据，当已知点数大于3且小于等于6时，可采用一次多项式曲面拟合；当已知点数大于6且小于等于 9时，可采用二次多项式曲面拟合；当已知点数大于9时，可采用三次多项式曲面拟合。

class CCAL : public CObject

{

public:

doublea1,a2,a3,a4,a5,a6,a7,a8,a9;

KNOWDATE m_knowdate;

CMatrix mtxA, mtxB, mtxC;

CMatrix A;

public:

BOOL Calx0y0();

double X0;

double Y0;

BOOL CAL3();///三次多项式

BOOL CAL2C(); ///二次多项式

BOOL CAL1C();///一次多项式

CCAL();

virtual ~CCAL();

};

四 结束语

经大量计算表明，不一定三次多项式曲面拟合精度高于二次多项式曲面拟合，需根据具体情况具体分析。如果测区面积较大或地形起伏较大,似大地水准面的形状很复杂,这样把该测区视为一个整体进行拟合,存在很多局限且精度很低。为此,可进行分区拟合,以提高拟合精度。

注：文章内的图表及公式请以PDF格式查看