基于加權Voronoi選取控制點的HJ星影像精校正方法

劉興權，李全文

（中南大學地球科學與信息物理學院，湖南 長沙 410083）

0 引言

遙感影像的畸變往往是多種因素綜合影響的結果，包括遙感平臺的姿態穩定性、成像相機、非均質大氣信號的非直線傳輸及地形起伏等因素導致的影像像元的偏移等［1，2］。對此，可以根據不同遙感平臺和影像成像環境選擇不同的校正模型完成幾何精校正工作［3］。如可以采用影像攜帶的星歷參數文件消除衛星平臺等的系統誤差，對于地形起伏較大區域的遙感影像校正可以加入DEM進行正射校正。通常多從數學算法上改進影像校正模型［4－8］以提高影像校正精度，如采用嚴密的數學模型模擬影像的誤差分布情況，進而反算準確的遙感影像。不管是采用描述衛星姿態參數的多項式有理函數參數，還是采用DEM等其他輔助數據進行校正，都離不開GCP（Ground Control Point）的采集［9］，因此如何提高GCP的選取精度對于獲得高質量的遙感影像精校正結果具有重要意義。

GCP的分布范圍、分布密度、分布形態、精度都直接影響影像的校正結果，尤其是GCP的分布形態［10］。許多學者在對遙感影像進行精校正時，都會強調 GCP分布要均勻［2，11］，但是，關于 GCP“均勻分布”卻沒有一種可以度量的描述，均勻性判斷往往依賴于影像校正者的主觀判斷。另外，不同的遙感影像其畸變嚴重程度和畸變性質都不一樣，對于理想的GCP分布，GCP的采樣應能體現不同影像的畸變特征和畸變分布情況，并且能較好地體現影像的主要畸變類型。在對遙感影像進行幾何精校正有關畸變權重因子的研究中，有學者認為不同區域的高程和坡度會影響GCP在臨近區域的代表性，因此從校正區域的DEM數據獲取高程和坡度作為GCP采樣的權重因子設計GCP采樣策略［12］；該研究認為高程和坡度導致GCP不均一性控制范圍，其實質是影像對應區域的地形起伏造成影像像元偏移的結果，因此，解決該問題的一個優化方法可以使用原有的DEM高程數據進行正射校正。還有學者對沿海區域以線性和面狀方式獲取GCP，以研究沒有控制點的近海區域受兩種GCP采樣方式的影響程度［13］，所研究的GCP采樣方式適用于沿海，難以獲取GCP區域的影像幾何精校正，不適合內陸遙感影像。另外，王江浩等將通用Kriging模型融入Polynomial中設計了基于通用Kriging模型的GCP采樣方法［9］，取得了較好的效果。Cheng等［1］在使用 Kriging判斷遙感影像誤差分布趨勢基礎上提出了使用Polynomial各向異性空間建模方法，從改善校正模型角度提高了影像的校正精度。本文借鑒上述方法，在初步評估遙感影像的誤差規律基礎上，從影像自身畸變規律出發，探尋GCP選取的影響因素，并根據影像自身的誤差趨勢提取誤差權重因子和方向權重因子，不依賴于校正影像以外的第三方數據源獲取權重因子。最后，用適當的權重因子構建Voronoi圖，據此采集、修正GCP樣本點完成校正工作。

本文實驗數據是HJ星二級數據產品，以全局校正模型Polynomial多項式作為實驗模型。對于同一個區域，二級HJ星遙感影像和TM影像具有相同的投影坐標。針對HJ星影像質量畸變性質比TM影像畸變性質更為復雜這一問題，本文主要引入加權Voronoi，以初始GCP誤差和方向誤差作為Voronoi權重因子，以加權Voronoi圖元面積比率相等作為度量GCP分布均衡性的依據，據此獲取分布均衡合理的GCP樣本點。本文從理論和方法上給出了GCP“均勻分布”的可度量定義，可為相關的數據處理工作提供有益參考。

1 Voronoi簡介

1.1 Voronoi圖

