單張影像空間后方交會程序

1、攝影測量實習(xí)報告單張影像空間后方交會程序?qū)嵙?xí)時間 2013.5.20-2013.5.24學(xué)生班級10測繪（1）班學(xué)生姓名 房新明學(xué)生學(xué)號1072143138所在院系礦業(yè)工程學(xué)院指導(dǎo)老師張會戰(zhàn)邵亞琴一.實習(xí)目的1 .深入理解單張影像空間后方交會的原理，體會在有多余觀測情況下，用 最小二乘平差方法編程實現(xiàn)解求影像外方位元素的過程。2 .利用visual c+編寫一個完整的單張影像空間后方交會程序，通過對提 供的試驗數(shù)據(jù)進(jìn)行計算，輸出像片的外方位元素并進(jìn)行評定精度。3.通過編寫程序?qū)崿F(xiàn)單張影像空間后方交會計算，掌握非線性方程線性化 的過程、相應(yīng)數(shù)據(jù)讀入與存儲的方法以及迭代計算的特點，鞏固各類基礎(chǔ)課程

2、及 計算機(jī)課程的學(xué)習(xí)內(nèi)容，培養(yǎng)上機(jī)調(diào)試程序的動手能力，通過對實驗結(jié)果的分析， 增強(qiáng)綜合運用所學(xué)知識解決專業(yè)實際問題的能力。4、掌握空間后方交會的定義和實現(xiàn)算法(1)定義：空間后方交會是以單幅影像為基礎(chǔ)，從該影像所覆蓋地面范圍內(nèi) 若干控制點的已知地面坐標(biāo)和相應(yīng)點的像坐標(biāo)量測值出發(fā)，根據(jù)共線條件方程， 解求該影像在航空攝影時刻的外方位元素 xs, ys, zs,小，k。(2)算法：由于每一對像方和物方共腕點可列出2個方程，因此若有3個已知地面坐標(biāo)的控制點，則可列出 6個方程，解求6個外方位元素的改正數(shù) xs, ys, azs, 小， k。實際應(yīng)用中為了提高解算精度，常有多余觀 測方程，通常是在影像

3、的四個角上選取 4個或均勻地選擇更多的地面控制點，因 而要用最小二乘平差方法進(jìn)行計算。5、了解空間后方交會的基本過程(1)獲取已知數(shù)據(jù)。從攝影資料中查取影像比例尺1/m,平均攝影距離(航空攝影的航高)、內(nèi)方位元素x0, y0, f;獲取控制點的空間坐標(biāo)xt, yt, zt(2)量測控制點的像點坐標(biāo)并進(jìn)行必要的影像坐標(biāo)系統(tǒng)誤差改正，得到像點 坐標(biāo)。(3)確定未知數(shù)的初始值。單像空間后方交會必須給出待定參數(shù)的初始值， 在豎直航空攝影且地面控制點大體對稱分布的情況下，xs0和ys0為均值，zs0為航高，小、k的初值都設(shè)為00或者k的初值可在航跡圖上找出或根據(jù)控 制點坐標(biāo)通過坐標(biāo)正反變換求出。(4)計

4、算旋轉(zhuǎn)矩陣ro利用角元素近似值計算方向余弦值，組成 r陣。(5)逐點計算像點坐標(biāo)的近似值。利用未知數(shù)的近似值按共線條件式計算控 制點像點坐標(biāo)的近似值(x), (y)。(6)逐點計算誤差方程式的系數(shù)和常數(shù)項，組成誤差方程式。(7)計算法方程的系數(shù)矩陣ata與常數(shù)項atl,組成法方程式。(8)解求外方位元素。根據(jù)法方程，解求外方位元素改正數(shù)，并與相應(yīng)的近 似值求和，得到外方位元素新的近似值。(9)檢查計算是否收斂。將所求得的外方位元素的改正數(shù)與規(guī)定的限差比較， 通常對小，k的改正數(shù)4小， k給予限差，通常為0.1,當(dāng)3個 改正數(shù)均小于0.1時，迭代結(jié)束。否則用新的近似值重復(fù)(4)(8)步驟的計 算

