基于FAST測點(diǎn)的大幅寬HJ星圖像幾何精糾正方法

李全文， 趙衛(wèi)林， 楊曉梅， 劉興權(quán)， 張濤

( 1.中南大學(xué)地球科學(xué)與信息物理學(xué)院，長沙 410083； 2.中國科學(xué)院地理科學(xué)與資源研究所資源與環(huán)境信息系統(tǒng)國家重點(diǎn)實(shí)驗(yàn)室，北京 100101； 3.鄭州市市政勘測設(shè)計(jì)研究院，鄭州 450046)

0 引言

環(huán)境與災(zāi)害監(jiān)測預(yù)報(bào)小衛(wèi)星星座(HJ-1A /1B星，以下簡稱HJ星)是我國第一個(gè)用于環(huán)境與災(zāi)害監(jiān)測預(yù)報(bào)的專業(yè)小衛(wèi)星星座，也是我國又一個(gè)多星多載荷民用對(duì)地觀測系統(tǒng)，目前在生態(tài)環(huán)境監(jiān)測評(píng)估、區(qū)域環(huán)境災(zāi)害監(jiān)測[1]、全球環(huán)境影響分析等多方面得到高度重視和廣泛應(yīng)用。劉睿等[2]對(duì)HJ圖像進(jìn)行的數(shù)據(jù)評(píng)價(jià)表明，整景HJ圖像存在較大的整體幾何誤差，圖像在X和Y這2個(gè)方向的誤差分別超過10和20個(gè)像元，且在2個(gè)方向的誤差沒有規(guī)律可循，在實(shí)際應(yīng)用中需要考慮幾何精糾正問題。

目前，鑒于大幅寬HJ圖像整體幾何誤差復(fù)雜的特點(diǎn)，在以HJ圖像為數(shù)據(jù)源、涉及較高精度的分析應(yīng)用中，多是截取小區(qū)域圖像，在較少的地面控制點(diǎn)(ground control point，GCP)條件下完成幾何精糾正后再進(jìn)行圖像解譯。如劉睿等[3-6]都是先對(duì)截取的小范圍HJ圖像數(shù)據(jù)預(yù)處理后再進(jìn)行分析應(yīng)用的。一般來說，可依照傳統(tǒng)方法選取較少的GCP控制點(diǎn)，采用全局糾正模型完成小范圍HJ圖像的幾何精糾正，使糾正結(jié)果滿足誤差在1～2個(gè)像元的精度要求。但如何保證大幅寬HJ圖像的工程化應(yīng)用，將整景圖像的無規(guī)律幾何形變控制在一定精度范圍內(nèi)，快速實(shí)現(xiàn)整景HJ圖像的幾何精糾正，則是個(gè)研究新課題。

全局糾正模型中，多項(xiàng)式模型簡單易用，一般可以糾正圖像的X，Y方向平移、比例尺變形、旋轉(zhuǎn)和傾斜[7]。采用二次多項(xiàng)式模型對(duì)平坦地區(qū)的IKONOS Geo圖像[8]和IRS-1C圖像[9]做幾何精糾正可得到很好的糾正結(jié)果。采用三次多項(xiàng)式模型對(duì) SPOT HRV[10]，TM[11]和IKONOS Geo[12-14]等衛(wèi)星圖像做幾何精糾正，在地形起伏較大區(qū)域也顯示有好的糾正結(jié)果。但以上糾正模型只適用于覆蓋區(qū)域較小、幾何變形單一的圖像精糾正。有關(guān)試驗(yàn)表明，直接將上述全局糾正模型應(yīng)用于整景HJ星圖像，會(huì)使糾正后的圖像出現(xiàn)某些局部糾正精度比較好、某些局部糾正精度卻很差的不均衡情況。因此，對(duì)于幾何形變復(fù)雜的大幅寬圖像，應(yīng)該選擇局部糾正模型(如選擇基于三角網(wǎng)分割的局部模型)完成圖像幾何精糾正[15]?；贒elaunay三角分割的“橡皮拉伸”(rubber sheeting)局部糾正模型適用于地面控制點(diǎn)豐富的情況，在控制點(diǎn)足夠多的條件下，可以達(dá)到很高的糾正精度[15]。然而地面控制點(diǎn)的分布和選取方式也會(huì)對(duì)圖像精糾正的結(jié)果造成非常大的影響[16]； 而手動(dòng)選點(diǎn)方式過于依賴主觀經(jīng)驗(yàn)判斷，也不適合進(jìn)行大量GCP的選取。如采用人工選點(diǎn)方法對(duì)HJ圖像進(jìn)行幾何精糾正，參照劉盛等[17]得出的精糾正選點(diǎn)密度為0.06～0.09個(gè)/km2的結(jié)論，1景HJ圖像大概需要選擇3 000個(gè)以上GCP，這對(duì)于30 m空間分辨率的HJ圖像而言，在肉眼信息識(shí)別和精度控制上是不具有可行性的。因此，快速獲取大量可靠的GCP及采用局部糾正模型是提高HJ圖像幾何精糾正精度的關(guān)鍵[18]。

