王炎，田雨

（中交上海航道勘察設(shè)計研究院有限公司，上海 200120）

河口海岸地區(qū)經(jīng)濟(jì)發(fā)達(dá)，人類活動頻繁，大型工程密集，地形演變事關(guān)堤岸、航運(yùn)及港口安全，是地理學(xué)和工程學(xué)研究的主要方向[1-4]。由于該地區(qū)環(huán)境復(fù)雜，常規(guī)測繪方法效率低、成本高，通常伴隨大量多余觀測，利用當(dāng)前的測量設(shè)備需耗費(fèi)大量時間才能完全覆蓋河口地區(qū)[5-7]。ASTER GDEM（The Advanced Spaceborne Thermal Emission and Reflection Radiometer Global Digital Elevation Mode）數(shù)據(jù)的發(fā)布使得公眾可以免費(fèi)獲取全球范圍內(nèi)分辨率30 m的地形數(shù)據(jù)[8-9]，已廣泛應(yīng)用于工程開發(fā)、地質(zhì)與災(zāi)害評估等領(lǐng)域[10]，但受采集方式的限制，ASTER GDEM數(shù)據(jù)缺乏對水下地形的描述。若能充分利用該地區(qū)的海圖數(shù)據(jù)與ASTER GDEM數(shù)據(jù)進(jìn)行融合，將是快速獲取水陸一體連續(xù)地形的有效途徑。

ASTER GDEM數(shù)據(jù)與海圖數(shù)據(jù)的融合主要存在兩個問題：一為空間參考系不統(tǒng)一[11]；二為原始數(shù)據(jù)精度不一致[12]。這會導(dǎo)致融合邊界處的高程值與周圍高程值存在明顯差異。一些學(xué)者對如何處理該高程差異開展了相關(guān)研究。如許家琨等[13]總結(jié)了我國海岸帶地形測量涉及到的垂直基準(zhǔn)面并探討了各基面間的轉(zhuǎn)換關(guān)系；Kuuskivi等[14]對融合邊界處不連續(xù)數(shù)據(jù)進(jìn)行平滑處理，但會改變已知高程，給融合結(jié)果帶來不必要的精度損失；Alca^ntara等[15]提取出SRTM與歷史地形圖的等高線信息并合并，通過內(nèi)插重新生成高程模型，但融合結(jié)果空間分辨率較低；凌峰等[12]對融合數(shù)據(jù)的高程差異空間相關(guān)性進(jìn)行分析，有效減小了兩種數(shù)據(jù)間的相對誤差，但缺乏絕對精度分析；受制于多源數(shù)據(jù)獲取機(jī)理與數(shù)據(jù)特性的差異，地形數(shù)據(jù)間的無縫融合仍是國內(nèi)外研究的難點(diǎn)[10]，特別是ASTER GDEM數(shù)據(jù)精度遠(yuǎn)低于海圖數(shù)據(jù)精度，有關(guān)其融合過程中的誤差校正研究尚未見報道。

本文以長江口南支測區(qū)為研究地區(qū)，分析了ASTER GDEM數(shù)據(jù)與海圖數(shù)據(jù)間的差異并對該地區(qū)ASTER GDEM數(shù)據(jù)與海圖數(shù)據(jù)進(jìn)行融合研究。針對ASTER GDEM原始數(shù)據(jù)精度較低的缺陷，利用曲面擬合法對其進(jìn)行誤差校正，不僅減少了融合數(shù)據(jù)間的高程差異，而且也增強(qiáng)了數(shù)據(jù)后期使用的可靠性。本方法可以實(shí)現(xiàn)區(qū)域尺度的水陸一體連續(xù)地形模型的快速構(gòu)建，為河口地區(qū)的野外測量任務(wù)規(guī)劃提供決策依據(jù)，從而減少多余觀測，提高測量效率，降低工作成本。

1 研究方法

1.1 海圖與ASTER GDEM數(shù)據(jù)差異分析

海圖是地圖的一種，是按照水上航行需求繪制的地圖，符合GB 12319—1998《中國海圖圖式》，包含了岸形、島嶼、礁石、水深、航標(biāo)和無線電導(dǎo)航臺等內(nèi)容[16]。海圖數(shù)據(jù)重點(diǎn)表示水深信息和影響艦船航行的海部要素，岸線以上的地形要素主要根據(jù)地形圖編繪，僅有少量地面控制點(diǎn)含有高程信息。