5、，直到滿足要求為止。二、實習(xí)要求1、認(rèn)真復(fù)習(xí)單張影像空間后方交會的原理及有關(guān)內(nèi)容；2、編寫程序，計算左片號23和右片號24的外方位元素并進(jìn)行精度評定；3、驗證數(shù)據(jù)部分外方位元素和旋轉(zhuǎn)矩陣 r已經(jīng)給出，請驗證自己編寫的程序是 否正確；三、實習(xí)環(huán)境1、硬件環(huán)境：windows操作系統(tǒng)2、軟件環(huán)境：vc+ 四、實習(xí)原理1、以單幅影像為基礎(chǔ)，從該影像所覆蓋地面范圍內(nèi)若干控制點的已知地面 坐標(biāo)和相應(yīng)點的像坐標(biāo)量測值出發(fā)，根據(jù)共線條件方程，求解該影象在航空攝影時刻的像片外方位元素 xs, ys , zs，/，k共線條件方程如下： x-x0=-f*a1(x-xs)+b1(y-ys)+c1(z-zs)/a3(

6、x-xs)+b3(y-ys)+c3(z-zs) y-y0=-f*a2(x-xs)+b2(y-ys)+c2(z-zs)/a3(x-xs)+b3(y-ys)+c3(z-zs)x,y為像點的像平面坐標(biāo)；xs,ys，偽影像的外方位元素；xs, ys ,zs為攝站點的物方空間坐標(biāo)；x,y,z為物方點的物方空間坐標(biāo)；旋轉(zhuǎn)矩陣為r ;由于此共線條件方程是非線性方程，先對其進(jìn)行線性化，像點觀測值一般視為等權(quán)，即p=i;矩陣形式：v=ax-l , p=i;通過間接平差，為提高精度，增加多余觀測方程，根據(jù)最小二乘平差原理，可計 算出外方位元素的改正數(shù)。經(jīng)過迭代計算，每次迭代用未知數(shù)的近似值與上次迭 代計算的改正數(shù)

