攝影后方交會0215

1、 攝影測量空間后方交會礦業(yè)工程學(xué)院_10級GIS_215 - 10 - 攝影測量單像空間后方交會10_GIS_2班 已知攝影機(jī)主距f=153.24mm，四對點的像點坐標(biāo)與相應(yīng)的地面坐標(biāo)列入下表：點號像點坐標(biāo)地面坐標(biāo)x(mm)y(mm)x(m)y(m)z(m)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.31求六個方位元素及其中誤差運用C+編寫程序，具體步驟如下圖：首先將

2、原始數(shù)據(jù)存儲到.txt文本下，供后期編程過程中調(diào)用程序代碼如下：/ 空間后方交會.cpp : Defines the entry point for the console application./#include #include #include #include #include #include#include mtxfunction.h /矩陣求逆運算、轉(zhuǎn)置、相乘 using namespace std;/*聲明結(jié)構(gòu)類型GCP，用來保存控制點坐標(biāo)*/struct GCP double x; /影像控制點x坐標(biāo) double y; /影像控制點y坐標(biāo) double Xt; /地面控制

3、點Xt坐標(biāo) double Yt; /地面控制點Yt坐標(biāo) double Zt; /地面控制點Zt坐標(biāo);/*/ 聲明讀取文本數(shù)據(jù)函數(shù)int Inputdata(int &GCPN, GCP *data, double &m, double &f, double &x0, double &y0);/*/int main() system(title 單像空間后方交會); /聲明變量 int GCPN=100; /控制點對數(shù)GCPN double m =50000, f = 153.24, x0 = 0, y0 = 0; /比例尺m，主距f，外方位元素x0,y0 GCP data4; /控制點坐標(biāo)矩陣

4、 double Xs = 0, Ys = 0, Zs = 0, fai = 0, omiga = 0, k = 0; /內(nèi)個外方位元素 double mtxX6; /X矩陣 X=Xs,Ys,Zs, mtxX3 = 1; double mtxR9; /旋轉(zhuǎn)矩陣R double mtxresult166, mtxresult26; /1保存ATA及其逆矩陣 和2保存ATL，相乘得到X結(jié)果 int i, n = 0; double sum = 0, m0; /*/ / 讀取數(shù)據(jù)gotokaishi: ; /goto標(biāo)記 if(Inputdata(GCPN, data, m, f, x0, y0) c

5、out 讀取“原始數(shù)據(jù).txt”出現(xiàn)錯誤！n; cout 請將“原始數(shù)據(jù).txt”置于程序所在目錄中，并檢查數(shù)據(jù)格式。; cout n數(shù)據(jù)格式如下所示n(第一行為控制點個數(shù)和m，第二行為f,x0,y0，其余行為控制點坐標(biāo)x,y,Xt,Yt,Zt)n; cout 4t0.0n; cout 153.24t0.0t0.0n; cout -86.15t-68.99t36589.41t25273.32t2195.17n; cout -53.40t 82.21t37631.08t31324.51t728.69n; cout -14.78t-76.63t39100.97t24934.98t2386.50n;

6、 cout 10.46t 64.43t40426.54t30319.81t757.31n; cout ch; if(Y = ch | y = ch) system(cls); goto gotokaishi; cout 程序即將退出，謝謝使用！n ; system(PAUSE); return 0; /讀取數(shù)據(jù)不正確，返回1 /計算Xs0,Ys0,Zs0 for(int i = 0; i GCPN; +i) Xs += datai.Xt / 4; Ys += datai.Yt / 4; Zs = f * m; arrayXtg; /地面坐標(biāo)歸算后值Xt arrayYtg; /地面坐標(biāo)歸算后值Y

7、t arrayZtg ; /地面坐標(biāo)歸算后值Zt arrayxtg; /影像坐標(biāo)歸算后值xt arrayytg; /影像坐標(biāo)歸算后值yt double mtxA2 * 506; /系數(shù)矩陣A 4個a11a16，a21a26 double mtxAT62 * 50; /A的轉(zhuǎn)置 double mtxL2 * 50; /L矩陣 L=lx,ly cout 主距f= f 比例尺m= m endl; cout 像點坐標(biāo)ttt 地面點坐標(biāo)n xtt ytt Xttt Yttt Ztn; for(int i = 0; i GCPN; +i) cout datai.x datai.y datai.Xt dat

8、ai.Yt datai.Zt 0.00001 | mtxX4 0.00001 | mtxX5 0.00001) & n 100) /*每次迭代錢將ATA和ATL清空*/ for(int i = 0; i 6; +i) for(int j = 0; j 6; +j) mtxresult1ij = 0; mtxresult2i = 0; +n;/統(tǒng)計迭代次數(shù) cout n* *; printf(n第%d次迭代n, n); cout Xs= Xs tYs= Ys tZs= Zs n= fai t= omiga t= k endl; /*旋轉(zhuǎn)矩陣R*/ mtxR0 = cos(fai) * cos(k

9、) - sin(fai) * sin(omiga) * sin(k); mtxR1 = -cos(fai) * sin(k) - sin(fai) * sin(omiga) * cos(k); mtxR2 = -sin(fai) * cos(omiga); mtxR3 = cos(omiga) * sin(k); mtxR4 = cos(omiga) * cos(k); mtxR5 = -sin(omiga); mtxR6 = sin(fai) * cos(k) + cos(fai) * sin(omiga) * sin(k); mtxR7 = -sin(fai) * sin(k) + cos