在大幅寬圖像的角點(diǎn)檢測方法中，從運(yùn)行效率和精度上考慮，加速分段測試特征(features from accelerated segment test，F(xiàn)AST)算法比尺度不變特征變換(scale-invariant feature transform，SIFT)和最小核值相似區(qū)(smallest univalue segment assimilating nucleus，SUSAN)算法更適合于HJ星大幅寬圖像獲取大量的備選GCP。針對(duì)當(dāng)前算法獲取的大量備選GCP仍存在誤點(diǎn)以及用于誤點(diǎn)篩查的傳統(tǒng)均方根誤差(root-mean-square error，RMSE)閾值選取過于主觀的問題，本文在不依賴數(shù)字高程模型(digital elevation model，DEM)和有理多項(xiàng)式參數(shù) (rational poly-nomial coefficients，RPC)[19]條件下，使用相關(guān)分析方法分析評(píng)估全體候選GCP的質(zhì)量，并完成誤點(diǎn)篩查和獲取可靠的GCP； 最后使用基于Delaunay三角分割的“橡皮拉伸”局部糾正模型對(duì)HJ圖像完成幾何精糾正。另外，本文還建立了一套適合于大幅寬圖像幾何糾正質(zhì)量評(píng)估的多指標(biāo)精度評(píng)價(jià)體系，用于評(píng)價(jià)糾正方法的可靠性和可行性。

1 基于FAST選點(diǎn)的局部幾何精糾正

1.1 FAST自動(dòng)測點(diǎn)與篩查

1.1.1 算法原理

|I(c′)-I(c)|≤t,

(1)

式中：I(c′)為以3像元為半徑的圓形窗口中的點(diǎn)的灰度值；t為閾值。只要滿足式(1)的連續(xù)弧長的數(shù)目≥9，則該中心點(diǎn)為角點(diǎn)。實(shí)際上，I(c′)到I(c)的距離可以不是1，也就是I(c)的“周圍”可以定義為變化的圓或者根據(jù)待糾正圖像的特點(diǎn)將圓半徑設(shè)置大一些，那么當(dāng)I(c′)到I(c)的距離增大時(shí)，圓周線上的所有點(diǎn)不一定能滿足式(1)的條件； 因此實(shí)際糾正時(shí)，需要設(shè)定一個(gè)角點(diǎn)的響應(yīng)值。對(duì)于一個(gè)候選點(diǎn)是否為角點(diǎn)，可以通過定義一個(gè)角點(diǎn)響應(yīng)函數(shù)(corner response function，CRF)來判斷，即

N=∑x?[c(p)]|I(x)-I(p)|>εd，

(2)

式中：I(x)為圓周線上任意一點(diǎn)的圖像灰度值；I(p)為中心像元點(diǎn)的圖像灰度值；p為中心像元點(diǎn)即候選點(diǎn)；εd為給定的一個(gè)極小閾值。通過上述的角點(diǎn)響應(yīng)函數(shù)，可以累加出圓周上滿足式(1)的像元點(diǎn)的個(gè)數(shù)N。如果N大于給定的閾值，就可以確定該候選點(diǎn)為角點(diǎn)。本文初始閾值取17，用此閾值可以較快地排除部分偽角點(diǎn)?？梢酝ㄟ^對(duì)比結(jié)果選擇合適的閾值。完成角點(diǎn)檢測后，即可使用HJ星二級(jí)數(shù)據(jù)已有的初始坐標(biāo)進(jìn)行關(guān)聯(lián)匹配。

1.1.2 誤點(diǎn)篩查

使用FAST算法匹配選點(diǎn)能夠生成大量的GCP(可達(dá)7 000多個(gè))，其中仍可能存在誤差較大的點(diǎn)。在幾何變形性質(zhì)較為單一的情況下，通常根據(jù)某個(gè)糾正模型的RMSE的大小來判斷GCP是否為誤點(diǎn)，對(duì)于小范圍內(nèi)的圖像糾正，該方法能滿足誤點(diǎn)篩查要求。使用FAST算法獲取的大量地面GCP，也可以使用RMSE閾值作為誤點(diǎn)的參考。實(shí)驗(yàn)發(fā)現(xiàn)，通過RMSE閾值選取的GCP只是潛在誤點(diǎn)； 在這些潛在的誤點(diǎn)中，有時(shí)較大的RMSE對(duì)應(yīng)的GCP匹配度可能很高，而較小的RMSE對(duì)應(yīng)的GCP匹配度可能更低。因此，設(shè)置合適的RMSE閾值篩查誤點(diǎn)并反饋FAST匹配選點(diǎn)是保證HJ圖像幾何糾正質(zhì)量的重要環(huán)節(jié)。

FAST算法在自動(dòng)匹配進(jìn)行迭代選點(diǎn)時(shí)也進(jìn)行了初步的偽點(diǎn)篩查，篩查方法基于候選點(diǎn)和鄰域點(diǎn)的差值的統(tǒng)計(jì)閾值進(jìn)行，因此偽點(diǎn)檢測結(jié)果受閾值設(shè)置影響較大，同時(shí)也受圖像質(zhì)量影響。而現(xiàn)有軟件的初次誤點(diǎn)檢測方法是在獲取原始GCP后，當(dāng)被檢驗(yàn)的點(diǎn)和其他大多數(shù)點(diǎn)的精度不能很好吻合時(shí)，該檢驗(yàn)點(diǎn)將被視為誤點(diǎn)被剔除[20]； 但所得到的GCP仍包含不滿足要求的誤點(diǎn)(特別是對(duì)大幅寬、幾何變形較大圖像的處理更是如此)，因此在獲取GCP后仍需要進(jìn)行誤點(diǎn)篩查的后處理工作。對(duì)誤點(diǎn)的篩查往往依賴于個(gè)人主觀判斷，缺乏整體的數(shù)據(jù)規(guī)律分析。以具有代表性的ERDAS軟件的誤點(diǎn)篩查而言，在獲取GCP后，使用全局模型反算GCP的RMSE，通常會(huì)主觀地把大于某個(gè)RMSE值的點(diǎn)作為誤點(diǎn)剔除； 由于該方法對(duì)RMSE的計(jì)算結(jié)果偏優(yōu)[21]，導(dǎo)致誤點(diǎn)剔除的閾值篩選具有很大的不確定性，從而使圖像幾何糾正精度偏低。采用局部模型進(jìn)行糾正時(shí)，ERDAS軟件沒有提供誤點(diǎn)評(píng)價(jià)的方法，故需要借助全局模型進(jìn)行輔助性判斷。對(duì)此，為降低RMSE閾值選擇的主觀性，本文引入誤點(diǎn)和潛在誤點(diǎn)的相關(guān)性分析方法，從數(shù)據(jù)的整體規(guī)律上確定誤點(diǎn)篩查的方法，以此確定合適的RMSE閾值篩查誤點(diǎn)，提高GCP的準(zhǔn)確性。在閾值范圍內(nèi)，實(shí)際誤點(diǎn)占潛在誤點(diǎn)的百分比最大，在保證GCP數(shù)量足夠多且分布均勻的情況下，可以將該RMSE閾值對(duì)應(yīng)的GCP刪除，或者手動(dòng)修改誤差較大的GCP。

