吳厚清　蔡以雷

（浙江省測(cè)繪資料檔案館　浙江杭州　310000）

四參數(shù)仿射變換在大數(shù)據(jù)整合中的實(shí)踐與探討

吳厚清蔡以雷

（浙江省測(cè)繪資料檔案館浙江杭州310000）

摘要：根據(jù)布爾莎七參數(shù)模型轉(zhuǎn)換原理，靈活運(yùn)用四參數(shù)模型對(duì)未知坐標(biāo)系的空間數(shù)據(jù)進(jìn)行了轉(zhuǎn)換糾正，采用迭代方式保證其轉(zhuǎn)換成果的精度，能滿足測(cè)繪與地理信息部門在大數(shù)據(jù)時(shí)代對(duì)獲取的數(shù)據(jù)進(jìn)行整合利用的需求，并為未知坐標(biāo)系的數(shù)據(jù)整合提供新的思路。

關(guān)鍵詞：數(shù)據(jù)整合；布爾莎七參數(shù)；坐標(biāo)轉(zhuǎn)換

1　引言

在大數(shù)據(jù)時(shí)代，測(cè)繪與地理信息部門即是數(shù)據(jù)的生產(chǎn)者，也是數(shù)據(jù)的整合者。在工作實(shí)踐中，往往要探討對(duì)不同部門和行業(yè)獲取的數(shù)據(jù)進(jìn)行整合，數(shù)據(jù)整合的前提是對(duì)不同的數(shù)據(jù)進(jìn)行統(tǒng)一的坐標(biāo)定位，也即必須進(jìn)行坐標(biāo)轉(zhuǎn)換。

2　常用坐標(biāo)系統(tǒng)

來(lái)自不同部門和行業(yè)的數(shù)據(jù)，通常具有不同的坐標(biāo)系，基本的有大地坐標(biāo)系和平面直角坐標(biāo)系二大類，而大地坐標(biāo)系又分地心坐標(biāo)系和參心坐標(biāo)系，采用不同的橢球幾何參數(shù)。常見(jiàn)的有1984世界大地坐標(biāo)系 （即WGS-84坐標(biāo)系）、1954年北京坐標(biāo)系、1980西安坐標(biāo)系、2000國(guó)家大地坐標(biāo)系等。

目前世界各國(guó)采用最廣泛的高斯-克呂格投影和墨卡托投影（utm）均是正形投影（等角投影），即該投影在小區(qū)域范圍內(nèi)使平面圖形與橢球面上的圖形保持相似。為了限制長(zhǎng)度變形，根據(jù)國(guó)際測(cè)量協(xié)會(huì)規(guī)定，將全球按一定經(jīng)差分成若干帶。我國(guó)采用6度分帶或3度分帶。

高斯平面直角坐標(biāo)系一般以中央經(jīng)線投影為縱軸x，赤道投影為橫軸y，兩軸交點(diǎn)即為各帶的坐標(biāo)原點(diǎn)。為了避免橫坐標(biāo)出現(xiàn)負(fù)值，在投影中規(guī)定將坐標(biāo)縱軸西移500公里當(dāng)作起始軸。為了區(qū)別某一坐標(biāo)系統(tǒng)屬于哪一帶，通常在橫軸坐標(biāo)前加上帶號(hào)，如（4231898m，21655933m），其中21即為帶號(hào)。城建坐標(biāo)多采用三度帶的高斯-克呂格投影。同一坐標(biāo)系下的大地坐標(biāo)（即經(jīng)緯度坐標(biāo)b，l）與其對(duì)應(yīng)的高斯平面直角坐標(biāo)（x，y）有嚴(yán)格的轉(zhuǎn)換關(guān)系。

3　坐標(biāo)系統(tǒng)轉(zhuǎn)換

關(guān)于坐標(biāo)轉(zhuǎn)換，首先要搞清楚轉(zhuǎn)換的嚴(yán)密性問(wèn)題，即在同一個(gè)橢球里的坐標(biāo)轉(zhuǎn)換都是相對(duì)嚴(yán)密的，而在不同的橢球之間的轉(zhuǎn)換是相對(duì)近似的。例如，由1954北京坐標(biāo)系的大地坐標(biāo)轉(zhuǎn)換到1954北京坐標(biāo)系的高斯平面直角坐標(biāo)是在同一參考橢球體范疇內(nèi)的坐標(biāo)轉(zhuǎn)換，其轉(zhuǎn)換過(guò)程是嚴(yán)密的。由WGS-84的大地坐標(biāo)轉(zhuǎn)換到1954北京坐標(biāo)系，就屬于不同橢球體間的轉(zhuǎn)換。

