高程擬合上機(jī)報(bào)告

1、地藻新號(hào)易技未號(hào)比地藻新號(hào)易技未號(hào)比toeGPS水準(zhǔn)擬合上機(jī)實(shí)驗(yàn)報(bào)告班級(jí)：測(cè)繪1202姓名：*學(xué)號(hào):*一、實(shí)驗(yàn)?zāi)康纳钊肜斫釭PS水準(zhǔn)原理利用MATLAB實(shí)現(xiàn)幾何內(nèi)插法擬合GPS水準(zhǔn)高程二、實(shí)驗(yàn)設(shè)備50個(gè)GPS水準(zhǔn)數(shù)據(jù)，選擇10個(gè)為已知點(diǎn)，其余點(diǎn)為未知點(diǎn);MATLAB 軟件。三、實(shí)驗(yàn)原理曲面擬合法是用于GPS的分布在一定區(qū)域的時(shí)候，且可以選 擇數(shù)學(xué)曲面擬合該區(qū)域的似大地水準(zhǔn)面，構(gòu)造適當(dāng)?shù)臄?shù)學(xué)模型， 計(jì)算該區(qū)域內(nèi)的高程異常值，然后求出正常高。具體思想是：已知測(cè)區(qū)的若干已知水準(zhǔn)點(diǎn)，并且用GPS測(cè)定 這些點(diǎn)的高程，利用公式求得這些點(diǎn)的高程異常，已知點(diǎn)的平面 坐標(biāo)是已知的，所以利用其平面坐標(biāo)的（x,y）

2、和高程異常值E 構(gòu)造出來(lái)的數(shù)學(xué)模型擬合最為接近該測(cè)區(qū)的似大地水準(zhǔn)面，然后 內(nèi)插求出未知點(diǎn)的高程異常值，進(jìn)而求出正常高。詳細(xì)內(nèi)容如下：測(cè)站點(diǎn)的大地高H和正常高h(yuǎn)之間有如下關(guān)系：h = H Y 多項(xiàng)式函數(shù)擬合的基本思想是在小區(qū)域GPS網(wǎng)內(nèi)，將大地水 準(zhǔn)面看成曲面（或平面），將高程異常表示為平面坐標(biāo)（dB,dL） 的函數(shù)，通過(guò)網(wǎng)中起算點(diǎn)（既進(jìn)行了 GPS測(cè)量又進(jìn)行了幾何水 準(zhǔn)聯(lián)測(cè)的點(diǎn)）已知點(diǎn)的高程異常確定測(cè)區(qū)的似大地水準(zhǔn)面形狀，求出其余各點(diǎn)的高程異常，然后根據(jù)式求出其他點(diǎn)的正常高，其數(shù)學(xué)模型為： TOC o 1-5 h z E = f(dB , dL) + 式中f(x , y)是擬合的似大地水準(zhǔn)面，

3、是擬合誤差f(dB,dL)= a + adB + a dL + a dB * dB + a * dL * dL + a dB * dL +.012345-，IdB=B - B 0,dL=L-L L 一 L 0 n其中：n為GPS網(wǎng)中點(diǎn)的數(shù)量，(B,L)為已知點(diǎn)的大地坐標(biāo)，a ,a ,a ,a ,a ,a .為擬合待定參數(shù)；x,y為各GPS點(diǎn)的平面坐標(biāo)近 012345似值，一般取平面坐標(biāo)減去全部點(diǎn)平面坐標(biāo)的均值。本次試驗(yàn)采用二次多項(xiàng)式擬合，故：取式中的一、二次項(xiàng)將大地水準(zhǔn)面擬合為:f(dB,dL) = a + adB + a dL + a dB * dB + a dL * dL + a dB *

4、 dL012345即得二次曲面擬合模型:dB& = a a a a a& = a a a a adLa5 dB * dBdL * dLdB * dL每一個(gè)起算點(diǎn)可以組成一個(gè)上式，若共存在m個(gè)這樣的點(diǎn)， 則可列出m個(gè)方程: = a + a dB + a dL + a dB2 + a dL + a dB dL +8 TOC o 1-5 h z 0112131415111E = a + a dB + a dL + a dB 2 + a dL + a B L +801222324252 22E = a + a B + a L + a dB2 + a dL + a dB dL +8 m 01 m 2 m