Voronoi圖可以看做是Delaunay的偶圖，在二維平面上的Voronoi圖是根據已知點集對平面施行的一種分割。其數學定義［14］如下：給定二維平面上n個點構成的集合，n≥3，由所給出的對平面的分割，稱為以Pi（i＝1，2，…，n）為生成元（或母點）的 Voronoi圖，其中d（p，pi）為p和pi間的Euclidean距離區域V（pi）稱為pi的 Voronoi區域，記B（pj，pi）＝｝，則B（pj，pi）就是pi和pj的平均分割線。顯然，Voronoi圖元是由平面上所有到pi的距離比到S中其他點的距離都小的點組成的集合，如線段上的點p到p1和p2的距離是相等的，所有在上側的點p到p1的距離都小于p到其他pi距離。如圖1所示，實線部分得到的是Voronoi圖，虛線部分得到的是Delaunay三角剖分圖，二者為對偶關系，即Voronoi圖可以和Delaunay圖相互轉換。實線包圍的Voronoi圖記為V（p），其具有很好的區域代表性。

圖1 Voronoi圖和Delaunay剖分圖Fig.1 Voronoi diagram and Delaunay tessellation

在Voronoi圖中，每個頂點至少有3條邊與之相連，根據歐拉定理，Voronoi邊和頂點具有e≥3v／2的關系，頂點和Voronoi面具有e≤3n－3的關系［15］，即n個點集合S的Voronoi圖至多有2n－5個頂點和3n－6條邊。結合圖1觀察，可以發現點集S中的點pi的每一個最近臨近點pj可以確定V（pi）的一條邊。

1.2 加權Voronoi圖

普通Voronoi圖中，每個母點位置一樣，假設在二維平面上n個點的集合S＝｛p1，p2，…，pn｝，n≥3，對每個生成元pi（i＝1，2，…，n）賦以非負實數權重wi（i＝1，2，…，n），稱 D（p，pi）＝d（p，pi）／wi為p和pi間的加權距離，稱d（p，pj）｝為點pi的權重為wi的Voronoi區域。將V（pi，wi）（i＝1，2，…，n）及其邊界稱為以pi（本實驗中pi指GCP采樣點）為生成元、wi為權重的加權Voronoi圖［16］。由定義可知，Voronoi圖是在加權Voronoi圖所有生成元的權重均相等時的特殊情況。圖2為相同母點下普通Voronoi圖和加權Voronoi圖，圖2b中wi表示各個母點的權重因子，在實際應用中如何合理確定wi是有效應用加權Voronoi圖的關鍵。

圖2 普通Voronoi圖和加權Voronoi圖Fig.2 The difference between normal Voronoi diagram and weighted Voronoi diagram

2 基于加權Voronoi選取GCP

2.1 實驗數據

實驗所使用到的整景TM遙感影像的編號為lt51210412010014bjc00。從編號為hj＿20111124＿454＿84＿hj1b－ccd2的HJ星影像區域裁取和TM等同區域的影像，兩個遙感影像時間間隔1a，影像對應區域地物的微小變化對矯正方法的驗證影響很小。在中尺度的應用中，可以直接應用TM影像進行信息解譯工作以及進行相關的土地利用現狀分析和土地利用動態變化研究等［17］。TM影像可見光波段和近紅外波段的空間分辨率和HJ星影像的空間分辨率相同，因此，選擇TM影像作為校正實驗基準影像。一景TM影像幅寬為185km，如圖3b所示，裁剪的等幅寬HJ影像如圖3a所示。

圖3 實驗影像數據Fig.3 The RS image data for the study

2.2 構建影像畸變權重因子

應用加權Voronoi選擇均勻GCP樣本點的關鍵是確定能反映遙感影像畸變性質的適當畸變權重。因此，對于使用全局矯正模型的矯正方法，需要考慮不同的遙感影像不同數據級別影像的變形性質，據此確定能反映影像整體變形性質的權重因子。這里根據HJ星影像實驗數據的等級和影像方向畸變的單向性的變形性質，確定用于構建加權Voronoi權重因子，以此獲取具有代表性的GCP樣本。

在遙感影像幾何精校正中，根據校正原理的不同常用到的校正模型可分為兩種：一種是全局校正模型，效果較好的是polynomial多項式模型；另一種是局部校正模型，如Affine模型和Rubber shitting模型。HJ星遙感影像的變形性質評估結果顯示，其幾何畸變往往是以一種或者相似幾種變形為主體，而相對整幅影像其它畸變類型對應的圖像區域所占面積很小，因此，GCP的采樣應能反映影像的整體畸變特征，不應直接按照規則空間格網采樣（Spatial Coverage Sample，SCS）方式進行。SCS的GCP采樣方式相當于默認所有GCP臨近區域的變形性質一樣，實際上不同GCP所在區域的幾何變形性質都有差別，大部分GCP應分布于最能體現遙感影像主要畸變類型的區域，即相對于遙感影像的主要畸變面積占優的區域。因此，在空間分布上，面積占優的影像主體畸變區的GCP可能要比其他畸變性質嚴重，而面積相對很小的區域的分布密度要大。