ASTER GDEM數(shù)據(jù)是美國國家航天局（NASA）與日本經(jīng)濟(jì)產(chǎn)業(yè)省（METI）根據(jù)Terra衛(wèi)星所獲取的150萬景近紅外影像，基于同軌立體攝影測量原理制作完成的，屬于地球觀測系統(tǒng)（EOS）計劃的一部分[17]，覆蓋范圍為 83°N—83°S，是迄今為止可為公眾免費(fèi)下載的最完整、最可靠的全球數(shù)字高程數(shù)據(jù)[8]。圖1展示了長江口地區(qū)ASTER GDEM數(shù)據(jù)圖像，從圖上可以看出，因?yàn)楣鈱W(xué)傳感器的采集特性，水陸交界處存在大量異常數(shù)據(jù)，因此在使用ASTER GDEM數(shù)據(jù)前，需結(jié)合當(dāng)?shù)氐匦翁卣鲗Ξ惓?shù)據(jù)進(jìn)行探測并剔除。

圖1 長江口地區(qū)ASTER GDEM數(shù)據(jù)圖像Fig.1 ASTER GDEM imagery of Yangtze estuary

表1列出了ASTER GDEM與海圖數(shù)據(jù)的基本信息對照，兩者在空間參考系方面存在較大差異?？臻g坐標(biāo)系方面，海圖數(shù)據(jù)采用的CGCS2000坐標(biāo)系與ASTER GDEM采用的WGS-84坐標(biāo)系的基本定義是一致的，參考橢球體非常接近，僅有扁率這一常數(shù)有細(xì)微差別，但在當(dāng)前測量精度水平下這種微小差值可以忽略；在垂直基準(zhǔn)方面，海圖數(shù)據(jù)采用理論最低潮面作為高程基準(zhǔn)，與ASTER GDEM數(shù)據(jù)所采用的EGM96高程基準(zhǔn)存在系統(tǒng)誤差，不可避免的存在陸地高程與水深數(shù)據(jù)不能有效銜接的缺陷，故需要通過水準(zhǔn)測量，建立兩種數(shù)據(jù)源垂直基準(zhǔn)間的轉(zhuǎn)換數(shù)學(xué)模型[18]。

表1 海圖與ASTER GDEM數(shù)據(jù)基本信息表Table 1 Basic information of bathymetric map and ASTER GDEM data

1.2 基于曲面擬合算法的ASTER GDEM誤差校正

ASTER GDEM數(shù)據(jù)在我國范圍內(nèi)都是通過無控制點(diǎn)測量所得，垂直精度約為±20 m，即使在精度最高的平原地區(qū)也達(dá)到±5 m[19-20]，因此無法直接用于高精度水陸一體連續(xù)地形模型的構(gòu)建，而必須對高程誤差進(jìn)行校正后再與海圖數(shù)據(jù)進(jìn)行融合。

傳統(tǒng)的誤差校正方法首先統(tǒng)計兩種數(shù)據(jù)重疊區(qū)域的高程差值，建立回歸方程，然后在誤差均方根最小的情況下對融合數(shù)據(jù)進(jìn)行變換[21]。該方法對參與融合的高程數(shù)據(jù)均做相同處理，實(shí)際上，由傳感器構(gòu)造特點(diǎn)，環(huán)境變化等因素產(chǎn)生的誤差組分通常具有空間相關(guān)性[22]，傳統(tǒng)方法沒有考慮到這一點(diǎn)，給融合結(jié)果帶來不必要的精度損失。

為了提高融合結(jié)果的數(shù)據(jù)精度，一種更好的辦法就是利用高精度地面控制點(diǎn)對ASTER GDEM數(shù)據(jù)進(jìn)行誤差校正。兩種高程數(shù)據(jù)間的關(guān)系可以表達(dá)為：

式中：H（x，y）表示坐標(biāo)（x，y）處的地面控制點(diǎn)高程值；A（x，y）表示坐標(biāo)（x，y）處的 ASTER GDEM 數(shù)據(jù)高程值，E（x，y）表示坐標(biāo)（x，y）處的 ASTER GDEM高程誤差?？紤]到部分誤差組分具有空間相關(guān)性，本文采用二次曲面擬合方法建立高程誤差E（x，y）與空間位置（x，y）之間的相關(guān)關(guān)系，進(jìn)而根據(jù)該相關(guān)關(guān)系擬合出整個區(qū)域的誤差曲面，實(shí)現(xiàn)區(qū)域內(nèi)ASTER GDEM高程數(shù)據(jù)的誤差校正。

二次曲面擬合法公式如下：

為了進(jìn)一步提高求解精度，參與擬合計算的點(diǎn)數(shù)應(yīng)多于未知參數(shù)個數(shù)，通過參數(shù)平差原理進(jìn)行求解[23]，誤差方程如下：