不同橢球體間的坐標(biāo)轉(zhuǎn)換在局部地區(qū)的采用的常用辦法是相似變換法，即利用部分分布相對(duì)合理高等級(jí)公共點(diǎn)求出相應(yīng)的轉(zhuǎn)換參數(shù)。一般而言，比較嚴(yán)密的是用布爾莎七參數(shù)的相似變換法，即x平移，y平移，z平移，x旋轉(zhuǎn)，y旋轉(zhuǎn)，z旋轉(zhuǎn)，尺度變化k。要求得七參數(shù)就需要在一個(gè)地區(qū)需要3個(gè)以上的已知點(diǎn)，如果區(qū)域范圍不大，最遠(yuǎn)點(diǎn)間的距離不大于30km（經(jīng)驗(yàn)值），這可以用三參數(shù)，即x平移，y平移，z平移，而將x旋轉(zhuǎn)，y旋轉(zhuǎn)，z旋轉(zhuǎn)，尺度變化k視為0，所以三參數(shù)只是七參數(shù)的一種特例。很多研究表明，使用七參數(shù)模型換算時(shí)，大地高對(duì)于其平面坐標(biāo)換算的結(jié)果即平面坐標(biāo)（x，y）影響不大；如果不考慮高程的影響，對(duì)于不同橢球體下的高斯平面直角坐標(biāo)可采用四參數(shù)的相似變換法，即四參數(shù)（x平移，y平移，尺度變化m，旋轉(zhuǎn)角度α）。如果用戶要求的精度低于20米，在一定范圍（2＇*2＇）內(nèi)，就直接可以用二參數(shù)法（δb，δl）或（δx，δy）修正。但在實(shí)際操作中，這也取決于選取的公共點(diǎn)是否合理，并保證其足夠的精度。

4四　參數(shù)仿射變換實(shí)例

對(duì)未知坐標(biāo)系數(shù)據(jù)的轉(zhuǎn)換，在無(wú)需考慮高程的影響下，可采用四參數(shù)仿射變換，即四參數(shù)（x平移，y平移，尺度變化k，旋轉(zhuǎn)角度α）來(lái)轉(zhuǎn)換未知坐標(biāo)系統(tǒng)的數(shù)據(jù)，同時(shí)可適當(dāng)增加選取公共點(diǎn)，來(lái)保證其足夠的精度。在浙江省地理空間數(shù)據(jù)交換和共享平臺(tái)數(shù)據(jù)整合實(shí)踐中，獲取了浙江省統(tǒng)計(jì)局的第三次經(jīng)濟(jì)普查工作中產(chǎn)生的浙江省范圍內(nèi)的建筑物地址數(shù)據(jù)（不屬于2000坐標(biāo)系數(shù)據(jù)），這些地址數(shù)據(jù)都存在相應(yīng)的普查區(qū)工作邊界，而普查區(qū)工作邊界與浙江省縣及鄉(xiāng)鎮(zhèn)級(jí)以上的行政區(qū)域境界存在高度的相關(guān)性，利用邊界之間的相關(guān)性來(lái)獲取全省不同區(qū)域多個(gè)同名控制點(diǎn)，在實(shí)際操作中，為提高轉(zhuǎn)換精度，將浙江省區(qū)域范圍按第三次經(jīng)濟(jì)普查工作劃分的普查工作區(qū)為基本單元，運(yùn)用四參數(shù)仿射變換并采用迭代方式，分別求得每個(gè)基本單元的轉(zhuǎn)換參數(shù)，就可以得到我們需要的2000坐標(biāo)系成果。

4.1技術(shù)路線

由于矢量數(shù)據(jù)都是由點(diǎn)組成，只要利用點(diǎn)數(shù)據(jù)就可以根據(jù)四參數(shù)的仿射變換的計(jì)算公式進(jìn)行轉(zhuǎn)換。其中最關(guān)鍵的是找出同名控制點(diǎn)在兩個(gè)不同坐標(biāo)系中的不同坐標(biāo)。四參數(shù)的仿射變換的計(jì)算公式如下：x1=a*x0+b*y0+c，y1=m*x0+n* y0+d，其中（a,b,c,m,n,d）就是需要求的參數(shù)，其實(shí)中間是有內(nèi)在關(guān)系可以省略掉兩個(gè)參數(shù)的，最后需要求一個(gè)旋轉(zhuǎn)角度、兩個(gè)偏移參數(shù)，再加一個(gè)尺度變化。