本實驗使用的全局多項式校正模型不需要描述衛星軌道參數的RPC文件，該模型簡單易用，一般可以糾正影像XY方向平移、比例尺變形、旋轉和傾斜［18］。二維的多項式模型對平坦地區的IKONOS Geo影像［19］和IRS－1C影像［20］做幾何精校正效果很好。三維多項式模型對SPOT－HRV［21］、Landsat－TM［22］和IKONOS Geo影像的幾何精校正［23－25］，在較大地形起伏區域也顯示較好的效果，因此本文實驗選擇多項式模型進行校正實驗。

構建加權Voronoi權重因子，包括原始誤差權重WE和方向誤差權重WD。WE權重定義為基準影像同名點到被校正HJ星影像同名點的歐氏距離，其數學表達如下：

式中：XH和XT分別表示HJ星影像和TM影像同名點的X坐標值，YH和YT分別表示Y坐標值。

TM影像和HJ星二級數據產品具有相同的地理坐標投影，由于HJ星遙感影像相對于TM影像存在整體偏移，因此，WD權重因子定義為TM影像到HJ星影像同名點GCP的矢量連線按照逆時針方向到基準方向夾角的絕對值，這里的基準方向可以看做是經過TM同名原點與笛卡爾坐標系Y軸平行的方向。WD的數學表達為：。

在獲取WD和WE后，可以計算用于構建Voronoi圖的總體權重因子W：

式中：α和β表示加權系數，WiE表示對WE進行Kriging插值后自然分層的第i層面積；類似的，WiD表示對WD進行Kriging插值后自然分層的第i層面積。可以根據初始影像的變形性質確定α和β的大小。本實驗中，HJ星數據為二級數據，已經帶有地理坐標和投影信息，相對于TM基準影像發生了單向偏移。通過對原始獲取的GCP絕對值誤差的Kriging插值可以判定HJ星相對于TM影像存在整體單向的偏移，因此，這里的加權系數取近似值1。另外，為了快速獲取W，也可以使用Erdas的Autosync獲取待校正HJ星影像和TM基準影像的GCP坐標。

2.3 建立加權Voronoi圖

為了獲取WE和WD，開始可以按照SRS（Simple Random Sample）或者SCS方式對GCP進行采樣，這兩種采樣方式獲得的GCP對于遙感影像精校正而言并不是最優的分布，但是二者的通用Kriging插值可以揭示影像潛在的誤差分布趨勢。使用這一原理，可以從GCP誤差分布趨勢中獲取用于構建加權Voronoi圖的誤差權重和方向誤差權重。這里反映GCP均勻分布的標志是用GCP采樣點構建的Voronoi圖的圖元面積之比接近于1，即各個Voronoi圖元的面積比率接近相等。據此可以對SRS或SCS方式采樣獲得的GCP進行調整以獲取均一的Voronoi圖元，操作方法如圖4所示。

圖4 生成均勻GCP方法Fig.4 The method to acquire uniform distribution GCP

GCP的采樣可以在ENVI、Erdas、ArcMap中進行，根據權重模型公式可以使用ArcMap的Python編程計算得到誤差權重因子，如果在ArcMap采樣則需要導入到其他校正軟件中。本實驗在ArcMap中采樣，使用ArcMap中的通用Kriging插值法獲取WE和WD，并通過腳本程序使用柵格數據計算得到W，插值結果如圖5所示。

圖5 GCP絕對誤差和方向誤差權重分布Fig.5 The absolute error and direction error′s distribution map of GCP

在得到W后，即可以根據原始SRS或者SCS采樣結果按權重構建初始Voronoi圖，依據Voronoi圖元面積比率相等原則對原始SRS或SCS采樣獲得的GCP進行修改，直至得到均一的Voronoi圖為止。

3 結果與精度分析