顧客滿意是指顧客在接受企業(yè)提供的服務(wù)或者商品之后的消費(fèi)感知。顧客的滿意受到多方位的影響，包括產(chǎn)品的外觀、質(zhì)量、價(jià)格以及銷售人員的服務(wù)態(tài)度等等。顧客的愉悅感越高，滿意度就越高。簡而言之，就是說顧客滿意就是指在消費(fèi)的過程中感覺良好，能夠滿足其消費(fèi)的需求與期望，進(jìn)而將其自身的感受通過直接表達(dá)或暗示的方式對(duì)外傳遞，肯定其消費(fèi)過程的一種方式。而滿意度是消費(fèi)者根據(jù)消費(fèi)體驗(yàn)做出的感知評(píng)價(jià)。顧客滿意是由多方面的因素產(chǎn)生的，不僅包括企業(yè)產(chǎn)品與服務(wù)方面，還包括消費(fèi)者的自身標(biāo)準(zhǔn)。因此要實(shí)現(xiàn)顧客滿意是一項(xiàng)復(fù)雜的工程。

1.2 Delaunay三角網(wǎng)分割

匹配選好同名點(diǎn)并剔除偽點(diǎn)后，即可進(jìn)行幾何糾正。鑒于HJ圖像具有較大的整體幾何誤差，本文采用了Delaunay局部糾正方法。在圖像幾何糾正中，通過Delaunay原則構(gòu)建控制點(diǎn)(GCP)三角網(wǎng)，三角網(wǎng)中每個(gè)三角形由3個(gè)GCP作為三角形頂點(diǎn)。構(gòu)建好三角網(wǎng)后，使用多項(xiàng)式建立原始圖像和目標(biāo)圖像之間的數(shù)學(xué)轉(zhuǎn)換關(guān)系。由于這樣的數(shù)學(xué)轉(zhuǎn)換是通過每個(gè)三角形的3個(gè)GCP獨(dú)立運(yùn)算完成，而不是使用同一個(gè)轉(zhuǎn)換模型對(duì)整景圖像進(jìn)行轉(zhuǎn)換，因此這種局部模型可以很好地對(duì)復(fù)雜圖像進(jìn)行分割糾正，獲得良好的糾正結(jié)果[22]。該方法適合在能夠獲取大量可靠控制點(diǎn)的條件下使用。

Delaunay的特點(diǎn)是，任意一個(gè)三角形外接圓都不包括該三角形之外的頂點(diǎn)(圖1)。這使得Delaunay三角網(wǎng)中的每個(gè)三角形的內(nèi)角盡可能接近等角，從而具有很好的局部適應(yīng)性。

圖1 由13個(gè)點(diǎn)構(gòu)成的Delaunay三角網(wǎng)

最為簡單和常用的三角形內(nèi)部糾正模型是一次多項(xiàng)式模型，即

(3)

該模型不需要額外的參數(shù)，由3個(gè)已知頂點(diǎn)條件即可求解式(3)中3個(gè)未知參數(shù)，完成轉(zhuǎn)換計(jì)算。由于采用“線性橡皮拉伸”(linear rubber sheeting)不一定能得到平滑的糾正結(jié)果，因此根據(jù)需要還可以使用非線性轉(zhuǎn)換模型。

基于FAST算法計(jì)算得到大量的GCP，用GCP構(gòu)建Delaunay三角網(wǎng)，將整景HJ圖像劃分為相對(duì)較小的區(qū)域； 可以認(rèn)為每一個(gè)三角形所覆蓋的區(qū)域地形相對(duì)單一，內(nèi)部幾何形變也相對(duì)簡單[23]，故使用一次多項(xiàng)式校正模型對(duì)每個(gè)三角形進(jìn)行糾正轉(zhuǎn)換。因此，本文采用“線性橡皮拉伸”局部模型對(duì)HJ星圖像進(jìn)行幾何精糾正。

2 精糾正評(píng)價(jià)指標(biāo)