7、之和作為新的近似值，重復(fù)計算，求出新的改正數(shù)，這樣反復(fù)趨近，直到改正數(shù)小于某個限值為止。2、精度評定mi=m0*，qii 其中 m0=v【vm / (2n-6)3、公式推導(dǎo)一一 .干:0 干：0 千：一 一0fx(xs,ys,zs,q%m)=%(xs-xs)+e(ys-ys)+=hzs-zs) + xsyszs0 0 0二(：-0) ( - - 0) (1 -1 0) fx(x；,y：,z0j0,-0j0) cto干j 0mxs,yszsfixlxxs-xs) +0,fy. ys汗0(ys-ys0) 7zs(zs-z：)干0 c ：fc ：:fy (:-* y ( , - .0) y (.0)

8、 f(xo y0 z ： ,00):().：. ()() y( s, s，s，, )同理fy,xsdxs千ydys-ys千y dzs:zs干y0干：汗;0d d df；/的所以干x0dxs 二 dys，s干十.:z0-dzsffd :清生d.，小0emfy0 dxsxs以dys-ys干y0-zsdzs千y0：fy0yd- fy0 =01、基本原理共線條件線性方程式為:f dxs tdzs （ f2x- ）d4 - ?do + ydm - （ x - x計）=0fyxyy2dys 工dzs - - d -（ f 一 ）d - xd -（ y - y計）=0 zzff誤差方程式為:f2= dxs

9、+工dzs -（ f + ）d - do + yd _（ x _ x計）=vx z zf f工 dys = dzs -2；d4-（f+y ）d。- xdm-（ y - y計）=vy計算改正值：dxs,dys,dzs,d* ,ds計計算改后的外方位元素：xs xs dxs ys ys dys zs zs dzs二 d . 二 d . 二 d六、 程 序 過 程 框 圖開始七.實習(xí)數(shù)據(jù)1、模擬像片一對：左片23號 右片24號2、相片比例尺：1/300003、航攝機(jī)主距：f=150mm4、每張相片有四個控制點5、各片像點坐標(biāo)及地面坐標(biāo)片號占 八、號像點坐標(biāo)（mm）地面坐標(biāo)（m）xyxyz231-91

10、.596-74.859100000.0137500.011.003-94.23081.446100000.0142500.0036.00795.207-75.512106000.0137500.0042.00996.79783.077106000.0142500.0056.00244-102.695-79.618103000.0137500.0090.006-99.90481.754103000.0142500.0031.001086.890-77.540109000.0137500.007.001288.90476.257109000.0142500.005.00八.驗證數(shù)據(jù)已知四對點的影像

11、坐標(biāo)和地面坐標(biāo)及航攝儀內(nèi)方位元素f=153.42mm, x0=y0=0影像坐標(biāo)地面坐標(biāo)x(mm)y(mm)x (mm)y (mm)z (mm)1-86.15-68.9936589.4125273.322195.172-53.4082.2137631.0831324.51728.693-14.78-76.6339100.9724934.982386.50410.4664.4340426.5430319.81757.313、結(jié)算結(jié)果0.067530.003990.99772 -0.002210.001840.99990xs=39795.450.99771ys=27476.46 r= -0.0675

12、3zs=7572.69-0.00412九.實現(xiàn)程序輸入文件形式如下:c+源程序如下：#include #include #include #include #include using namespace std;const int n=6;void inverse (double cnn);templatevoid transpose (t1*mat1,t2*mat2,int a,int b); templatevoid multi(t1*mat1,t2 * mat2,t2 * result,int a,int b,int c);templatevoid input (t*mat,int a

13、,int b);templatevoid output(t*mat,char*s,int a,int b);int main()ofstream outfile;cout.precision(5);double x0=0.0, y0=0.0; double fk=0.15324;/內(nèi)方位元素double m=39689; /估算比例尺double b45=0.0,r33,xg61,at68,ata66,atl61;input (b,4,5);/從文件中讀取控制點的影像坐標(biāo)和地面坐標(biāo)，存入數(shù)組bdouble xs=0.0, ys=0.0, zs=0.0,q=0.0,w=0.0,k=0.0;dou

14、ble x,y,z,l81,a86;/確定未知數(shù)的出始值for(int i=0;i4;i+)xs=xs+bi2;ys=ys+bi3;zs=zs+bi4;xs=xs/4; ys=ys/4; zs=zs/4+m*fk;int f=0;do/迭代計算f+;/組成旋轉(zhuǎn)矩陣r00=cos(q)*cos(k)-sin(q)*sin(w)*sin(k);r01=-cos(q)*sin(k)-sin(q)*sin(w)*cos(k);r02=-sin(q)*cos(w);r10=cos(w)*sin(k);r11=cos(w)*cos(k);r12=-sin(w);r20=sin(q)*cos(k)+cos(

15、q)*sin(w)*sin(k);r21=-sin(q)*sin(k)+cos(q)*sin(w)*cos(k);r22=cos(q)*cos(w);/計算系數(shù)陣和常數(shù)項for(int i=0,k=0,j=0;i=6.0/206265.0|xg40=6.0/206265.0|xg50=6.0/2062 65.0);cout”迭代次數(shù)為：fendl;/精度評定double axg81,v81,vt18,vtv11,m0,d66;multi(a,xg ,axg,8,6,1);for( i=0;i8;i+)計算改正數(shù)vi0=axgi0-li0;transpose (/vt,1,8);multi(vt

16、,v,vtv,1,8,1);m0=vtv00/2;for(i=0;i6;i+)for(int j=0;j6;j+) dij=m0*ataij;/屏幕輸出誤差方程系數(shù)陣、常數(shù)項、改正數(shù)output(a,誤差方程系數(shù)陣a為：,8,6);output(l,常數(shù)項 l 為：,8,1);output(xg,改正數(shù)為：,6,1);outfile.open(aim.txt,ios:app);/打開并添加 aim.txt 文件outfile.precision(10);/以文件的形式輸出像片外方位元素、旋轉(zhuǎn)矩陣、方差陣outfile 一、像片的外方位元素為：endlendl;outfilesetw(10)xs