4.2程序?qū)崿F(xiàn)

4.2.1獲取同名控制點(diǎn)

首先利用普查區(qū)工作邊界跟浙江省縣及鄉(xiāng)鎮(zhèn)級(jí)以上的境界的高度相關(guān)性的特點(diǎn)，各自取得界線相交點(diǎn)作為同名控制點(diǎn)在不同坐標(biāo)系中的坐標(biāo)數(shù)據(jù)組，控制點(diǎn)部分截圖示例如圖1

4.2.2獲取旋轉(zhuǎn)角度

根據(jù)對(duì)角兩個(gè)控制點(diǎn)數(shù)據(jù)獲取旋轉(zhuǎn)角度參數(shù)。

＇〈summary〉

＇獲取旋轉(zhuǎn)角度

＇〈/summary〉

＇〈paramname="fromCoordPoint1"〉源點(diǎn) 1〈/ param〉

＇〈paramname="toCoordPoint1"〉目標(biāo)點(diǎn) 1〈/ param〉

＇〈paramname="fromCoordPoint2"〉源點(diǎn) 2〈/ param〉

＇〈paramname="toCoordPoint2"〉目標(biāo)點(diǎn) 2〈/ param〉

＇〈returns〉返回旋轉(zhuǎn)角度〈/returns〉

PrivateFunctionGetRotation（fromPoint1As IPoint,toPoint1AsIPoint,fromPoint2AsIPoint, toPoint2AsIPoint）AsDouble

Dima,bAsDouble

a=（toPoint2.y-toPoint1.y）*（fromPoint2.xfromPoint1.x）-（toPoint2.x-toPoint1.x）*（fromPoint2.y-fromPoint1.y）

b=（toPoint2.x-toPoint1.x）*（fromPoint2.xfromPoint1.x）+（toPoint2.y-toPoint1.y）*（fromPoint2.y-fromPoint1.y）

If（Math.Abs（b）〉0）Then

GetRotation=Math.Tan（a/b）

Else

GetRotation=Math.Tan（0）

EndIf

EndFunction

4.2.3獲取縮放比例

＇〈summary〉

＇根據(jù)對(duì)角兩個(gè)控制點(diǎn)數(shù)據(jù)和旋轉(zhuǎn)角度獲取尺度變化參數(shù)。

＇〈/summary〉

＇〈paramname="fromCoordPoint1"〉源點(diǎn) 1〈/ param〉

＇〈paramname="toCoordPoint1"〉目標(biāo)點(diǎn) 1〈/ param〉

＇〈paramname="fromCoordPoint2"〉源點(diǎn) 2〈/ param〉

＇〈paramname="toCoordPoint2"〉目標(biāo)點(diǎn) 2〈/ param〉

＇〈paramname="rotation"〉旋轉(zhuǎn)角度〈/param〉

＇〈returns〉返回縮放比例〈/returns〉

PrivateFunctionGetScale（fromPoint1AsIPoint, toPoint1AsIPoint,fromPoint2AsIPoint,toPoint2As IPoint,rotationAsDouble）AsDouble

Dima,bAsDouble

a=toPoint2.x-toPoint1.x

b=（fromPoint2.x-fromPoint1.x）*Math.Cos（rotation）-（fromPoint2.y-fromPoint1.y）*Math.Sin（rotation）

If（Math.Abs（b）〉0）Then

GetScale=a/b

Else

GetScale=0

EndIf

EndFunction

4.2.4獲取X、Y偏移量

＇〈summary〉

＇根據(jù)一個(gè)控制點(diǎn)數(shù)據(jù)、旋轉(zhuǎn)角度和尺度變化獲取X、Y兩個(gè)方向的偏移量參數(shù)。

＇〈/summary〉

＇〈paramname="fromCoordPoint1"〉源點(diǎn) 1〈/ param〉

＇〈paramname="toCoordPoint1"〉目標(biāo)點(diǎn) 1〈/ param〉

＇〈paramname="rotation"〉旋轉(zhuǎn)角度〈/param〉