在對比了Polynomial的1次、2次、3次校正結果后，取效果較好的polynomial 3次方模型作為幾何精校正方法的對比實驗。對原始的GCP采樣計算加權Voronoi圖，結果如圖6a所示；按加權Voronoi圖均一性原則對原始GCP進行修正，重新構建Voronoi圖，結果如圖6b所示。

從圖6a可以看到，原始GCP的Voronoi圖元彼此差異非常大，說明每個GCP所控制的影像區域很不均一，不能有效反映影像的全局變形性質，缺乏整體代表性。對比修正后的GCP的Voronoi圖可以發現，每個圖元大小基本一致，對應的每個GCP的控制范圍均一，能較好地反映校正影像的整體變形性質，具有很好的全局代表性，個別圖元稍微偏小，不影響GCP的整體均一性分布。

圖6 GCP加權Voronoi圖Fig.6 The weighted Voronoi diagram of GCP

對比圖6a和圖6b發現，14、23樣本點對應的Voronoi單元比周圍Voronoi單元異常大，靠近邊界的GCP對應的Voronoi單元相對小得多。因此，需要對原始的GCP進行修正，修正后GCP的加權Voronoi圖每個圖元大小基本保持均等。由于在影像中原本屬于理想選取GCP的位置可能會因為缺乏相應的易于識別的地物點，因而只能盡可能移動GCP靠近理想位置，修正GCP的初步工作可以在ArcMap中使用Python腳本完成，之后可以導入到相應的校正軟件中繼續進行HJ星遙感影像幾何精校正。原始GCP對應的Voronoi圖元和修正后圖元的均勻性參數如表1所示。很明顯，修正后的圖元的最小值和最大值比修正前的圖元最小最大值更趨近于均值，修正后圖元的標準差較修正前的Voronoi圖元的標準差小許多，說明修正后的圖元的大小基本相同，總體滿足均一性條件。

表1 GCP對應Voronoi圖元均一性對比Table 1 The contradistinction of Voronoi primitivesbetween original GCP and modified GCP

實驗中使用GCP的RMSE衡量影像的校正精度，為了更好地評價不同GCP選取方法對影像校正精度的影響差異，在覆蓋整幅影像均勻選取16個獨立GCP檢驗點用于精度檢驗。實驗中將本文方法和以SCS采樣的GCP進行影像精校正方法對比（表2）。從表2可以發現，在相同的GCP數量情況下，以SCS方法選取的GCP進行幾何精校正的總RMSE為35.9，而按照加權Voronoi方式選取的GCP對HJ星遙感影像進行幾何精校正的總RMSE為23.4，實驗影像空間分辨率為30m，按30m對應一個像元算，校正精度優于0.8個像元，結果顯示本文方法比SCS方法提高34%的精度。

表2 采樣SCS和加權Voronoi方法進行GCP采樣的影像校正精度Table 2 The HJ image′s correction precision with SCS method and weighted Voronoi method

4 結論

對HJ星遙感影像進行幾何精校正時，在簡單隨機采樣或者空間規則覆蓋采樣基礎上，采用Kriging插值方式獲取初始遙感影像的幾何誤差評價，進而獲取影像誤差權重因子，在此基礎上引入加權Voronoi圖用于評價原始GCP分布的合理性，即在遙感影像數據預處理過程中，對HJ星二級數據影像進行幾何精校正的GCP分布均勻性給出了一個可以度量的描述方式。在有關HJ星影像數據的研究應用中，對影像精校正有較高精度要求時，除了考慮校正模型的適用性和使用額外的輔助數據之外，本文在不依賴其他輔助數據條件下，從GCP的選取方式上給出一個提高HJ星幾何精校正精度的優化方法，使傳統的GCP采樣的校正精度提高了34%左右。本文研究方法可為相關遙感影像的數據預處理工作提供有益的參考。

［1］ CHENG K S，YEH H C，TSAI C H.An anisotropic spatial modeling approach for remote sensing image rectification［J］.Remote Sensing of Environment，2000，73（1）：46－54.

［2］ 楊樹文.薛重生.航片二次幾何校正的應用研究［J］.遙感技術與應用，2002，17（3）：154－157.

［3］ EASTMAN R D，LE MOIGNE J，NETANYAHU N S.Research issues in image registration for remote sensing［J］.Computer Vision and Pattern Recognition，2007，1（1）：1－8.

