【文章內(nèi)容簡(jiǎn)介】 是縮放參數(shù)），Dx、Dy、Dz為平移參數(shù)，ωx、ωy、ωz為旋轉(zhuǎn)參數(shù)（弧度）。有些文獻(xiàn)上縮放參數(shù)為μ，其單位是ppm，它與K的關(guān)系如下：布爾莎模型的特點(diǎn)：計(jì)算簡(jiǎn)單；它是對(duì)相似變換的近似描述，并不是嚴(yán)格的相似變換；不能找到一套布爾莎七參數(shù)用于描述其逆變換。就是說(shuō)，如果已知坐標(biāo)系1至坐標(biāo)系2的布爾莎七參數(shù)，是無(wú)法得到坐標(biāo)系2至坐標(biāo)系1的布爾莎七參數(shù)的。要求逆變換必須使用轉(zhuǎn)換矩陣的逆陣： 并不是所有GPS軟件都使用布爾莎模型，這里將介紹DSNP+模型。它是DSNP公司3SPack軟件采用的七參數(shù)模型。介紹它的意義在于：它是一種嚴(yán)格的三維相似變換。姿態(tài)儀定位、遙感圖像通用構(gòu)象方程都可以使用該模型。規(guī)定坐標(biāo)系1至坐標(biāo)系2的DSNP+七參數(shù)意義如下： 根據(jù)(234)求解七參數(shù)的過(guò)程中，往往會(huì)發(fā)現(xiàn)法方程是病態(tài)的。原因在于(234)式在旋轉(zhuǎn)的時(shí)候是以原點(diǎn)為旋轉(zhuǎn)中心的，原點(diǎn)離參與計(jì)算的點(diǎn)都很遠(yuǎn)，所以旋轉(zhuǎn)角的微小誤差將會(huì)引起最終旋轉(zhuǎn)結(jié)果較大的偏差。那么這種距離對(duì)最終的法方程是否病態(tài)有多大影響呢？假定坐標(biāo)系1有三個(gè)坐標(biāo)已知的點(diǎn)，分別為 上面的定量結(jié)果只是一個(gè)特例，并不適用于一般情況，但是有理由相信：N的條件數(shù)與r2成正比，隨著旋轉(zhuǎn)距離r的增大，N的條件數(shù)快速增大，法方程迅速轉(zhuǎn)為病態(tài)。所以必須盡量減小旋轉(zhuǎn)距離r。如果不解算七參數(shù)中的平移參數(shù)Dx、Dy、Dz，則 這說(shuō)明：如果不解算平移參數(shù)，則法方程不會(huì)轉(zhuǎn)入病態(tài)。 當(dāng)解算參數(shù)里含有平移參數(shù)時(shí)，就必須解決法方程的病態(tài)問(wèn)題。下面是一種解決方法。 求解七參數(shù)前，先對(duì)坐標(biāo)系2中的坐標(biāo)做平移，形成坐標(biāo)139。、239。公式表達(dá)如下：大部分情況下，按(250)式進(jìn)行間接平差就能很好的解決法方程病態(tài)問(wèn)題。不過(guò)需要注意平差后的平移參數(shù)Dx39。、Dy39。、Dz39。并不是最終想要的可以通過(guò)下式獲得最終需要的平移參數(shù)：在傳統(tǒng)的七參數(shù)轉(zhuǎn)換中，往往僅考慮單向轉(zhuǎn)換。如：一般只需要考慮WGS84至北京54的坐標(biāo)轉(zhuǎn)換。求解七參數(shù)時(shí)，往往將WGS84中的坐標(biāo)當(dāng)成已知量，北京54中的坐標(biāo)當(dāng)成未知的隨機(jī)變量，通過(guò)最小二乘法求解七參數(shù)。這種情況下，七參數(shù)模型顯然偏重于北京54坐標(biāo)，如果用這套參數(shù)反向計(jì)算，即將北京54坐標(biāo)轉(zhuǎn)換為WGS84坐標(biāo)，則轉(zhuǎn)換出來(lái)的結(jié)果就不會(huì)很理想。 舉例說(shuō)明：人的身高H和體重W是有一定線性關(guān)系的，即W=aH+b。為了求解這個(gè)模型中的參數(shù)a、b，量取了若干人的身高和體重。如果僅需要通過(guò)H計(jì)算W，那么最小二乘法的約束條件就是[10]：如果僅需要通過(guò)W計(jì)算H，那么最小二乘法的約束條件就是：如果需要H和W相互轉(zhuǎn)換，則約束條件(253)和(254)都是不公平的，必須提供新的約束條件。本文提出了一種約束方法：這個(gè)方法是對(duì)(253)和(254)的簡(jiǎn)單綜合。對(duì)于該方法而言在組成誤差方程的時(shí)候，一組數(shù)據(jù)（Hi，Wi）將生成兩個(gè)誤差方程，而不是通常的一個(gè)，其余的計(jì)算與普通的間接平差沒(méi)有區(qū)別，是切實(shí)可行的。以求解坐標(biāo)系1至坐標(biāo)系2的布爾莎七參數(shù)為例，將會(huì)有三種最小二乘約束： （1）僅需要通過(guò)坐標(biāo)1計(jì)算坐標(biāo)2假定某個(gè)已知點(diǎn)在坐標(biāo)系1和坐標(biāo)系2的坐標(biāo)分別為 ，通過(guò)下式可以計(jì)算出該點(diǎn)在坐標(biāo)系2的坐標(biāo)： 可以稱(257)這種最小二乘約束方式為“盡量靠近坐標(biāo)系2”（2）僅需要通過(guò)坐標(biāo)2計(jì)算坐標(biāo)1假定某個(gè)已知點(diǎn)在坐標(biāo)系1和坐標(biāo)系2的坐標(biāo)分別為和，通過(guò)下式可以計(jì)算出該點(diǎn)在坐標(biāo)系1的坐標(biāo)：可以稱(259)這種最小二乘約束方式為“盡量靠近坐標(biāo)系1”。（3）需要坐標(biāo)1和坐標(biāo)2相互轉(zhuǎn)換假定某個(gè)已知點(diǎn)在坐標(biāo)系1和坐標(biāo)系2的坐標(biāo)分別為 和，通過(guò)公式(256)和(258)可以計(jì)算出該點(diǎn)在坐標(biāo)系2的計(jì)算坐標(biāo)和。按下面的約束條件進(jìn)行七參數(shù)求解：可以稱(260)這種最小二乘約束方式為“盡量靠近坐標(biāo)系1和2”。已知A、B、C三點(diǎn)在坐標(biāo)系2下的空間直角坐標(biāo)（），求坐標(biāo)系1至坐標(biāo)系2的布爾莎七參數(shù)，分別以(257)、(259)、(260)為約束進(jìn)行求解。首先，求出A、B、C三點(diǎn)的平均坐標(biāo)如下：然后對(duì)A、B、C做平移，平移結(jié)果如下： 最小二乘約束為(259)模型(258)是非線性的，組成誤差方程的時(shí)候需要通過(guò)求偏導(dǎo)來(lái)完成線性化，這是一個(gè)比較復(fù)雜的過(guò)程。好在這里只需要數(shù)值解，所以可以通過(guò)數(shù)值微分的方法來(lái)組成法方程。待求參數(shù)取初始值為則誤差方程為具體形式如下：法方程為，式中。注意這是一個(gè)非線性的模型，不能求解一次法方程就得到七參數(shù)，必須迭代計(jì)算[10][11]。最終的計(jì)算結(jié)果為： 最小二乘約束為(257) 待求參數(shù)取初始值為 將上式中的平移參數(shù)Dx39。、Dy39。、Dz39。代入(252)可算得所需的平移參數(shù)。最終的計(jì)算結(jié)果如下： 最小二乘約束為(260) 誤差方程是前面兩者誤差方程的組合，即誤差方程V =Bδ x L中殘差記將上面計(jì)算出來(lái)的三套參數(shù)(264)、(269)、(271)分別代入模型(256)、(258)可以計(jì)算出δδ2。 上表中，每行的最小值都由下劃線標(biāo)出。第1套參數(shù)計(jì)算出的δ1最小，第2套參數(shù)計(jì)算出來(lái)的δ2最小，第3套參數(shù)計(jì)算出來(lái)的δ1+δ2最小，這些與約束條件(259)、(257)、(260)是相符合的。GPS測(cè)得的坐標(biāo)是WGS84大地坐標(biāo)L、B、H，直接使用大地坐標(biāo)是很不方便的。如：給定A、B點(diǎn)的大地坐標(biāo)，如何計(jì)算A至B的方位和距離？解決的辦法就是將A、B點(diǎn)按一定的規(guī)則轉(zhuǎn)換到一個(gè)平面內(nèi)，這個(gè)轉(zhuǎn)換規(guī)則就是地圖投影。地圖投影有很多方式[12]，如：等角投影、等距投影、等積投影。每種方式都有若干的計(jì)算方法，在此就不一一介紹了。本文著重介紹在工程中應(yīng)用十分廣泛的高斯——克呂格投影?！藚胃裢队?，假想有一個(gè)橢圓柱面橫套在地球橢球體外面，并與某一條子午線（此子午線稱為中央子午線或軸子午線）相切，橢圓柱的中心軸通過(guò)橢球體中心，然后用一定投影方法，將中央子午線兩側(cè)各一定經(jīng)差范圍內(nèi)的地區(qū)投影到橢圓柱面上，再將此柱面展開(kāi)即成為投影面。此投影為高斯投影。高斯投影是等角投影的一種?！藚胃裢队? 在投影面上，中央子午線和赤道的投影都是直線，并且以中央子午線和赤道的交點(diǎn)O作為坐標(biāo)原點(diǎn)，以中央子午線的投影為N軸（北坐標(biāo)），以赤道的投影為E軸（東坐標(biāo)）。高斯投影的特點(diǎn)：中央子午線無(wú)變形；小范圍內(nèi)無(wú)角度變形，圖形保持相似；離中央子午線越遠(yuǎn)，變形越大。高斯投影的正算公式（由大地經(jīng)緯度計(jì)算投影坐標(biāo)）上式中l(wèi)為經(jīng)度差——投影點(diǎn)大地經(jīng)度減去中心經(jīng)度的差值。a0為赤道至投影點(diǎn)大地緯度B的子午線弧長(zhǎng)。注意：a1~a6中的N是卯酉圈曲率半徑，不要和北坐標(biāo)N混淆了。高斯投影的反算公式（由投影坐標(biāo)計(jì)算大地經(jīng)緯度）Bf為底點(diǎn)緯度：赤道至Bf的子午線弧長(zhǎng)等于投影點(diǎn)的北坐標(biāo)N。Nf是底點(diǎn)緯度對(duì)應(yīng)的卯酉圈曲率半徑。子午線弧長(zhǎng)的計(jì)算是高斯投影正反算的基礎(chǔ)，所以很重要。文獻(xiàn)[13][14]對(duì)計(jì)算公式做了適當(dāng)簡(jiǎn)化，但沒(méi)有從根本上解決編程復(fù)雜、精度不易擴(kuò)充的缺陷，本文就此進(jìn)行了深入研究。，點(diǎn)P在赤道上，Q是子午線上的一點(diǎn)，其大地緯度為B。弧PQ的長(zhǎng)度就是待求的子午線弧長(zhǎng)，顯然它是大地緯度B的函數(shù)，記作X(B)。附錄A中的CEllipse::LatToNorth函數(shù)就是使用公式(290)、(291)、(292)編制的計(jì)算子午線弧長(zhǎng)的代碼。這段代碼非常簡(jiǎn)單，而且計(jì)算精度易于擴(kuò)充。高斯投影反算時(shí)，需要求底點(diǎn)緯度Bf，使其滿足X (B f)=N（N是投影北坐標(biāo)），亦即求解方程f (B )=X (B )N=0的根。文獻(xiàn)[15]就此提出了牛頓迭代法，實(shí)際計(jì)算表明該方法收斂很快，是個(gè)非常好的方法。由此可知函數(shù)X (B)在[0,π]之間的圖形（）。求解底點(diǎn)緯度的牛頓迭代公式為 （1）若N0則記bNeg為T(mén)RUE（表示N和底點(diǎn)緯度均為負(fù)數(shù)），否則bNeg為FALSE； （2）若bNeg為T(mén)RUE則將N的相反數(shù)賦予N，保證N≥0； （5）按(295)迭代計(jì)算，得到底點(diǎn)緯度Bf； （6）若bNeg為T(mén)RUE則底點(diǎn)緯度為負(fù)，需要在第（5）步的Bf前加負(fù)號(hào)。 附錄A中的CEllipse::NorthToLat函數(shù)就是使用上述步驟進(jìn)行底點(diǎn)緯度計(jì)算的。 可見(jiàn)迭代公式(295)對(duì)迭代初始值的取法略顯繁瑣。下面的迭代公式可以取任意的迭代初始值：不過(guò)(296)比(295)收斂慢得多，計(jì)算量也因此增大。、插值GPS測(cè)定的高程是大地高，而實(shí)際工程中往往需要正常高（海拔高程）。兩者之間的換算關(guān)系很簡(jiǎn)單：正常高=大地高－高程異常。但是高程異常的變化異常復(fù)雜，很難用一套公式精確的描述。因此，在知道有限個(gè)位置處的高程異常后，進(jìn)行插值計(jì)算，獲取待定點(diǎn)的高程異常就顯得非常的重要。上式中L0 、 B0是某一特定點(diǎn)的大地經(jīng)緯度，為已知值。如果事先知道若干點(diǎn)的大地經(jīng)緯度和高程異常，則使用多項(xiàng)式模型前要計(jì)算矩陣A?？梢匀0 、 B0為這些點(diǎn)的平均經(jīng)緯度或其它的固定值，然后使用間接平差的方法來(lái)計(jì)算系數(shù)aij。平差后一般都會(huì)有殘差，所以多項(xiàng)式模型屬于高程異常擬合。距離加權(quán)平均是一種加權(quán)平均法[17]，權(quán)一般由距離來(lái)定，所以稱之為距離加權(quán)平均。假定已知n個(gè)點(diǎn)的大地經(jīng)緯度和高程異常ζ，則任意一點(diǎn)的高程異常計(jì)算公式如下：，高程異常已知的點(diǎn)都在橫豎線的交叉點(diǎn)上規(guī)則排列，這就是格網(wǎng)模型，交叉點(diǎn)被稱之為節(jié)點(diǎn)。給定任一點(diǎn)，都可以找到相鄰的幾個(gè)節(jié)點(diǎn)，然后按一定的公式內(nèi)插計(jì)算出該點(diǎn)的高程異常。，已知A、 B、 C、 D點(diǎn)的高程異常ζA、 ζB 、ζC 、ζD，需要求P點(diǎn)的高程異常ζ。雙線性內(nèi)插的做法就是：AB段線性內(nèi)插出a點(diǎn)的高程異常，CD段線性內(nèi)插出c點(diǎn)的高程異常，然后ac段線性內(nèi)插出P點(diǎn)的高程異常。同理，也可以使用bd段線性內(nèi)插出P點(diǎn)的高程異常，兩種方法計(jì)算出來(lái)的結(jié)果相同。最終的計(jì)算結(jié)果為[18]上式中ai是與ζA、 ζB 、ζC 、ζD有關(guān)的系數(shù)。可見(jiàn)雙線性內(nèi)插是一種二階的多項(xiàng)式模型。根據(jù)已知數(shù)據(jù)可以組成線性方程組：，O點(diǎn)的高程異常已知，網(wǎng)格橫向間距為m，豎向間距為n。根據(jù)導(dǎo)數(shù)的數(shù)值算法，有[19]：將方程組(2101)中B1里的偏導(dǎo)數(shù)全部使用數(shù)值方法計(jì)算，則組成新的線性方程組： 對(duì)于高程異常的擬合、插值，本文介紹了兩大類模型：距離加權(quán)平均模型和多項(xiàng)式模型。距離加權(quán)平均模型的計(jì)算非常簡(jiǎn)單，但是權(quán)的取值非常重要。為了更加真實(shí)的描述實(shí)際情況，可能需要多次嘗試才能選出合適的定權(quán)方法。多項(xiàng)式模型的計(jì)算也比較簡(jiǎn)單，但是它的階數(shù)不能太高，否則容易產(chǎn)生龍格（Runge）現(xiàn)象[8]；階數(shù)也不能太低，否則無(wú)法反映高程異常的復(fù)雜變化。對(duì)于地形起伏較小的平坦地區(qū)或沿海地區(qū)，高程異常變化不大，使用雙線性多項(xiàng)式甚至一次多項(xiàng)式都可以很好的進(jìn)行擬合[16]。對(duì)于高程異常變化比較大的地區(qū)，可以建立格網(wǎng)，進(jìn)行分塊擬合。一般每小塊之間使用雙三次多項(xiàng)式擬合就足夠了。建立格網(wǎng)降低了多項(xiàng)式模型的階數(shù)，消除了龍格現(xiàn)象，同時(shí)能夠描述復(fù)雜的高程異常變化，擬合效果比較好，但是格網(wǎng)的建立工作比較繁重，有時(shí)限制了它在工程中的應(yīng)用。 GPS坐標(biāo)轉(zhuǎn)換GPS獲得的坐標(biāo)一般都是WGS84大地坐標(biāo)，而工程中需要的是施工坐標(biāo)和海拔高程。GPS接收機(jī)一般都有坐標(biāo)轉(zhuǎn)換功能，可以實(shí)現(xiàn)這個(gè)轉(zhuǎn)換。但是對(duì)于不同的接收機(jī)這個(gè)配置是不同，需要時(shí)間去學(xué)習(xí)、驗(yàn)證，一旦接收機(jī)沒(méi)有所需的轉(zhuǎn)換功能就比較難于處理了。所以，編制水下地形測(cè)量軟件時(shí)，最好有坐標(biāo)轉(zhuǎn)換功能，以便靈活的接各種型號(hào)的GPS接收機(jī)。一般將WGS84大地坐標(biāo)轉(zhuǎn)換為施工坐標(biāo)的過(guò)程如下：將WGS84大地坐標(biāo)轉(zhuǎn)換為空間直角坐標(biāo)；使用七參數(shù)將第1步的空間直角坐標(biāo)轉(zhuǎn)換到用戶橢球（如：克拉索夫斯基橢球）；將第2步獲得的用戶空間直角坐標(biāo)轉(zhuǎn)換為用戶橢球大地坐標(biāo)；使用第3步的用戶大地經(jīng)緯度在用戶橢球上進(jìn)行投影，獲得投影坐標(biāo)；如果投影坐標(biāo)系與施工坐標(biāo)系不重合，還需要對(duì)投影坐標(biāo)做二維相似變換，變?yōu)槭┕て矫孀鴺?biāo)；計(jì)算高程異常；大地高減去第6步的高程異常得到施工需要的正常高（海拔高程）。這里的大地高可以是第1步中的WGS84大地高，也可以是第3步求得的用戶橢球大地高，具體選哪個(gè)要根據(jù)模型來(lái)定。實(shí)際工程中，GPS坐標(biāo)轉(zhuǎn)換的步驟和方法多種多樣，如：參考文獻(xiàn)[20]提到的EE2橢球。通過(guò)軟件定義一個(gè)標(biāo)準(zhǔn)的轉(zhuǎn)換是不能完全解決所有實(shí)際問(wèn)題的，因此，在編制水下地形測(cè)量軟件時(shí)必須提供坐標(biāo)轉(zhuǎn)換接口，以便針對(duì)不同的工程需求，進(jìn)行不同的坐標(biāo)轉(zhuǎn)換。第3章GPS坐標(biāo)傳遞GPS測(cè)得的坐標(biāo)僅僅是天線相位中心的位置，要將這個(gè)坐標(biāo)傳遞給船體的任一點(diǎn)（如：測(cè)深儀換能器）就必須建立船體坐標(biāo)系。GPS坐標(biāo)的傳遞，就是船體坐標(biāo)系與施工坐標(biāo)系的轉(zhuǎn)換關(guān)系。U軸從船艉指向船艏，V軸由船左指向船右，W軸是鉛垂線，指向天空。W軸的零點(diǎn)：水平面與靜止船體的相交平面就是W=0。如果測(cè)量船上僅有一臺(tái)GPS，那么測(cè)深的時(shí)候必須將測(cè)深儀換能器與GPS天線安裝在同一鉛垂線上，保證兩者的平面坐標(biāo)相同。至于高程只需加上一個(gè)常數(shù)即可。測(cè)量船上裝有兩臺(tái)GPS，天線2的相位中心在船體坐標(biāo)系下的坐標(biāo)固定為u1，v1，w1和u2，v2，w2。天線2的實(shí)時(shí)施工坐標(biāo)為x1，y1，h1和x2，y2，h2。U軸在施工坐標(biāo)系的方位角α等于天線1至2在施工坐標(biāo)系下的