在幾何精糾正后圖像的誤差評(píng)估中，小范圍圖像的幾何誤差相對(duì)簡單，一般使用總RMSE和中誤差評(píng)價(jià)圖像糾正的優(yōu)劣。對(duì)于大幅寬、幾何誤差較大的HJ圖像的糾正結(jié)果，單憑誤差統(tǒng)計(jì)難以表征圖像的誤差分布情況和局部誤差特性，不足以說明糾正后圖像內(nèi)部各區(qū)域誤差的大小和分布都滿足要求。因此，本文在相關(guān)研究[22-25]基礎(chǔ)上，建立了HJ圖像精糾正的精度評(píng)價(jià)指標(biāo)體系，從整體和局部2個(gè)方面評(píng)價(jià)HJ圖像的幾何糾正精度。評(píng)價(jià)方法包括： 誤差統(tǒng)計(jì)、RMSE的空間插值、殘差的標(biāo)準(zhǔn)偏差橢圓和殘差的空間自相關(guān)指數(shù)(Moran指數(shù))。

誤差統(tǒng)計(jì)內(nèi)容包括基于FAST測算得到的GCP、檢核后的GCP、總體中誤差、總RMSE和控制點(diǎn)密度。統(tǒng)計(jì)方法是從整體給出HJ圖像幾何精糾正結(jié)果的定量描述，獲取實(shí)驗(yàn)結(jié)果的總體特征。

對(duì)RMSE進(jìn)行空間插值則可以直觀地掌握HJ圖像幾何糾正的誤差大小和分布范圍，以彌補(bǔ)統(tǒng)計(jì)法對(duì)局部誤差評(píng)估的不足。獲取GCP后，采集了95個(gè)獨(dú)立檢核點(diǎn)，以HJ星整景圖像作為插值范圍，以檢核點(diǎn)的RMSE為插值屬性，使用反距離加權(quán)(inverse distance weighting，IDW)和Kriking方法[26]對(duì)HJ圖像進(jìn)行插值評(píng)估。

誤差統(tǒng)計(jì)從總體上給出對(duì)圖像幾何糾正定量誤差的描述； 插值能較直觀地揭示誤差分布范圍，但得不到誤差偏向的定量描述。使用標(biāo)準(zhǔn)偏差橢圓(standard deviation ellipse)則能度量檢核點(diǎn)的RMSE在空間中的方向特征和離散水平[25]； 如果校正后圖像的RMSE值很小且分布均一，RMSE的標(biāo)準(zhǔn)偏差橢圓的形狀就會(huì)接近標(biāo)準(zhǔn)圓，橢圓的偏轉(zhuǎn)方向應(yīng)該很小[27]。另外，本文使用Moran指數(shù)[28]對(duì)幾何糾正后HJ圖像檢核點(diǎn)的RMSE的鄰近相關(guān)性進(jìn)行評(píng)估，判斷糾正模型對(duì)HJ圖像進(jìn)行幾何精糾正的適合程度。Moran指數(shù)取值范圍為[-1，1](取0值表示不相關(guān)，-1表示強(qiáng)負(fù)相關(guān)，1表示強(qiáng)正相關(guān))。幾何糾正精度很高圖像的RMSE的分布應(yīng)該是彼此獨(dú)立、隨機(jī)分布的，對(duì)應(yīng)的Moran指數(shù)應(yīng)該趨向于0，檢驗(yàn)P值置信度應(yīng)該很低。

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

3.1 實(shí)驗(yàn)圖像

HJ星分為A星和B星。A，B星分別攜載2臺(tái)設(shè)計(jì)原理完全相同的CCD相機(jī)，雙CCD相機(jī)具有相同的幾何一致性。單顆衛(wèi)星的單個(gè)CCD成像幅寬為360 km，雙CCD成像幅寬為710 km[29]。本文實(shí)驗(yàn)是基于福建全省HJ圖像數(shù)據(jù)的工程處理進(jìn)行的，實(shí)驗(yàn)涉及A，B星數(shù)據(jù)。劉睿等[2]結(jié)合TM數(shù)據(jù)對(duì)HJ星的A，B星進(jìn)行的評(píng)估表明，A，B星的1A圖像具有較大的不規(guī)則幾何誤差。本文使用云覆蓋較少的HJ_20111124_454_84_HJ1A-CCD2和HJ_20111224_454_84_HJ1B-CCD2(即A，B星圖像)(圖2)進(jìn)行幾何誤差驗(yàn)證。

圖2 HJ A(上)和B(下)星圖像相對(duì)TM圖像的誤差

HJ星A，B星圖像相對(duì)TM圖像分別有高達(dá)701 m和811 m的異向絕對(duì)誤差，二者相對(duì)TM圖像的誤差都較大。以編號(hào)HJ_20111124_454_84_HJ1B-CCD2的圖像做實(shí)驗(yàn)說明，HJ圖像的覆蓋范圍如圖3所示。

在中尺度應(yīng)用中，可直接應(yīng)用TM圖像進(jìn)行土地利用現(xiàn)狀和動(dòng)態(tài)變化研究[30]。TM圖像的可見光波段和近紅外波段的空間分辨率與HJ圖像的一致，因此以TM圖像作為糾正的基準(zhǔn)圖像。TM圖像幅寬為185 km，遠(yuǎn)小于HJ圖像的710 km，因此，需要使用多景TM圖像拼接出覆蓋HJ圖像的基準(zhǔn)圖像。如圖3所示，基準(zhǔn)圖像用9景TM圖像鑲嵌而成(黃色方框內(nèi)為1景HJ圖像的覆蓋范圍)。

圖3 實(shí)驗(yàn)區(qū)域HJ圖像相對(duì)TM圖像覆蓋范圍

3.2 精度檢核數(shù)據(jù)