式中：a0，a1，…，a5即為擬合參數(shù)，首先根據(jù)控制點(diǎn)位坐標(biāo)（x，y）和對應(yīng)高程誤差 E（x，y）對其進(jìn)行求解，再反過來將擬合參數(shù)和其余ASTER GDEM像元坐標(biāo)代入式（3）即可求得相應(yīng)位置的高程誤差，從而生成誤差曲面，對區(qū)域內(nèi)ASTER GDEM高程數(shù)據(jù)進(jìn)行校正。

1.3 方法實(shí)現(xiàn)

基于上述方法構(gòu)建水陸一體連續(xù)地形模型的技術(shù)流程如圖2所示。第一步要對多源數(shù)據(jù)進(jìn)行空間基準(zhǔn)的統(tǒng)一，平面基準(zhǔn)統(tǒng)一至CGCS2000坐標(biāo)系，垂直基準(zhǔn)統(tǒng)一至1985國家高程基準(zhǔn)（似大地水準(zhǔn)面）。郭海榮[24]的研究表明，1985國家高程基準(zhǔn)點(diǎn)與EGM96高程基準(zhǔn)之間有35.7 cm的系統(tǒng)垂直偏差，且自東向西明顯增大，海圖數(shù)據(jù)提供有對應(yīng)年份的平均海面高度，以此建立兩種垂直基準(zhǔn)間的數(shù)學(xué)轉(zhuǎn)換模型。

圖2 基于海圖和ASTER GDEM數(shù)據(jù)融合構(gòu)建水陸一體連續(xù)地形模型技術(shù)流程圖Fig.2 Flowchartofthemethodthatintegratingbathymetric maps and ASTER GDEM data to build continuous terrain of water-land integration

隨后對ASTER GDEM數(shù)據(jù)進(jìn)行水陸分離。通常而言，ASTER GDEM中水面數(shù)據(jù)為一定值，可通過設(shè)置閾值進(jìn)行提取。但鏡面反射等因素導(dǎo)致的異常數(shù)據(jù)會對分離結(jié)果造成干擾，因此還要借助準(zhǔn)確的岸線信息。航行圖中定義的岸線位置為平均大潮高潮位時水陸分界的痕跡線，在堤防工程完善的地區(qū)也可借助其他遙感影像獲得岸線的位置和形態(tài)[25]。結(jié)合實(shí)地情況，融合航行圖與遙感影像提供的岸線位置，以此作為分界線，對ASTER GDEM進(jìn)行水陸分離。水面數(shù)據(jù)作為無效數(shù)據(jù)進(jìn)行剔除。陸地數(shù)據(jù)采用二次曲面擬合法進(jìn)行誤差校正并進(jìn)行精度評估。

海圖提供了所在區(qū)域的水深信息，但多為離散點(diǎn)和等深線，需要進(jìn)一步內(nèi)插處理生成DEM格網(wǎng)來填補(bǔ)ASTER GDEM數(shù)據(jù)的水域空白。常規(guī)的插值算法產(chǎn)生的產(chǎn)品質(zhì)量取決于采樣點(diǎn)的密度，在數(shù)據(jù)稀疏區(qū)會形成假值，不適于水下地形模型的建立[26]。樣條函數(shù)插值算法雖然能通過已知點(diǎn)建立光滑的曲面，但在地形急劇變化地帶，光滑的特性反而掩蓋了地形本來的不連續(xù)特征[27]。目前，張力樣條插值法（Tension Spline）是構(gòu)建海底地形模型的標(biāo)準(zhǔn)插值算法[28-30]，通過張力因子的設(shè)置進(jìn)一步控制了插值超限的現(xiàn)象，也是眾多插值算法中精度最高的。

最后，將新生成的陸地DEM與水下DEM進(jìn)行鑲嵌合并處理，生成水陸一體連續(xù)地形模型。

2 實(shí)驗(yàn)結(jié)果與分析

下面以長江口地區(qū)南支水道為試驗(yàn)區(qū)域，檢驗(yàn)上述方法的實(shí)際效果，在進(jìn)行數(shù)據(jù)融合前，首先利用已有資料將海圖數(shù)據(jù)配準(zhǔn)到ASTER GDEM數(shù)據(jù)，垂直基準(zhǔn)轉(zhuǎn)換至1985國家高程基準(zhǔn)，并將離散化的水深點(diǎn)、等深線按照張力樣條法插值生成水下DEM。

圖3 地形模型構(gòu)建效果Fig.3 Effect of terrain model construction

