李 杰， 鄭文霞

(湖北省地質(zhì)局 第一地質(zhì)大隊(duì),湖北 大冶 435100)

中國幅員遼闊,各地方使用的大地坐標(biāo)系各不相同,由于歷史原因造成有1954年北京坐標(biāo)系、1980西安坐標(biāo)系、地方獨(dú)立坐標(biāo)系等。根據(jù)《國土資源部國家測(cè)繪地理信息局關(guān)于加快使用2000國家大地坐標(biāo)系的通知》(國測(cè)資發(fā)[2017]30號(hào))的要求,于2018年6月底之前完成國土資源數(shù)據(jù)轉(zhuǎn)換到國家2000坐標(biāo)系,由于各坐標(biāo)系統(tǒng)使用的橢球參數(shù)不同、橢球定位定向不同,因此沒有一個(gè)嚴(yán)格的數(shù)學(xué)轉(zhuǎn)換模型進(jìn)行全面的轉(zhuǎn)換,只能進(jìn)行空間相似變換的方式求解出變換參數(shù),再進(jìn)行坐標(biāo)轉(zhuǎn)換。常用的轉(zhuǎn)換模型有布爾莎、平面四參數(shù)、三維七參數(shù)、二維七參數(shù)等模型。

病態(tài)問題是指輸入數(shù)據(jù)有微小誤差引起輸出數(shù)據(jù)相對(duì)誤差很大,病態(tài)方程組是方程組的系數(shù)矩陣和常數(shù)項(xiàng)矩陣有微小誤差,造成方程組的解不穩(wěn)定；而方程組的病態(tài)程度可以用矩陣的條件數(shù)來衡量[1]。

本文主要探討布爾莎模型及其產(chǎn)生病態(tài)方程組的原因和求解方法,并用VC++語言編程實(shí)現(xiàn)參數(shù)的求解過程。

1 布爾莎模型及其適用條件

布爾莎模型是假設(shè)了三個(gè)旋轉(zhuǎn)參數(shù)εX、εY、εZ均為微小轉(zhuǎn)角后的坐標(biāo)轉(zhuǎn)換模型,其模型[2]如下：

(1)

式中：△X0、△Y0、△Z0為三個(gè)平移參數(shù)；m為比例因子。

由于假設(shè)了三個(gè)旋轉(zhuǎn)參數(shù)為微小角度,就限制了布爾莎模型的使用范圍,一般在獨(dú)立坐標(biāo)系、地方坐標(biāo)系中等旋轉(zhuǎn)角較大的坐標(biāo)轉(zhuǎn)換時(shí),該模型就不適用了。

另外,布爾莎模型適宜進(jìn)行省級(jí)及以上大范圍坐標(biāo)轉(zhuǎn)換,對(duì)于局部小范圍的轉(zhuǎn)換容易產(chǎn)生病態(tài)方程組,使得轉(zhuǎn)換參數(shù)求解和轉(zhuǎn)換結(jié)果不穩(wěn)定。

2 布爾莎模型的誤差方程及其程序?qū)崿F(xiàn)

2.1 誤差方程組成

(2)

式中：a1=a4εX；a2=a4εY；a3=a4εZ；a4=1+m。

2.2 誤差方程程序設(shè)計(jì)的實(shí)現(xiàn)

根據(jù)公式(2)的誤差方程,可編程實(shí)現(xiàn)系數(shù)矩陣和常數(shù)項(xiàng)矩陣元素的計(jì)算。

(1) 誤差方程系數(shù)項(xiàng)計(jì)算的程序?qū)崿F(xiàn)關(guān)鍵代碼。

int k=0;

for (i=0;i

{

*(pA+k++) = 1;

*(pA+k++) = 0;

*(pA+k++) = 0;

*(pA+k++) = 0;

*(pA+k++) = 0-*(coordS_h+i);

*(pA+k++) =*(coordS_y+i);

*(pA+k++) =*(coordS_x+i);

*(pA+k++) = 0;

*(pA+k++) = 1;

*(pA+k++) = 0;

*(pA+k++) =*(coordS_h+i);

*(pA+k++) = 0;

*(pA+k++) = 0-*(coordS_x+i);

*(pA+k++) =*(coordS_y+i);

*(pA+k++) = 0;

*(pA+k++) = 0;

*(pA+k++) = 1;

*(pA+k++) = 0-*(coordS_y+i);

*(pA+k++) =*(coordS_x+i);

*(pA+k++) = 0;

*(pA+k++) =*(coordS_h+i);

}

(2) 常數(shù)項(xiàng)計(jì)算的程序?qū)崿F(xiàn)關(guān)鍵代碼