橡皮拉伸模型的控制點(diǎn)殘差和RMSE與其他常用糾正模型有所不同。在橡皮拉伸模型中，通過圖像匹配得到的所有GCP點(diǎn)都用于Delaunay三角網(wǎng)的構(gòu)建，三角網(wǎng)上的3個(gè)頂點(diǎn)被認(rèn)為是精確的，即X和Y殘差及RMSE都是0。所以，用橡皮拉伸模型進(jìn)行誤差評(píng)估時(shí)，不能直接用原始GCP的X和Y殘差及RMSE； 需要選擇不用于構(gòu)建Delaunay三角網(wǎng)的同名檢核點(diǎn)，用于糾正結(jié)果的誤差評(píng)價(jià)。實(shí)驗(yàn)中，在HJ圖像和TM參考圖像中分別選取95個(gè)同名點(diǎn)作為檢核點(diǎn)(圖4)。

圖4 HJ圖像糾正檢核點(diǎn)分布

3.3 糾正步驟

幾何糾正的步驟可以分為： 基于FAST選取候選點(diǎn)； 使用多項(xiàng)式和相關(guān)系數(shù)評(píng)估候選點(diǎn)精度，并篩除誤點(diǎn)； 使用線性橡皮拉伸局部糾正模型進(jìn)行糾正，并進(jìn)行精度分析。糾正的關(guān)鍵步驟如圖5所示。

圖5 基于FAST算法的HJ圖像糾正關(guān)鍵步驟

實(shí)驗(yàn)中FAST選點(diǎn)可以用Matlab實(shí)現(xiàn)，也可以參考現(xiàn)有軟件完成，將得到的點(diǎn)導(dǎo)出作為候選點(diǎn)。

3.3.1 FAST選點(diǎn)參數(shù)配置

基于FAST算法的測點(diǎn)要求待匹配的2景圖像之間有足夠的重疊度。本文鑲嵌的TM基準(zhǔn)圖像對(duì)HJ圖像的重疊度達(dá)到100%。此外，F(xiàn)AST測點(diǎn)參數(shù)的設(shè)置很關(guān)鍵，直接會(huì)影響到GCP選點(diǎn)的數(shù)量和質(zhì)量。實(shí)驗(yàn)中，將初始的FAST選點(diǎn)閾值設(shè)置為20； 由于HJ星圖像過大，故將匹配點(diǎn)的搜索增長步長設(shè)置為40，迭代次數(shù)設(shè)置為2。鑒于多個(gè)波段之間的匹配運(yùn)算過于復(fù)雜耗時(shí)，本文只對(duì)HJ圖像和TM 圖像中的1個(gè)波段進(jìn)行匹配，取匹配效果最好的第3波段作為算法的匹配波段，匹配起始位置為1，最小匹配精度設(shè)為0.85，共得到4 635個(gè)候選GCP。GCP分布情況如圖6所示。

圖6 圖像匹配測算得到的GCP分布

3.3.2 誤點(diǎn)篩查

經(jīng)圖像匹配得到的K個(gè)GCP中仍存在不滿足匹配條件的誤點(diǎn)。本文采用相關(guān)分析方法分析RMSE和多項(xiàng)式測算的誤點(diǎn)數(shù)量與真實(shí)誤點(diǎn)數(shù)量的關(guān)系，據(jù)此評(píng)價(jià)所得候選點(diǎn)的整體質(zhì)量，并根據(jù)評(píng)價(jià)結(jié)果修正FAST測點(diǎn)匹配參數(shù)。本文在現(xiàn)有軟件的誤點(diǎn)篩查原理[20]基礎(chǔ)上進(jìn)行改進(jìn)，采用三次多項(xiàng)式的RMSE閾值篩查誤點(diǎn)。以0.1為步長，將RMSE閾值從大到小依次降低，統(tǒng)計(jì)每個(gè)閾值下的潛在誤點(diǎn)數(shù)M和實(shí)際的誤點(diǎn)數(shù)N。將統(tǒng)計(jì)結(jié)果以閾值RMSE為自變量、M和N為因變量，繪制出趨勢圖(圖7)。

圖7 RMSE與潛在誤點(diǎn)數(shù)M和實(shí)際誤點(diǎn)數(shù)N關(guān)系圖

計(jì)算RMSE閾值和誤點(diǎn)數(shù)N的相關(guān)系數(shù)結(jié)果如表1所示。RMSE和N的相關(guān)系數(shù)為-0.097 2，二者呈負(fù)相關(guān)關(guān)系，雙尾檢驗(yàn)的概率值為0(小于0.01，說明在0.01的置信水平下RMSE和N是顯著相關(guān)的)。

① 在0.01的置信水平下RMSE和N是顯著相關(guān)的(雙尾檢驗(yàn))。

從圖7可以看出，RMSE的閾值從最大的8.9～7.7得到的潛在誤點(diǎn)對(duì)應(yīng)的實(shí)際誤點(diǎn)的數(shù)量都是0，這說明大的RMSE閾值不一定能找到匹配不好的GCP。RMSE閾值從7.7～6.4變化范圍內(nèi)，實(shí)際誤點(diǎn)數(shù)量和潛在誤點(diǎn)數(shù)量呈線性增加，表明在該閾值區(qū)間隨著RMSE閾值的降低，潛在誤點(diǎn)M和實(shí)際誤點(diǎn)N都在線性增加。在5.9～5.2的RMSE區(qū)間內(nèi)，M增加得很快，而N基本保持不變。在極限條件下，取RMSE為0，得到的實(shí)際誤點(diǎn)是28個(gè)，對(duì)應(yīng)的潛在誤點(diǎn)數(shù)為4 635個(gè)，但實(shí)際沒有必要在4 635個(gè)點(diǎn)中去尋找可能的28個(gè)誤點(diǎn)。因此，可以取RMSE=5.8作為確定FAST誤點(diǎn)的閾值。實(shí)驗(yàn)中，找到11個(gè)誤點(diǎn)，刪除1誤點(diǎn)，手動(dòng)更正10誤點(diǎn)，總共得到4 634個(gè)用于幾何糾正的GCP。以上分析說明，如果僅以RMSE在8.9～7.7這個(gè)變化范圍內(nèi)找不到誤點(diǎn)而認(rèn)為幾何糾正后圖像中沒有誤點(diǎn)，就會(huì)造成GCP的誤用。以RMSE=5.8作為選擇誤點(diǎn)的閾值較為適宜，經(jīng)過檢核后的GCP可以保證足夠的匹配精度。