處理結(jié)果如圖3所示，圖3（a）為ASTER GDEM原始數(shù)據(jù)圖像，長江口地區(qū)受徑流、潮流、波浪和復(fù)雜地形等因素影響，多變的動力條件[31]、復(fù)雜的泥沙輸運(yùn)以及大面積水體造成的鏡面反射效應(yīng)，都影響了光學(xué)傳感器的數(shù)據(jù)質(zhì)量，圖中水面附近的高亮度格網(wǎng)為異常數(shù)據(jù)，圖3（b）是本次試驗(yàn)的融合結(jié)果，原始數(shù)據(jù)中的大量異常數(shù)據(jù)得到有效修復(fù)，深水航道與江中淺灘清晰可辨，區(qū)域內(nèi)高程值集中在46～-46 m之間，與實(shí)際情況相符。

為對本方法的效果進(jìn)行定量分析，采用中誤差（RMSE）為評價指標(biāo)比較校正前后的精度差異，當(dāng)研究區(qū)域地形起伏不大且控制點(diǎn)位分布均勻時，二次曲面擬合法有著較高的精度。如圖4所示，本次實(shí)驗(yàn)收集地面控制點(diǎn)13個，其中10個點(diǎn)進(jìn)行擬合計算求解參數(shù)，另外3個點(diǎn)用作校核計算，驗(yàn)證擬合精度。點(diǎn)位詳細(xì)信息見表2。比較控制點(diǎn)水準(zhǔn)高程與對應(yīng)點(diǎn)位處ASTER像元值，原始數(shù)據(jù)中誤差（RMSE）為6.282 m。經(jīng)過二次曲面擬合計算后的點(diǎn)位信息見表3。擬合計算后的點(diǎn)位中誤差（RMSE）為0.445 m?？梢钥吹交诙吻鏀M合的誤差校正可以有效提高原始數(shù)據(jù)高程精度。

圖4 研究區(qū)域控制點(diǎn)點(diǎn)位分布Fig.4 The distribution of regional control points in study area

表2 控制點(diǎn)點(diǎn)位信息表Table 2 Detail information of control point m

表3 擬合像元值與水準(zhǔn)高程對比Table 3 The comparison between geoid elevation and pixel value with fitting

為了進(jìn)一步觀察海圖與ASTER GDEM數(shù)據(jù)的融合效果，分別在分流口、南支主槽選取貫穿水陸的典型剖面，如圖5所示。由剖面圖（圖5（b）、（d）、（f））可以看出，兩條剖面所處區(qū)域的原始數(shù)據(jù)均有明顯的“跳點(diǎn)”現(xiàn)象，地形多處突變與實(shí)際嚴(yán)重不符?；诤D和ASTER GDEM數(shù)據(jù)融合構(gòu)建的水陸一體連續(xù)地形模型的剖面圖（圖5（a）、（c）、（e））則較真實(shí)地反映了水陸地形，地形剖面變化平緩，水陸交界處無明顯跳躍。結(jié)果表明，本方法有效改善地形數(shù)據(jù)融合過程中存在的高程差異問題，提高ASTER GDEM高程精度，能夠快速構(gòu)建水陸一體連續(xù)地形模型。

圖5 典型剖面地形對比Fig.5 Typical profile topography comparison

3 結(jié)語

以含有準(zhǔn)確高程信息的地物數(shù)據(jù)作為輔助，可在地形數(shù)據(jù)融合過程中提供一定的附加空間信息和限定條件，從而有助于獲得更為準(zhǔn)確的融合結(jié)果。本文采用高精度地面控制點(diǎn)作為輔助信息，以長江口測區(qū)為例，對海圖數(shù)據(jù)與ASTER GDEM數(shù)據(jù)進(jìn)行差異分析與融合試驗(yàn)。實(shí)現(xiàn)了區(qū)域尺度的水陸一體連續(xù)地形模型的快速構(gòu)建。該地形模型可用作布置野外測量任務(wù)的決策依據(jù)，有效克服了河口地區(qū)野外測量工作中存在的測量數(shù)據(jù)冗余、測量目標(biāo)不直觀等缺陷。

本文的融合方法具有原始數(shù)據(jù)易獲取、融合結(jié)果精度高等優(yōu)點(diǎn)，但融合結(jié)果可能與研究區(qū)地形特點(diǎn)相關(guān)。本次研究區(qū)域地形平坦，數(shù)據(jù)波動較小，二次曲面擬合法能取得較好效果，在地形急劇變化的地區(qū)，該方法有待進(jìn)一步驗(yàn)證。此外，由于地面控制點(diǎn)數(shù)量較少，研究區(qū)域相對較大，高程誤差的空間相關(guān)關(guān)系有待進(jìn)一步驗(yàn)證。即便如此，在其他地形數(shù)據(jù)的情況下，基于海圖與ASTER GDEM數(shù)據(jù)融合的方法，可以快速構(gòu)建所在區(qū)域的全面地形信息，從而為相關(guān)研究提供更精確的基礎(chǔ)地形資料。