int j=0;

for (int i=0;i< PointCount; i++) 常數(shù)項(xiàng)數(shù)組計(jì)算

{

*(pL+j++)=*(coordT_x+i)-*(coordS_x+i);

*(pL+j++)=*(coordT_y+i)-*(coordS_y+i);

*(pL+j++)=*(coordT_h+i)-*(coordS_h+i);

}

3 布爾莎模型的法方程組成

3.1 法方程系數(shù)及常數(shù)項(xiàng)計(jì)算的程序?qū)崿F(xiàn)

CMatrix_cal為矩陣計(jì)算類,mat_A為誤差方程系數(shù)矩陣,mat_L為誤差方程常數(shù)項(xiàng)矩陣,mat_AT為轉(zhuǎn)置矩陣,mat_ATA為法方程系數(shù)陣,mat_ATL為法方程常數(shù)項(xiàng)矩陣。程序?qū)崿F(xiàn)如下：

CMatrix_cal mat_AT(7,PointCount*3);

mat_AT=mat_A.Transpose();

CMatrix_cal mat_ATA(7,7);

mat_ATA=mat_AT*mat_A;

CMatrix _cal mat_ATL(7,1);

mat_ATL=mat_AT*mat_L;

3.2 法方程系數(shù)矩陣產(chǎn)生病態(tài)的原因

由公式(2)可知法方程系數(shù)矩陣是一個(gè)非奇異的對(duì)稱方陣,其矩陣元素之間相差數(shù)量級(jí)很大[3],如果進(jìn)行局部小范圍坐標(biāo)轉(zhuǎn)換,公共控制點(diǎn)的縱橫坐標(biāo)差異不大,則矩陣元素行向量相關(guān)性很強(qiáng)(近似線性相關(guān)),導(dǎo)致法方程系數(shù)矩陣接近奇異陣,矩陣求逆結(jié)果很不穩(wěn)定[4],這樣,公共點(diǎn)有很小的誤差將直接導(dǎo)致法方程常數(shù)項(xiàng)的變化,系數(shù)矩陣的病態(tài)程度嚴(yán)重程度,將影響法方程解的誤差大小[5]。

要判定法方程系數(shù)矩陣的病態(tài)程度,可利用矩陣的條件數(shù)來衡量[6]。

3.3 法方程系數(shù)矩陣的條件數(shù)

矩陣的條件數(shù)的大小,決定了矩陣是“良態(tài)”還是“病態(tài)”,一般矩陣的條件數(shù)遠(yuǎn)>1時(shí)(經(jīng)驗(yàn)值為100),為中度病態(tài),條件數(shù)>1 000時(shí)為較嚴(yán)重病態(tài),條件數(shù)越大表明病態(tài)越嚴(yán)重[7]。

而矩陣的條件數(shù)與矩陣的范數(shù)有關(guān),由矩陣范數(shù)的定義,可以用矩陣的行范數(shù)、列范數(shù)和2-范數(shù)來計(jì)算,其中：

(3)

(4)

(5)

這里給出行范數(shù)的VC++關(guān)鍵程序代碼如下：

double max_sq=0;

for (int v=1;v<=7;v++)

{

for (int w=0;w<7;w++)

{

sq=sq+abs(*(ATA+qw));

qw++;

}

s[v]=sq;

if (s[v]>=max_sq)

max_sq=s[v];

sq=0;

}

同樣,也可以用列范數(shù)和2-范數(shù)編程實(shí)現(xiàn)條件數(shù)的計(jì)算,這里不再贅述。

4 布爾莎模型的病態(tài)法方程解算

對(duì)于布爾莎模型可能產(chǎn)生的病態(tài)法方程(非嚴(yán)重病態(tài)),可利用迭代改善法[8]編程求解,具體求解步驟為：

① 利用“全選主元高斯消去法”求解未知參數(shù)(即七參數(shù))的近似解X(1)；

② 代入X(1),求余向量R=B-AX(1)的值；

③ 求解AQ=R,求取Q向量的值,然后X(2)=X(1)+Q；

④ 使X(1)=X(2),重復(fù)②至④,給定一個(gè)微小量ε,