根據(jù)圖7獲得的線性方程推算，在極限條件下，RMSE閾值取0時(shí)，得到的真實(shí)誤點(diǎn)數(shù)為28個(gè)； 對(duì)于K=4 635個(gè)候選點(diǎn)而言，誤點(diǎn)率為0.626%(足夠低)，因此原始的算法匹配參數(shù)是可取的，所得到的候選點(diǎn)經(jīng)過修正后可以用于幾何糾正。

3.4 精度評(píng)價(jià)

3.4.1 糾正誤差統(tǒng)計(jì)及其插值分布

對(duì)95個(gè)檢核點(diǎn)分別計(jì)算X和Y方向上的殘差、RMSE和GCP密度，得到的結(jié)果如表2所示。

從表2可以看出，檢核點(diǎn)的總RMSE在1.4個(gè)像元以內(nèi)，表明總體糾正精度非常高。但表2還無法說明圖像糾正結(jié)果的局部RMSE也能滿足要求，因此需要進(jìn)一步對(duì)比分析圖8得到的檢核點(diǎn)RMSE的Kriking和反距離加權(quán)(IDW)插值結(jié)果。

圖8(a)是Kriking插值結(jié)果，色調(diào)越接近白色，糾正誤差RMSE就越大，其值為0.49～1.48個(gè)像元(HJ圖像的像元大小為30 m)。該結(jié)果已經(jīng)很好地滿足了精度要求。從誤差的空間差異分布看，較大的RMSE分布于圖像的邊緣區(qū)域和中間的山地區(qū)域，其大小在可接受范圍內(nèi)。對(duì)比圖8(b)的IDW插值結(jié)果，淺灰色環(huán)及其內(nèi)部區(qū)域的RMSE≥1.5個(gè)像元，最亮區(qū)域及其第一環(huán)狀區(qū)域?yàn)?個(gè)像元≤RMSE≤3個(gè)像元。IDW高亮區(qū)域在整景圖像中僅占很少一部分，說明使用IDW插值和Kriking插值的結(jié)果大部分一致，糾正的誤差分布均勻平緩，初步說明HJ圖像精糾正結(jié)果良好。

(a)RMSE的Kriking插值 (b)RMSE的IDW插值

圖8 檢核點(diǎn)RMSE的插值結(jié)果

3.4.2 誤差的方向性和隨機(jī)性檢驗(yàn)

對(duì)檢核點(diǎn)的X和Y方向殘差計(jì)算標(biāo)準(zhǔn)偏差橢圓，結(jié)果如圖9所示。所得橢圓很接近于標(biāo)準(zhǔn)圓形，說明檢核點(diǎn)殘差在各方向呈同質(zhì)性分布。比照圖中比例尺，橢圓的長、短半軸都在1.5個(gè)像元以內(nèi)。將X和Y殘差疊置在標(biāo)準(zhǔn)偏差橢圓上可以看出，大部分檢核點(diǎn)落在標(biāo)準(zhǔn)橢圓范圍內(nèi)，表明HJ圖像糾正結(jié)果的精度非常高。

圖9 檢核點(diǎn)殘差的標(biāo)準(zhǔn)偏差橢圓及X和Y殘差分布

結(jié)合RMSE在HJ圖像范圍內(nèi)的空間插值結(jié)果可以看出，偏差橢圓表現(xiàn)出的X和Y殘差分布的均一性與圖8(b)表現(xiàn)出的平坦插值結(jié)果非常吻合，說明HJ圖像幾何精糾正的誤差分布滿足應(yīng)用需求，誤差大小也在可接受范圍內(nèi)。

使用Moran指數(shù)分析檢核點(diǎn)的RMSE在空間分布上的相關(guān)性，結(jié)果如表3所示；RMSE的散點(diǎn)圖如圖10(a)所示；RMSE的空間聯(lián)系局部指標(biāo)(local indicators of spatial association，LISA)顯著性分布如圖10(b)所示。

表3 檢核點(diǎn)RMSE的Moran指數(shù)和顯著性檢驗(yàn)

①以P=0.05做顯著性對(duì)比值； ②表示Moran指數(shù)的期望值。

(a) Moran指數(shù)散點(diǎn)圖(b) 空間聯(lián)系局部指標(biāo)顯著水平