＇〈paramname="scale"〉尺度變化〈/param〉

＇〈returns〉返回X方向偏移量〈/returns〉

PrivateFunctionGetXTranslation（fromPoint1As IPoint,toPoint1AsIPoint,rotationAsDouble,scale1 AsDouble）AsDouble

GetXTranslation=（toPoint1.x-scale1*（fromPoint1.x*Math.Cos（rotation）-fromPoint1.y* Math.Sin（rotation）））

EndFunction

＇〈summary〉

＇得到Y(jié)方向偏移量

＇〈/summary〉

＇〈paramname="fromCoordPoint1"〉源點(diǎn) 1〈/ param〉

＇〈paramname="toCoordPoint1"〉目標(biāo)點(diǎn) 1〈/ param〉

＇〈paramname="rotation"〉旋轉(zhuǎn)角度〈/param〉

＇〈paramname="scale"〉尺度變化〈/param〉

＇〈returns〉返回Y方向偏移量〈/returns〉

PrivateFunctionGetYTranslation（fromPoint1As IPoint,toPoint1AsIPoint,rotationAsDouble,scale1 AsDouble）AsDouble

GetYTranslation=（toPoint1.y-scale1*（fromPoint1.x*Math.Sin（rotation）+fromPoint1.y* Math.Cos（rotation）））

EndFunction

4.2.5調(diào)用四參數(shù)變換模型

＇說(shuō)明：采用相似變換模型（四參數(shù)變換模型）

＇X=ax-by+c

＇Y=bx+ay+d

PrivateFunctionUnaffineTransformation（pPoint AsIPoint）AsIPoint

DimaAsDouble,bAsDouble

DimtPointAsIPoint

GetFourParampPoint＇得到四參數(shù)

a=m_Scale*Math.Cos（m_RotationAngle）b=m_Scale*Math.Sin（m_RotationAngle）SettPoint=NewPoint

tPoint.x=a*pPoint.x-b*pPoint.y+m_DX tPoint.y=b*pPoint.x+a*pPoint.y+m_DY SetUnaffineTransformation=tPoint EndFunction

4.3坐標(biāo)轉(zhuǎn)換的成果

通過(guò)對(duì)浙江省統(tǒng)計(jì)局全省普查邊界數(shù)據(jù)的轉(zhuǎn)換，進(jìn)而獲取普查工作區(qū)范圍內(nèi)的2000坐標(biāo)系的建筑物地址數(shù)據(jù)。經(jīng)檢測(cè)，經(jīng)坐標(biāo)轉(zhuǎn)換后的數(shù)據(jù)平均誤差在3米以內(nèi)，滿足同比例尺地形圖的精度要求。如下：圖2為轉(zhuǎn)換前局部數(shù)據(jù)截圖，其平均誤差在500多米，圖3為轉(zhuǎn)換后局部數(shù)據(jù)截圖，平均誤差在3米以內(nèi)。

5　結(jié)論

經(jīng)過(guò)浙江省統(tǒng)計(jì)局的第三次經(jīng)濟(jì)普查工作產(chǎn)生的數(shù)據(jù)轉(zhuǎn)換項(xiàng)目，我們可以更深的理解不同坐標(biāo)系之間的轉(zhuǎn)換其實(shí)是兩種不同的橢球參數(shù)之間的轉(zhuǎn)換。其中七參數(shù)的布爾莎模型，是比較嚴(yán)密的轉(zhuǎn)換模型，不管是四參數(shù)、三參數(shù)、二參數(shù)都是其特定情況的簡(jiǎn)化轉(zhuǎn)換公式。在轉(zhuǎn)換未知坐標(biāo)系數(shù)據(jù)過(guò)程中，靈活套用布爾莎模型的簡(jiǎn)化轉(zhuǎn)換公式，是一種比較經(jīng)濟(jì)實(shí)用的選擇。

參考文獻(xiàn)：

［1］陳俊勇.《測(cè)繪通報(bào)》2003年02期《世界大地坐標(biāo)系1984的最新精化》

［2］程鵬飛，文漢江，成英燕，王華.《測(cè)繪學(xué)報(bào)》2009年03期《2000國(guó)家大地坐標(biāo)系橢球參數(shù)與GRS80和WGS84的比較》

［3］徐登云，郝麗娟.《西部資源》2012年02期 《2000國(guó)家大地坐標(biāo)系與GRS80及WGS84的比較》