上述①中的“全選主元高斯消去法”很多參考書可找到,這里將其成員函數(shù)命名為GetRootGuass()。病態(tài)方程求解程序?qū)崿F(xiàn)的關(guān)鍵代碼[9]如下：

CMatrixTestDlg les(matrixCoef,matrixConst);

if (! les.GetRootGuass(matrixResult))

return false;

double* x = matrixResult.GetData();

q=1.0+eps;

while (q>=eps)

{

if (i==0)

return false;

i=i-1;

CMatrix_cal matrixE = matrixCoef*matrixResult;

CMatrix_cal matrixR = matrixConst - matrixE;

CMatrixTestDlg les(matrixCoef,matrixR);

CMatrix_cal matrixRR;

if (! les.GetRootGuass(matrixRR))

return false;

double * r = matrixRR.GetData();

q=0.0;

for ( k=0; k<=n-1; k++)

{

qq=fabs(r[k])/(1.0+fabs(x[k]+r[k]));

if (qq>q)

q=qq;

}

for ( k=0; k<=n-1; k++)

x[k]=x[k]+r[k];

}

5 算例

這里模擬一套原坐標(biāo)系和目標(biāo)坐標(biāo)系下公共控制點(diǎn)的坐標(biāo)數(shù)據(jù),作為算例如圖1和圖2。

圖1 原坐標(biāo)文件Fig.1 Original coordinate file

圖2 目標(biāo)坐標(biāo)文件Fig.2 Target coordinate file

程序主界面如圖3。

圖3 程序主界面
Fig.3 Main interface of program

根據(jù)程序計(jì)算得法方程系數(shù)陣見圖4。

圖4 法方程系數(shù)矩陣
Fig.4 Coefficient matrix of normal equation

行范數(shù)(條件數(shù))計(jì)算結(jié)果見圖5。由行范數(shù)(條件數(shù))為1.7×109量級(jí),可見其病態(tài)較嚴(yán)重。利用上述程序,七參數(shù)計(jì)算結(jié)果如圖6。

圖5 條件數(shù)計(jì)算Fig.5 Condition number calculation

圖6 七參數(shù)計(jì)算Fig.6 Seven parameters calculation

轉(zhuǎn)換后各控制點(diǎn)的殘差見圖7。若在A7點(diǎn)上X、Y、Z三個(gè)坐標(biāo)分量上均減去4cm(模擬誤差),重新計(jì)算參數(shù)和殘差,則計(jì)算的殘差值如圖8。

圖7 坐標(biāo)殘差計(jì)算結(jié)果
Fig.7 Calculation results of coordinate residuals

圖8 加入模擬誤差后的坐標(biāo)殘差計(jì)算結(jié)果
Fig.8 Calculation results of coodinate residual error after adding simulation error

由殘差值可見,雖然法方程系數(shù)陣病態(tài)嚴(yán)重,公共控制點(diǎn)存在較大誤差,但是利用求解病態(tài)方程組的方法,編程實(shí)現(xiàn)布爾莎模型的參數(shù)求解,求得的參數(shù)精度高(除A7點(diǎn),其余點(diǎn)基本上是毫米級(jí)轉(zhuǎn)換精度),坐標(biāo)轉(zhuǎn)換成果穩(wěn)定。

6 結(jié)語

布爾莎模型參數(shù)的求解,需要多次計(jì)算,逐步剔除有粗差的公共控制點(diǎn),求得較好的參數(shù),但要保障轉(zhuǎn)換范圍內(nèi)公共控制點(diǎn)的總體數(shù)量和點(diǎn)位分布情況,還要注意以下幾個(gè)問題：①該模型僅適合微小旋轉(zhuǎn)角的坐標(biāo)系之間的坐標(biāo)轉(zhuǎn)換；②該模型適合進(jìn)行大范圍(省級(jí)及以上)坐標(biāo)系之間的坐標(biāo)轉(zhuǎn)換；③利用病態(tài)方程組求解方法能有效地提高參數(shù)求解的精度和可靠性；④非公共控制點(diǎn)的轉(zhuǎn)換應(yīng)盡量位于公共控制點(diǎn)的范圍內(nèi)。

本文通過VC++編程實(shí)現(xiàn)了上述求解過程,在生產(chǎn)實(shí)踐中得到了驗(yàn)證,為提高成果轉(zhuǎn)換精度和工作效率,起到了很大的作用。