17、=xssetw(10)ys=yssetw(10)zs=zsendl;outfilesetw(20)航向傾角為：qsetw(10)”旁向傾角為：wsetw(10)像片旋角為：kendl;outfileendl二、旋轉(zhuǎn)矩陣 r 為:endlendl;for( i=0;i3;i+)for(int j=0;j3;j+)outfilesetw(25)rijsetw(25);outfileendl;outfileendl;outfilesetw(0)三、精度評定結(jié)果為：endl;outfile.precision(5);for(i=0;i6;i+)for(int j=0;j6;j+)outfilesetw

18、(14)dijsetw(14);outfileendl;outfile.close();return 0;templatevoid transpose(t1*mat1,t2*mat2,int a,int b) int i,j;for(i=0;ib;i+)for(j=0;ja;j+)mat2皿產(chǎn)mat1皿;return;templatevoid multi(t1*mat1,t2 * mat2,t2 * result,int a,int b,int c) int i,j,k;for(i=0;ia;i+)for(j=0;jc;j+)resultij=0;for(k=0;kb;k+)resultij+

19、=mat1ik*mat2kj;return;)template void input (t*mat,int a,int b) ifstream infile;infile.open(控制點坐標(biāo).txt);while(!infile.eof()for (int i=0;ia;i+)for(int j=0;jmatij;)infile.close();return;)templatevoid output(t*mat,char*s,int a,int b) coutsetw(15)sendl;for(int i=0;ia;i+)for(int j=0;jb;j+)coutsetw(13)matij

20、;coutendl;) return;)void inverse(double cnn) int i,j,h,k;double p;double qn12;for(i=0;in;i+)/構(gòu)造高斯矩陣for(j=0;jn;j+)qij=cij;for(i=0;in;i+)for(j=n;j12;j+)if(i+6=j)qij=1;elseqij=0;for(h=k=0;kn-1;k+,h+)/消去對角線以下的數(shù)據(jù)for(i=k+1;in;i+)if(qih=0)continue;p=qkh/qih;for(j=0;j0;k-,h-) 消去對角線以上的數(shù)據(jù) for(i=k-1;i=0;i-)if(

21、qih=0)continue;p=qkh/qih;for(j=0;j12;j+)q皿*=p;q皿產(chǎn)qk皿for(i=0;in;i+)/將對角線上數(shù)據(jù)化為1 p=1.o/qii；for(j=0;j12;j+)qij*=p;for(i=0;in;i+) / 提取逆矩陣for(j=0;jn;j+)cij=qij+6；程序的結(jié)果輸出如下：（包括文本輸出結(jié)果和熒屏輸出中間數(shù)據(jù)） cju serswa ngweideilrtopsi：.副是程，生林a n haufa n gjiaohui.exe-2.8325e-0051.8876e-00fc1.6113e_005-0.203$8-0.025079-0.06899-1.8705e-00fc-2.841te-0051.28e-005-3,05096-0.181260.0b615-2.2250-005-q.18935q.04b1870.bh221-1.5562e-00fc-2.22e-00s-1.20046-0060.01s369-0.198830.0534-2.9484e-0051.993e-00b2.9731e-006-5.154970.00307-0.07663-1.937ba-qbb-2,?s27e-005lt729w-。甌-b.02b278-q.19h620.b147