5、 3 m 4 m 5 m m m從而組成誤差方程：V= Bx - L 上式中，B二dB1dBdL dB上式中，B二dB1dBdL dB 2 dL2dL1 dB dL；222dB dL1dB dL2，dBmdLmdB 2mdL2 dB dLmmaaaa a123】45E2.E m2解X = ( BTPB)BtPL 解算出a即可求出網(wǎng)中其余點(diǎn)的高程異常，并利用式求出各 i未知點(diǎn)的正常高h(yuǎn)。四、實(shí)驗(yàn)步驟1.打開(kāi)MATLAB，單擊菜單欄New Function，進(jìn)入代碼編輯界面a.b.1Ifur.ction 妃13此3心 =Un.titledC jjffluL健) 3SUNnrLED Suiuvary

6、 af this iunrtian goes here用D at ad Led EtsplsnatLan gaas haraend2.程序如下%選定10個(gè)已知點(diǎn) a=4160726 469370 43.085 15.395;47377247377247362247471447435647613447623347756547691847781441548964159270414873541571774153413415794941515564155573415702747.677 19.758;37.002 9.137;56.859 28.900;43.904 15.993;78.045 50.

7、071;41.697 13.721;51.519 23.492;33.57 5.574;32.715 4.685;%求經(jīng)度B和緯度L的平均值SumB=0;SumL=0;for j=1:10SumB = SumB+a(j,1);SumL = SumL+a(j,2); endAB=SumB/10;AL=SumL/10;%確定誤差方程的系數(shù)矩陣xx=zeros(10,6);for i=1:10 x(i,1)=1.0;x(i,2)=a(i,1)-AB;x(i,3)=a(i,2)-AL;x(i,4)=(a(i,1)-AB)*(a(i,1)-AB);x(i,5)=(a(i,2)-AL)*(a(i,2)-A

8、L);x(i,6)=(a(i,1)-AB)*(a(i,2)-AL);endx %系數(shù)%確定誤差方程中的ll=zeros(10,1);for i=1:10l(i,1)=a(i,3)-a(i,4);end這些值存放在矩陣A中%求出未知數(shù) a0、a1、a2、a3、a4、a5、這些值存放在矩陣A中s=x*x;t=inv(s);u=x*l;A=t*u;%確定其他未知點(diǎn)的正常高format long b=4149569 479940 35.852 7.734;4152803 47981935.0997.007;4155546 47979533.2925.211;4159369 47911332.4114.

9、358;4157612 47986132.8574.746;4147678 48216640.73712.527;4155194 48215434.0265.838;4159612 48260333.7725.596;4151539 48229135.6897.470;4157232 48348533.4185.190;4150583 48520265.07636.737;4158508 48535539.00910.710;4160739 48563632.0663.779;4156273 486046105.53777.212;4154803 487445128.872100.495;41

10、57419 48853232.0453.628;4163806 48844737.03 8.681;4163116 49093034.7986.400;4146497 4869364153030 4875664143589 4881174146497 4869364153030 4875664143589 4881174148280 4879694150278 4881664156103 4904124146715 4905034149128 4908234154168 4910234155429 4919354144712 4923384148524 4945234150861 494172

11、4154013 4938284146430 4951244145664 4983574148161 4974194143343 5007314145375 5009434144379 5040354142039 5040854139017 504544 y=zeros(40,6); for i=1:40y(i,1)=1.0;yy49.149 20.761;84.576 56.159;44.983 16.582;58.323 29.881;86.432 57.977;36.14 7.665;59.754 31.249;88.67 60.139;100.893 72.35276.64 48.093

12、;49.964 21.421;39.018 10.372;52.246 23.604;71.39 42.776;33.804 5.158;32.367 3.633;32.36 3.630;64.272 35.468;58.95 30.142;31.911 3.044;37.269 8.389;37.935 9.068;=b(i,1)-AB;=b(i,2)-AL;y(i,4)=(b(i,1)-AB)*(b(i,1)-AL);y(i,5)=(b(i,2)-AL)*(b(i,2)-AL);y(i,6) = (b(i,1)-AB)*(b(i,2)-AL);endN=y*A;v=zeros(40,1);

13、for i=1:40v(i,1)=N(i,1)-b(i,4)-b(i,3);end2.結(jié)果如下v =8.0098565291196953.7551045948584730.130847681703031-4.944092169033343-2.63772906451608010.4596484575754620.572133290109377-5.2397394792212575.368969245265220-2.1083662852810196.590176539675184-3.781851677697794-6.696553673154273-0.8374135878394921.09

14、6102518306822-2.325603346542781-10.645678926353902-9.64424167725983511.9332740584981423.38278363983892015.7341627839806589.5779861520663236.962823616087121-0.57774405243056011.6196123472440088.4678810622394991.9062710807889260.29690344537365814.2236314372493819.2548537611357846.2241651914837142.15359498581366911.97646057913740912.9714370736124429.73766821150840715.94116584967581213.34952601232330914.66122957654646817.62035291755520221.474294417506947Figure 1File Edit View Insert T