從地學(xué)統(tǒng)計(jì)角度而言，如果HJ圖像經(jīng)過幾何精糾正后，糾正誤差很小、精度很高的話，那么檢核點(diǎn)的RMSE在空間上的分布應(yīng)該是隨機(jī)的、彼此獨(dú)立的。表3中的Moran指數(shù)為0.035 2，很接近于0值，說明檢核點(diǎn)的RMSE的空間關(guān)聯(lián)性非常弱。以0.05做顯著性對(duì)比值，全局Moran指數(shù)的P值為0.23，局部的P值為0.29，2個(gè)P值都明顯大于0.05，說明RMSE的空間相關(guān)性顯著水平非常低，即RMSE的空間分布呈不相關(guān)的獨(dú)立分布。在RMSE的顯著性分布圖(圖10(b))的檢核點(diǎn)中，絕大多數(shù)是黑色和藍(lán)色的點(diǎn)，分別對(duì)應(yīng)不顯著或顯著水平非常低的狀態(tài)。全局和局部的P值檢驗(yàn)結(jié)果都表明幾何精糾正的結(jié)果誤差分布均勻、可信度非常高。

以上對(duì)比結(jié)果同時(shí)說明，本文的RMSE閾值選取方法能反映GCP誤點(diǎn)的潛在規(guī)律，提高GCP的客觀代表性。局部的Delaunay三角分割方法能很好地將具有復(fù)雜、無規(guī)律幾何誤差的HJ圖像轉(zhuǎn)變?yōu)閮?nèi)部相對(duì)簡單的三角形，糾正后HJ圖像的X和Y方向的殘差和總RMSE都很小，圖像內(nèi)部沒有出現(xiàn)誤差突變和局部誤差很大的現(xiàn)象，圖像總體糾正誤差非常小，糾正結(jié)果滿足應(yīng)用要求。

4 結(jié)論

大幅寬HJ圖像存在較大的無規(guī)律整體幾何誤差，采用傳統(tǒng)手動(dòng)采點(diǎn)方式無法滿足對(duì)HJ圖像進(jìn)行幾何精糾正的精度需求。本文提出在使用加速分段測試特征(FAST)匹配算法獲取大量GCP、并采用較為客觀的RMSE閾值選取方法篩查GCP的基礎(chǔ)上，使用“線性橡皮拉伸”局部模型對(duì)HJ圖像進(jìn)行幾何精糾正，經(jīng)過殘差散點(diǎn)圖、殘差空間插值、標(biāo)準(zhǔn)偏差橢圓以及Moran指數(shù)等方法的檢驗(yàn)，得到以下結(jié)論：

1)基于FAST匹配測點(diǎn)的局部精糾正方法能獲取足夠的GCP用于大幅寬HJ圖像的幾何精糾正。

2)采用“線性橡皮拉伸”局部模型能克服HJ星大幅寬圖像的整體幾何變形復(fù)雜的問題。

3)使用本文方法對(duì)大幅寬HJ圖像進(jìn)行精糾正后的誤差很小、分布均勻、獨(dú)立性強(qiáng)，能確保圖像精糾正后的精度滿足在1.5個(gè)像元以內(nèi)的要求。

4)本文算法實(shí)現(xiàn)簡單，在需要使用到整景HJ圖像幾何精糾正的情況下，使用本文局部模型糾正方法能夠縮短精糾正的時(shí)間，降低進(jìn)行圖像精糾正的工作強(qiáng)度，減少主觀操作誤差，確保圖像質(zhì)量。

參考文獻(xiàn)(References)：

[1] 鄧鐘,池天河,張新,等.基于HJ星的災(zāi)害特征信息自動(dòng)反演系統(tǒng)[J].計(jì)算機(jī)工程,2012,38(11):231-233.

Deng Z,Chi T H,Zhang X,et al.Disaster characteristics information automatic inversion system based on HJ satellite[J].Computer Engineering,2012,38(11):231-233.

[2] 劉睿,孫九林,王卷樂,等.環(huán)境與災(zāi)害監(jiān)測預(yù)報(bào)小衛(wèi)星CCD數(shù)據(jù)質(zhì)量評(píng)價(jià)[J].地球科學(xué)進(jìn)展,2011,26(9):971-979.

Liu R,Sun J L,Wang J L,et al.Data quality evaluation of Chinese HJ CCD sensor[J].Advances in Earth Sciences,2011,26(9):971-979.

[3] 劉睿,孫九林,張金區(qū),等.中國北方草地覆被的HJ星NDVI校正研究[J].草業(yè)學(xué)報(bào),2011,20(1):189-198.

Liu R,Sun J L,Zhang J Q,et al.HJ NDVI cross-calibration for grassland in Northern China[J].Acta Prataculturae Sinica,2011,20(1):189-198.

[4] 李花,李衛(wèi)國,黃義德.利用HJ星遙感進(jìn)行水稻抽穗期長勢分級(jí)監(jiān)測研究[J].遙感信息,2010,12(6):55-58.

Li H,Li W G,Huang Y D.The monitoring rice growth condition in heading stage by HJ-A/B satellite images[J].Remote Sensing Information,2010,12(6):55-58.

[5] 王瓊,王克如,李少昆,等.HJ衛(wèi)星數(shù)據(jù)在棉花種植面積提取中的應(yīng)用研究[J].棉花學(xué)報(bào),2012,24(6):503-510.

Wang Q,Wang K R,Li S K,et al.Cotton planting area extraction based on multi-temporal HJ remote sensing[J].Cotton Science,2012,24(6):503-510.

[6] Yang F,Juanle W,Pengfei C,et al.Comparison of HJ-1A CCD and TM data and for estimating grass LAI and fresh biomass[J].Journal of Remote Sensing,2012,16(5):1000-1008,1012-1023.

[7] 鐘小明.林區(qū)高空間分辨率遙感影像幾何精校正算法研究[D].西安:西安科技大學(xué),2011.

Zhong X M.Research on geometric rectification of high resolution remote sensing image in forest area[D].Xi’an:Xi’an University of Science and Technology,2011.

