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