布爾莎模型病態(tài)方程組解算的程序設計與實現(xiàn)

李 杰， 鄭文霞

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

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

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

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

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

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

(1)

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

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

另外,布爾莎模型適宜進行省級及以上大范圍坐標轉(zhuǎn)換,對于局部小范圍的轉(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 誤差方程程序設計的實現(xiàn)

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

(1) 誤差方程系數(shù)項計算的程序?qū)崿F(xià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ù)項計算的程序?qū)崿F(xiàn)關鍵代碼

int j=0;

for (int i=0;i< PointCount; i++) 常數(shù)項數(shù)組計算

{

*(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ù)項計算的程序?qū)崿F(xiàn)

CMatrix_cal為矩陣計算類,mat_A為誤差方程系數(shù)矩陣,mat_L為誤差方程常數(shù)項矩陣,mat_AT為轉(zhuǎn)置矩陣,mat_ATA為法方程系數(shù)陣,mat_ATL為法方程常數(shù)項矩陣。程序?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ù)矩陣是一個非奇異的對稱方陣,其矩陣元素之間相差數(shù)量級很大[3],如果進行局部小范圍坐標轉(zhuǎn)換,公共控制點的縱橫坐標差異不大,則矩陣元素行向量相關性很強(近似線性相關),導致法方程系數(shù)矩陣接近奇異陣,矩陣求逆結(jié)果很不穩(wěn)定[4],這樣,公共點有很小的誤差將直接導致法方程常數(shù)項的變化,系數(shù)矩陣的病態(tài)程度嚴重程度,將影響法方程解的誤差大小[5]。

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

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

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

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

(3)

(4)

(5)

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

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ù)編程實現(xiàn)條件數(shù)的計算,這里不再贅述。

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

對于布爾莎模型可能產(chǎn)生的病態(tài)法方程(非嚴重病態(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),重復②至④,給定一個微小量ε,

上述①中的“全選主元高斯消去法”很多參考書可找到,這里將其成員函數(shù)命名為GetRootGuass()。病態(tài)方程求解程序?qū)崿F(xià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 算例

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

圖1 原坐標文件Fig.1 Original coordinate file

圖2 目標坐標文件Fig.2 Target coordinate file

程序主界面如圖3。

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

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

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

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

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

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

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

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

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

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

6 結(jié)語

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

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