[8] Hanley H B,Fraser C S.Geopositioning accuracy of IKONOS imagery:Indications from two dimensional transformations[J].The Photogrammetric Record,2001,17(98):317-329.

[9] Zoej M J V,Mansourian A,Mojaradi B,et al.2D geometric correction of IKONOS imagery using genetic algorithm[J].International Archives of Photogrammetry and Remote Sensing and Spatial Information Sciences,2002,34(B4):1-6.

[10]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.

[11]Pala V,Pons X.Incorporation of relief in polynomial-based geometric corrections[J].Photogrammetric Engineering and Remote Sensing,1995,61(7):935-944.

[12]Ahn C H,Cho S,Jeon J.Ortho-rectification software applicable for IKONOS high resolution images:GeoPixel-Ortho[C]//IEEE 2001 International Geoscience and Remote Sensing Symposium.Sydney,NSW:IEEE,2001:555-557.

[13]Fraser C,Hanley H,Yamakawa T.Three-dimensional geopositioning accuracy of IKONOS imagery[J].The Photogrammetric Record,2002,17(99):465-479.

[14]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 Journal of Photogrammetry and Remote Sensing,2002,57(1/2):24-38.

[15]蔡喜琴,曹建君,蔡迪花,等.中巴地球資源衛(wèi)星CCD影像幾何糾正方法比較[J].遙感技術(shù)與應(yīng)用,2006,21(4):396-398.

Cai X Q,Cao J J,Cai D H,et al.Compare and analysis of CBERS CCD image geometric rectification algorithms[J].Remote Sensing Technology and Application,2006,21(4):396-398.

[16]Wang J H,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.

[17]劉盛,劉慶忠,李國偉.基于RS和DGPS的SPOT5圖像幾何精校正研究[J].安徽農(nóng)業(yè)科學(xué),2009,37(25):12335-12337.

Liu S,Liu Q Z,Li G W.Study on geometric precision correction of SPOT5 images based on RS and DGPS[J].Journal of Amhui Agricultral Sciences,2009,37(25):12335-12337.

[18]Ma L X,Zhou T Y,Xü J.Remote sensing image geometric accurate rectification based on automatic matching and triangulation[J].Journal of Remote Sensing,2011,15(5):927-939.

[19]劉江,岳慶興,邱振戈.RPC校正方法研究[J].國土資源遙感,2013,25(1):61-65.

Liu J,Yue Q X,Qiu Z G.Research on the approach to RPC emendation[J].Remote Sensing for Land and Resources,2013,25(1):61-65.

[20]Leica.Leica Geosystems Geospatial Imaging[M].IMAGINE AutoSyncTMUser’s Guide,2007.

[21]Brovelli M A,Crespi M,Fratarcangeli F,et al.Accuracy assessment of high resolution satellite imagery orientation by leave-one-out method[J].ISPRS Journal of Photogrammetry and Remote Sensing,2008,63(4):427-440.

[22]Leica.Leica Geosystems Geospatial Imaging[M].ERDAS Field Guide,2010.

[23]王新剛.遙感影像幾何校正精度分析[D].北京:中國地質(zhì)大學(xué)(北京),2009.

Wang X G.Geometric correction precision analysis for remote sensing image[D].Beijing:China University of Geosciences(Beijing),2009.

[24]李立鋼,劉波,尤紅建,等.星載遙感影像幾何精校正算法分析比較[J].光子學(xué)報(bào),2006,35(7):1028-1034.

Li L G,Liu B,You H J,et al.The comprehensive comparison of several algorithms for precision rectification of satellite imagery[J].Acta Photonica Sinica,2006,35(7):1028-1034.

[25]王江浩,葛詠.遙感影像幾何校正的GCP殘差模擬分析[J].遙感技術(shù)與應(yīng)用,2011,26(2):226-232.

Wang J H,Ge Y.Simulation analysis of GCP residuals in the remote sensing image registration[J].Remote Sensing Technology and Application,2011,26(2):226-232.

[26]楊成生,張勤,張雙成,等.改進(jìn)的Kriging算法用于GPS水汽插值研究[J].國土資源遙感,2013,25(1):39-43.

Yang C S,Zhang Q,Zhang S C,et al.Research on GPS water vapor interpolation by improved Kriging algorithm[J].Remote Sensing for Land and Resources,2013,25(1):39-43.

[27]Mitchell A.ESRI GIS the ESRI Guide to GIS Analysis:Spatial Measurements and Statistics vol.2.[M].California:ESRI Press,2005.

[28]王勁峰,廖一蘭,劉鑫.空間數(shù)據(jù)分析教程[M].北京:科學(xué)出版社,2010.

Wang J F,Liao Y L,Liu X.Spatial Data Analysis Turtorial[M].Beijing:Science Press,2010.

[29]李傳榮,賈媛媛,胡堅(jiān),等.HJ-1光學(xué)衛(wèi)星遙感應(yīng)用前景分析[J].國土資源遙感,2008,20(3):45-46.

Li C R,Jia Y Y,Hu J,et al.An analysis of the prospects of HJ-1 optical satellites in remote sensing application[J].Remote Sensin for Land and Resources,2008,20(3):45-46.

[30]冉有華,李文君,陳賢章.TM圖像土地利用分類精度驗(yàn)證與評(píng)估——以定西縣為例[J].遙感技術(shù)應(yīng)用,2003,18(2):81-86.

Ran Y H,Li W J,Chen X Z.Verification and assessment of land use classification by using TM image:Taking Dingxi County as an example[J].Remote Sensing Technology and Application,2003,18(2):81-86.