wgs84转bj54的高程拟合算法公式

《wgs84转bj54的高程拟合算法公式》由会员上传分享，免费在线阅读，更多相关内容在应用文档-天天文库。

1、'原理是用方程h=b0+b1*x+b2*y+b3*x*x+b4*y*y+b5*x*y来表达曲面,h指的是高程异常值，比如WGS84到bj54的高程差，然后根据6或者6个以上的公共点求出b0,b1……b5,然后如果要求某点的高程值，输入它的x,y就可以得到高程异常值h，然后利用WGS84的BLH中的H加上高程异常值就可以得到54的高程.'这个程序经过2011年01月上旬的实战精度比较高，不过存在一个弱点，就是如果北坐标比较大，如2333444.555，应该先人为的去掉最高位，这样矩阵运算才不会出异常。这是因为矩阵运算的

2、算法不够完善。有空再解决它。'CodeByKiseigo2011.01.06OptionExplicitPrivateSubcmdCalc_Click()DimmatA()AsDoubleDimmatB()AsDoubleReDimmatA(6,5)AsDouble'7个已知点ReDimmatB(6,0)AsDoubleCallSetKnownValueAB(matA,matB)DimarrPara()AsDouble'b0,b1,b2……b6这6个参数CallCalcB0toB6(matA,matB,arrPara

3、)'计算b0,b1,b2……b6这6个参数DimHoutAsDoubleHout=calcHfit(11,3,arrPara)'计算某位置的高程，这里刚好取已知点来验算FrmMain.Caption=Format(Hout,"0.000")'结果得93.7,说明结果正确EndSub'求高程拟合(二次曲面拟合)的参数B0,B1,B2,B3,B4,B5,B6ByKiseigo2011.01.0621:53HelpedbyBluePan'输入matA(5,5)最少6行，也就是最少6个已知高程点'输入matB(5,0)最少6

4、个点,这里是高程值,matB(0)是第一个点'输出：B0toB6Out,下标从0取起，一维数组，下标0-5PublicFunctionCalcB0toB6(matA()AsDouble,matB()AsDouble,B0toB6Out()AsDouble)'假设方程是h=b0+b1*x+b2*y+b3*x*x+b4*y*y+b5*x*y;方程由BluePan提供DimmaxPtAsInteger'公共点个数,要求>=6个.6表示6个点。maxPt=UBound(matA,1)+1'步骤1：加1空行,加1空列.因为矩阵

5、运算是从1开始，麻烦CallRedimMatrisAFrom1Nor0(matA)CallRedimMatrisAFrom1Nor0(matB)'步骤2：计算AT*A矩阵DimmatAT()AsDouble'A的转置矩阵ReDimmatAT(UBound(matA,2),UBound(matA,1))CallMTrans(UBound(matAT,1),UBound(matAT,2),matA,matAT)'求A的转置矩阵DimATA()AsDouble'A的转置*AReDimATA(UBound(matAT,1),

6、UBound(matA,2))'方阵CallMMul(UBound(matAT,1),UBound(matAT,2),UBound(matA,2),matAT,matA,ATA)'计算ATA(A的转置*A)'步骤3：计算(A的转置*A)的逆矩阵DimATAinv()AsDouble'A的转置*A的逆矩阵ReDimATAinv(UBound(ATA,1),UBound(ATA,2))DimiAsIntegerDimjAsIntegerFori=0ToUBound(ATA,1)Forj=0ToUBound(ATA,2)A

7、TAinv(i,j)=ATA(i,j)NextjNextiCallMRinv(UBound(ATAinv),ATAinv)'输出ATAinv'步骤4：计算ATB(A的转置*B)DimATB()AsDouble'A的转置*BReDimATB(UBound(matAT,1),UBound(matB,2))CallMMul(UBound(matAT,1),UBound(matAT,2),UBound(matB,2),matAT,matB,ATB)'计算ATB(A的转置*B)'步骤5:计算(A的转置*A的逆矩阵)*(A的转置

8、*B)DimB0toB6OutWithZero()AsDoubleReDimB0toB6OutWithZero(UBound(ATAinv,1),UBound(ATB,2))AsDoubleCallMMul(UBound(ATAinv,1),UBound(ATAinv,2),UBound(ATB,2),ATAinv,ATB,B0toB6OutW