［4］ HONG G，ZHANG Y.Wavelet－based image registration technique for high－resolution remote sensing images［J］.Computers＆ Geosciences，2008，34（12）：1708－1720.

［5］ YU L，ZHANG D.HOLDEN E J.A fast and fully automatic registration approach based on point features for multi－source remote－sensing images［J］.Computers ＆ Geosciences，2008，34（7）：838－848.

［6］ BUNTING P，LABROSSE F，LUCAS R.A multi－resolution area－based technique for automatic multi－modal image registration［J］.Image and Vision Computing，2010，28（8）：1203－1219.

［7］ LIU X，TIAN Z，CHAI C，et al.Multiscale registration of remote sensing image using robust SIFT features in Steerable－Domain［J］.The Egyptian Journal of Remote Sensing and Space Science，2011，14（2）：63－72.

［8］ LIU X，TIAN Z，DING M.A novel adaptive weights proximity matrix for image registration based on R－SIFT［J］.AEU－International Journal of Electronics and Communications，2011，65（12）：1040－1049.

［9］ WANG J，GE Y，HEUVELINK G B M，et al.Effect of the sampling design of ground control points on the geometric correction of remotely sensed imagery［J］.International Journal of Applied Earth Observation and Geoinformation，2012，18：91－100.

［10］ 李立鋼.星載遙感影像幾何精校正方法研究及系統設計［D］.中國科學院研究生院（西安光學精密機械研究所），2006.

［11］ 張世利，劉健，余坤勇，等.基于ERDAS幾何校正及在閩江流域影像處理中應用［J］.福建林學院學報，2007，27（4）：365－369.

［12］ 張王菲，岳彩榮.徐天蜀.加權Voronoi圖在GCP選取中的應用［J］.測繪科學，2009，34（S2）：195－197.

［13］ 張鷹，張東，顧燕，等.GCP分布對海岸帶TM影像幾何校正精度的影響研究［J］.海洋學報（中文版），2007，29（1）：31－37.

［14］ 馬立玲，張有會.分區加權Voronoi圖的性質及其面積計算［J］.計算機科學，2011，38（2）：195－198.

［15］ SACK J R.URRUTIA J.Handbook of Computational Geometry［D］.Carleton University of Ottawa，1999.

［16］ 趙志輝，李平，黃曉芹.加權Voronoi圖的離散生成［J］.計算機應用與軟件，2007，24（1）：135－136，139.

［17］ 冉有華，李文君.陳賢章.TM圖像土地利用分類精度驗證與評估——以定西縣為例［J］.遙感技術與應用，2003，18（2）：81－86.

［18］ 鐘小明.林區高空間分辨率遙感影像幾何精校正算法研究［D］.西安科技大學，2011.

［19］ HANLEY H，FRASER C.Geopositioning accuracy of IKONOS imagery：Indications from two dimensional transformations［J］.Photogramm Rec.，2001，17（98）：317－329.

［20］ VALADAN Z M，FOOMANI J.Mathematical modelling and accuracy testing of IRS－1Cstereo pairs［C］.Joint Workshop of ISPRS WG 1／1，1／3and 1V／4，Sensors and Mapping from Space，Hanover，1999.

［21］ OKAMOTO A，FRASER C，HATTORL S，et al.An alternative approach to the triangulation of SPOT imagery［J］.International Archives of Photogrammetry and Remote Sensing，1998，32（4）：457－462.

［22］ PALA V，PONS X.Incorporation of relief in polynomial－based geometric corrections［J］.Photogramm Eng.Remote Sens.，1995，61（7）：935－944.

［23］ AHN C，CHO S，JEON J.Ortho－rectification software applicable for IKONOS high resolution images：GeoPixel－Ortho［C］.Geoscience and Remote Sensing Symposium，2001.

［24］ FRASER C，HANLEY H，YAMAKAWA T.Three－dimensional geopositioning accuracy of IKONOS imagery［J］.Photogramm Rec.，2002，17（99）：465－479.

［25］ VASSILOPOULOU S，HURNI L，DIETRICH V，et al.Orthophoto generation using IKONOS imagery and high－resolution DEM：A case study on volcanic hazard monitoring of Nisyros Island（Greece）［J］.ISPRS－J Photogramm Remote Sens.，2002，57（1－2）：24－38.