10、(fai) * sin(omiga) * cos(k); mtxR8 = cos(fai) * cos(omiga); printf(旋轉(zhuǎn)矩陣R:n); cout mtxR0 t mtxR1 t mtxR2 endl; cout mtxR3 t mtxR4 t mtxR5 endl; cout mtxR6 t mtxR7 t mtxR8 endl endl; /*像點坐標(biāo)計算值*/ /cout歸算得到控制點坐標(biāo)n; /coutXtgt Ytgt Ztgttendl; /*xtgt ytgt */ for(int i = 0; i GCPN; +i) /*歸算控制點坐標(biāo)*/ Xtgi = mtx

11、R0 * (datai.Xt - Xs) + mtxR3 * (datai.Yt - Ys) + mtxR6 * (datai.Zt - Zs); Ytgi = mtxR1 * (datai.Xt - Xs) + mtxR4 * (datai.Yt - Ys) + mtxR7 * (datai.Zt - Zs); Ztgi = mtxR2 * (datai.Xt - Xs) + mtxR5 * (datai.Yt - Ys) + mtxR8 * (datai.Zt - Zs); double x_x0, y_y0; /x-x0=-f(X/Z)見書38頁（2-5-8） x_x0 = -f * X

12、tgi / Ztgi; y_y0 = -f * Ytgi / Ztgi; /*計算L矩陣*/ mtxL2 * i = datai.x - (x_x0); mtxL2 * i + 1 = datai.y - (y_y0); /*計算A矩陣*/ mtxA2 * i0 = (mtxR0 * f + mtxR2 * (x_x0) / Ztgi; mtxA2 * i1 = (mtxR3 * f + mtxR5 * (x_x0) / Ztgi; mtxA2 * i2 = (mtxR6 * f + mtxR8 * (x_x0) / Ztgi; mtxA2 * i3 = (y_y0) * sin(omiga)

13、 - (x_x0) / f) * (x_x0) * cos(k) - (y_y0) * sin(k) + f * cos(k) * cos(omiga); mtxA2 * i4 = -f * sin(k) - (x_x0) / f) * (x_x0) * sin(k) + (y_y0) * cos(k); mtxA2 * i5 = (y_y0); mtxA2 * i + 10 = (mtxR1 * f + mtxR2 * (y_y0) / Ztgi; mtxA2 * i + 11 = (mtxR4 * f + mtxR5 * (y_y0) / Ztgi; mtxA2 * i + 12 = (m

14、txR7 * f + mtxR8 * (y_y0) / Ztgi; mtxA2 * i + 13 = -(x_x0) * sin(omiga) - (y_y0) / f) * (x_x0) * cos(k) - (y_y0) * sin(k) - f * sin(k) * cos(omiga); mtxA2 * i + 14 = -f * cos(k) - (y_y0) / f) * (x_x0) * sin(k) + (y_y0) * cos(k); mtxA2 * i + 15 = -(x_x0); /*輸出矩陣A*/ /*cout矩陣An; for(int i=0;i8;+i) for(

15、int j=0;j6;+j) coutmtxAij_; coutendl; */ /*解法方程*/ / 矩陣相關(guān)運算在頭文件 mtxfunction.h中 / 求A的轉(zhuǎn)置AT Mxt_ZhuanZhi(&mtxA00, &mtxAT00, 2 * GCPN, 6); / 求A與AT的乘積，保存到mtxresult1矩陣 Mtx_XiangChen(&mtxAT00, &mtxA00, &mtxresult100, 6, 2 * GCPN, 6); / 求ATA的逆矩陣 Mtx_QiuNi(&mtxresult100, 6); / 求AT與L的乘積，結(jié)果保存到mtxresult2矩陣 Mtx_X

16、iangChen(&mtxAT00, mtxL, &mtxresult20, 6, 2 * GCPN, 1); / 求ATA的逆矩陣與ATL的乘積。即X=(ATA)-1 ATL Mtx_XiangChen(*mtxresult1, &mtxresult20, &mtxX0, 6, 6, 1); /Xs=Xs+Xs mtxX即X矩陣 Ys，Zs同理 Xs += mtxX0; Ys += mtxX1; Zs += mtxX2; /=+ mtxX即X矩陣 ，同理 fai += mtxX3; omiga += mtxX4; k += mtxX5; /輸出每次迭代計算的改正數(shù)/ cout 改正數(shù)n ;

17、for(int i = 0; i 6; +i) printf( mtxX%d=%lfn, i, mtxXi); /*旋轉(zhuǎn)矩陣R*/ mtxR0 = cos(fai) * cos(k) - sin(fai) * sin(omiga) * sin(k); mtxR1 = -cos(fai) * sin(k) - sin(fai) * sin(omiga) * cos(k); mtxR2 = -sin(fai) * cos(omiga); mtxR3 = cos(omiga) * sin(k); mtxR4 = cos(omiga) * cos(k); mtxR5 = -sin(omiga); mt

18、xR6 = sin(fai) * cos(k) + cos(fai) * sin(omiga) * sin(k); mtxR7 = -sin(fai) * sin(k) + cos(fai) * sin(omiga) * cos(k); mtxR8 = cos(fai) * cos(omiga); /輸出結(jié)果到文件 FILE *fp_out; if(!(fp_out = fopen(空間后方交會計算結(jié)果.txt, w)cout 錯誤，結(jié)果未能保存到文件n; /*輸出到文件*/ fprintf(fp_out, n單像空間后方交會計算程序); fprintf(fp_out, n解算完成); fpr

19、intf(fp_out, n迭代次數(shù)為:%dn, n); fprintf(fp_out, n像片外方位元素的解n); fprintf(fp_out, 航向頃角=%lftn, fai); fprintf(fp_out, 旁向傾角=%lftn, omiga); fprintf(fp_out, 像片旋角=%lftn, k); fprintf(fp_out, tXs=t%lfntYs=t%lfntZs=t%lfn, Xs, Ys, Zs); fprintf(fp_out, n旋轉(zhuǎn)矩陣R:n); fprintf(fp_out, %lftt%lftt%lftn, mtxR0, mtxR1, mtxR2);

20、 fprintf(fp_out, %lftt%lftt%lftn, mtxR3, mtxR4, mtxR5); fprintf(fp_out, %lftt%lftt%lftn, mtxR6, mtxR7, mtxR8); /*輸出到屏幕*/ cout nn單像空間后方交會計算程序 ; cout n解算完成 ; cout n迭代次數(shù)為: n endl; printf(n像片外方位元素的解n); cout 航向頃角=t fai endl; cout 旁向傾角=t omiga endl; cout 像片旋角=t k endl; printf(tXs=t%lfntYs=t%lfntZs=t%lfn,

21、Xs, Ys, Zs); cout endl; printf(旋轉(zhuǎn)矩陣R:n); cout mtxR0 t mtxR1 t mtxR2 endl; cout mtxR3 t mtxR4 t mtxR5 endl; cout mtxR6 t mtxR7 t mtxR8 endl endl; /*計算單位權(quán)中誤差m0*/ for(i = 0; i 2 * GCPN; i+) sum += pow(mtxLi, 2); m0 = sqrt(sum / (2 * GCPN - 6); /*輸出到文件*/ fprintf(fp_out, n單位權(quán)中誤差m0=%lfn, m0); fprintf(fp_o

22、ut, n六個外方位元素精度:n); fprintf(fp_out, Xs的精度為:t%lf米n, m0 * sqrt(mtxresult100); fprintf(fp_out, Ys的精度為:t%lf米n, m0 * sqrt(mtxresult111); fprintf(fp_out, Zs的精度為:t%lf米n, m0 * sqrt(mtxresult122); fprintf(fp_out, 的精度為:t%lf秒n, m0 * sqrt(mtxresult133); fprintf(fp_out, 的精度為:t%lf秒n, m0 * sqrt(mtxresult144); fprin

23、tf(fp_out, 的精度為:t%lf秒n, m0 * sqrt(mtxresult155); fprintf(fp_out, nn解算完成 ); /*輸出到屏幕*/ cout 單位權(quán)中誤差m0= m0 endl; cout 六個外方位元素精度: endl; cout Xs的精度為:t m0 *sqrt(mtxresult100) 米 endl; cout Ys的精度為:t m0 *sqrt(mtxresult111) 米 endl; cout Zs的精度為:t m0 *sqrt(mtxresult122) 米 endl; cout 的精度為:t m0 *sqrt(mtxresult133)

24、 秒 endl; cout 的精度為:t m0 *sqrt(mtxresult144) 秒 endl; cout 的精度為:t m0 *sqrt(mtxresult155) 秒 endl; cout nn解算完成n; fclose(fp_out); fp_out = NULL; MessageBox(NULL, 解算完成，點擊確定打開結(jié)果保存文件！,空間后方交會, 0); system(空間后方交會計算結(jié)果.txt); /system(PAUSE); return 0;/*/ GCPN是控制點對數(shù)，data用來保存控制點數(shù)據(jù)，m是比例尺，f,x0,y0是內(nèi)方位元素/ 從程序所在文件夾讀取“原始

25、數(shù)據(jù).txtint Inputdata(int &GCPN, GCP *data, double &m, double &f, double &x0, double &y0) FILE *fp_data; if(!(fp_data = fopen(D:原始數(shù)據(jù).txt, r) return 1; /如果文件讀取失敗，返回1 cout 讀取到的數(shù)據(jù)為n; fscanf(fp_data, %d%lf, &GCPN, &m); cout 控制點個數(shù)： GCPN endl; cout 航攝比例尺： m endl; if(GCPN 4) return 2; /控制點低于4個（不足），返回2 fscanf(fp_data, %lf%lf%lf, &f, &x0, &y0); f /= 1000; /毫米單位換算為米 cout 主距f= f tx0= x0 ty0= y0 endl; if(m 0 | f 0) return 3; /如果m0或f0，則返回3，即文本文件格式不對 / data = new GCP GCPN; /分配控制點對數(shù)個GCP類型空間 / cout像點坐標(biāo)ttt地面點坐標(biāo)nxttyttXtttYtttZtn; for(int i = 0; i GCPN; +i) if(fscanf(fp_data, %lf%lf%lf%lf%lf, &d