張麗

(上海市測(cè)繪院，上海 200129)

1 引 言

地表形變是在地殼表層發(fā)生的一種空間位移變化，由于地殼運(yùn)動(dòng)、火山運(yùn)動(dòng)、地震活動(dòng)等自然因素，或者是地下水開采、礦山開采、城市建設(shè)等人為因素引發(fā)。近幾年地震頻發(fā)，有專家預(yù)測(cè)全球開始進(jìn)入板塊運(yùn)動(dòng)活躍期，因此，對(duì)地表進(jìn)行變形監(jiān)測(cè)顯得尤為重要。

長(zhǎng)期以來(lái)，研究者主要依靠傳統(tǒng)的地面水準(zhǔn)測(cè)量和GPS監(jiān)測(cè)網(wǎng)來(lái)進(jìn)行監(jiān)測(cè)，其監(jiān)測(cè)精度雖可高達(dá)亞毫米級(jí)，但不適宜大區(qū)域的長(zhǎng)期重復(fù)監(jiān)測(cè)，且要消耗大量人力物力。近年來(lái)，隨著SAR技術(shù)的發(fā)展，其全天時(shí)、全天候工作，且可以大范圍、高精度地從空間獲取地形信息和地表微小的形變信息的優(yōu)勢(shì)彌補(bǔ)了這一缺陷；然而該技術(shù)容易受相位去相關(guān)影響，在形變過(guò)大區(qū)域常常無(wú)法監(jiān)測(cè)，且只能獲取沿雷達(dá)視線向即距離向的地表形變。

而光學(xué)影像匹配技術(shù)是借助遙感影像的幅度信息，采用互相關(guān)技術(shù)尋找兩幅圖像的同名點(diǎn)，從而獲得同名點(diǎn)處對(duì)應(yīng)像素在方位向和距離向的偏移量。由于避免了相位運(yùn)算，該方法能很好地彌補(bǔ)了受相位去相關(guān)影響這一不足。因此，本文將該方法應(yīng)用到具體地表形變監(jiān)測(cè)案例中，并研究其可行性。

2 光學(xué)影像匹配

2.1 光學(xué)影像匹配技術(shù)原理

影像匹配技術(shù)的根本就是圖像配準(zhǔn)，在光學(xué)影像里，圖像配準(zhǔn)也稱為幾何精校正，它規(guī)避了成像的空間幾何過(guò)程，直接利用地面控制點(diǎn)(GCP，Ground Control Points)數(shù)據(jù)，對(duì)衛(wèi)星影像的幾何畸變過(guò)程進(jìn)行數(shù)學(xué)模擬，利用已知條件確定模型中的位置參數(shù)，建立原始的畸變圖像空間與校正空間之間的相互對(duì)應(yīng)關(guān)系，然后利用這種對(duì)應(yīng)關(guān)系把畸變圖像空間中的全部元素變換到校正圖像空間中去，從而實(shí)現(xiàn)幾何精校正。精校正后，光學(xué)遙感影像間的幅度信息仍存在差異，通過(guò)傅立葉變換，在頻率域內(nèi)采用高精度的相位相關(guān)技術(shù)即可提取偏移量。

2.2 光學(xué)影像匹配技術(shù)流程

透過(guò)原理可以看出，光學(xué)影像匹配技術(shù)主要分為正射校正和相位相關(guān)計(jì)算兩步，本文所采用的數(shù)據(jù)處理流程大致包括：

①利用外部DEM數(shù)據(jù)及影像的太陽(yáng)高度角，方位角模擬生成DEM陰影圖；

②從震前影像與DEM陰影圖上選取同名點(diǎn)，對(duì)震前影像進(jìn)行正射糾正；

③從震后影像與已校正的震前影像上選取同名點(diǎn)，對(duì)震后影像進(jìn)行正射糾正；

④對(duì)已校正后的震前震后影像做相位相關(guān)計(jì)算，獲取有效形變值。

(1)選取同名點(diǎn)

在上述兩次配準(zhǔn)過(guò)程中均涉及同名點(diǎn)的選取，通常選取同名點(diǎn)分為人工選點(diǎn)和自動(dòng)選點(diǎn)兩種方法。本文采取人工選點(diǎn)法，即通過(guò)目視在兩幅待配準(zhǔn)影像上，手動(dòng)選擇若干同名點(diǎn)。同名點(diǎn)選取時(shí)需要遵循以下原則：

①至少選3個(gè)點(diǎn)，若待糾正影像和一幅色彩較暗的DEM配準(zhǔn)時(shí)需要選15～30個(gè)點(diǎn)；

②同名點(diǎn)要連續(xù)選擇，但需遠(yuǎn)離形變區(qū)域，若整幅圖都有形變就在形變最小處選點(diǎn)；

③點(diǎn)選取時(shí)不能太靠近圖像邊界；

④同名點(diǎn)周圍最好有容易識(shí)別的地物，最好選擇無(wú)須校正的人工物。

(2)影像正射校正與重采樣

正射校正是對(duì)圖像空間和幾何畸變進(jìn)行校正生成多中心投影平面正射圖像的處理過(guò)程。它除了能糾正一般系統(tǒng)因素產(chǎn)生的幾何畸變外，還可以消除地形引起的幾何畸變。采用少量的地面控制點(diǎn)與相機(jī)或衛(wèi)星成像模型相結(jié)合，確立相機(jī)或傳感器平臺(tái)、圖像和地面三者之間的空間位置關(guān)系，建立正確的校正公式，產(chǎn)生精確的正射影像。

影像重采樣是影像數(shù)據(jù)重新組織過(guò)程中的灰度處理方法，影像采樣是按一定的間隔采集影像灰度數(shù)值，當(dāng)閾值不位于采樣點(diǎn)上的原始函數(shù)的數(shù)值時(shí)，就需要利用已采樣點(diǎn)進(jìn)行內(nèi)插，稱為重采樣。本文所采用的重采樣法為三次卷積內(nèi)插法(也叫辛克插值法)，使用內(nèi)插點(diǎn)周圍的16個(gè)像元值進(jìn)行距離加權(quán)計(jì)算柵格值，先在Y方向內(nèi)插四次(或X方向)，再在X方向(或Y方向)內(nèi)插四次，最終得到該像元的柵格值。

(3)相位相關(guān)法估算偏移量

在經(jīng)過(guò)正射校正和重采樣后，兩幅影像應(yīng)該具有相同的分辨率，統(tǒng)一的坐標(biāo)系，但因?yàn)榇嬖诘乇硇巫?，兩幅影像并不能完全重合。因此，還需要利用COSI-CORR軟件模塊提供的高精度相位相關(guān)算法進(jìn)行偏移量估計(jì)。該模塊包含了兩種算法，一種是基于頻率域(Frequential)變換的方法，另一種是基于統(tǒng)計(jì)學(xué)(Statistical)理論的方法。相位相關(guān)法是用于配準(zhǔn)圖像進(jìn)行平移變換的典型方法，它主要是依賴于傅立葉偏移理論。

假設(shè)有兩幅影像i1和i2，用(△x，△y)來(lái)代表兩幅影像上對(duì)應(yīng)點(diǎn)的形變差異，則有公式：

i2(x，y)=i1(x-△x，y-△y)

(1)

用I1和I2分別代表兩幅影像的傅立葉變換，則有以下關(guān)系式：

I2(wx，wy)=I1(wx，wy)e-j(wx△x+wy△y)

(2)

其中wx和wy分別代表行和列的頻率變化。定義兩幅影像的互能量譜為如下：

(3)

上式中，C(wx，wy)是i1(x，y)和i2(x，y)之間的互能量譜，*代表復(fù)數(shù)共軛相乘。

影像間的相對(duì)形變量可以通過(guò)互能量譜的二維斜率來(lái)計(jì)算，因此對(duì)式(3)進(jìn)行傅立葉逆變換：

F-1{ej(wx△x+wy△y)}=δ(wx△x+wy△y)

(4)

而在應(yīng)用中，C(wx，wy)僅為兩幅影像間理論上的互能量譜，由實(shí)際計(jì)算得到的互能量譜為Q(wx，wy)。定義一個(gè)目標(biāo)函數(shù)為φ(△x，△y)：

(5)

其中，W(wx，wy)表示加權(quán)矩陣。目標(biāo)函數(shù)φ(△x，△y)在其他位置都為零，只有在平移位置處不為零，因此這個(gè)位置就是影像間的相對(duì)偏移，當(dāng)目標(biāo)函數(shù)φ(△x，△y)達(dá)到最小值時(shí)的(△x，△y)就是所求的偏移量。

3 實(shí)驗(yàn)數(shù)據(jù)處理

北京時(shí)間2010年4月14日，在我國(guó)青海省玉樹縣發(fā)生了震級(jí)為里氏7.1級(jí)的大地震。此次地震的發(fā)震構(gòu)造是青藏高原東部規(guī)模巨大的甘孜-玉樹斷裂帶，該斷裂帶長(zhǎng)度約 500 km，總體走向呈285°～315°弧形展布，是一條大型左旋走滑斷裂帶。

3.1 實(shí)驗(yàn)數(shù)據(jù)參數(shù)

本實(shí)驗(yàn)選用了兩幅SPOT5衛(wèi)星影像。影像數(shù)據(jù)采集的時(shí)間分別為2009年11月5日和2010年4月15日，影像數(shù)據(jù)具體參數(shù)如表1所示。選用時(shí)間相距較近的衛(wèi)星影像可以最大限度排除人工建筑等因素的影響，協(xié)調(diào)時(shí)間、太陽(yáng)方位角、入射角的一致，確保了不會(huì)由于角度問題導(dǎo)致正射校正的偏差。

實(shí)驗(yàn)選用的SPOT5衛(wèi)星影像參數(shù) 表1

續(xù)表1

選用的SRTM數(shù)據(jù)絕對(duì)高程精度為 16 m，相對(duì)高程精度為 10 m左右，坐標(biāo)系為WGS84地理坐標(biāo)系，投影方式為UTM投影。

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

本文利用基于ENVI平臺(tái)開發(fā)的COSI-CORR軟件模塊，按照上述技術(shù)路線進(jìn)行數(shù)據(jù)處理。

(1)生成DEM陰影圖

由于已有SRTM數(shù)據(jù)與SPOT5影像分辨率存在較大差異，數(shù)據(jù)處理的第一步是先將DEM插值到空間分辨率 10 m，然后根據(jù)太陽(yáng)高度角和太陽(yáng)方位角生成一幅DEM陰影圖，將其作為參考影像，用于震前影像的正射校正中，如圖1所示。

(2)SPOT5影像匹配

兩次配準(zhǔn)時(shí)，在SPOT影像和DEM陰影圖上，通過(guò)人工選點(diǎn)法分別選取了18個(gè)和19個(gè)均勻分布的同名點(diǎn)。在選取過(guò)4個(gè)及以上同名點(diǎn)后，軟件模塊會(huì)自動(dòng)預(yù)測(cè)后續(xù)同名點(diǎn)位置，為提高配準(zhǔn)精度，在經(jīng)過(guò)重新選取和刪除誤差較大的點(diǎn)后，保證兩次的均方根誤為分別為0.48和0.43，均小于1個(gè)像元。

圖1 生成的DEM陰影圖

選取同名點(diǎn)后，需將其轉(zhuǎn)換為GCP點(diǎn)并進(jìn)行優(yōu)化，才能用于后續(xù)校正與重采樣。而COSI-CORR軟件模塊提供兩種優(yōu)化方法(相位相關(guān)法和統(tǒng)計(jì)相關(guān)法)，可以使同名點(diǎn)間配準(zhǔn)誤差值達(dá)到分米級(jí)甚至是厘米級(jí)。第一次正射校正時(shí)，由于參考影像為DEM陰影圖，其影像特點(diǎn)與SPOT影像存在很大差異，因而選用統(tǒng)計(jì)相關(guān)法。而第二次正射校正時(shí)，由于參考影像為校正過(guò)后的震前影像，兩幅影像特點(diǎn)相一致，因此選用相位相關(guān)法來(lái)提高計(jì)算精度。兩次優(yōu)化的結(jié)果表明，選取5次迭代計(jì)算能夠使誤配準(zhǔn)值即殘差偏移量達(dá)到收斂。點(diǎn)優(yōu)化過(guò)程的迭代結(jié)果如圖2所示。

圖2 5次迭代后GCP點(diǎn)處的殘差偏移量均值和標(biāo)準(zhǔn)差

(3)影像校正與重采樣

當(dāng)GCP點(diǎn)優(yōu)化結(jié)果達(dá)到收斂且滿足精度要求后，影像進(jìn)行正射校正與重采樣。此模塊采用的正射校正模型為嚴(yán)格軌道模型(Pushroom Sensor)，校正過(guò)程利用頭文件中衛(wèi)星狀態(tài)參數(shù)提高幾何校正的精度。在選用精確度最高的辛克插值法進(jìn)行重采樣后，得到校正后的震前、震后影像如圖3和圖4所示。

圖3 校正后的震前影像

圖4 校正后的震后影像

(4)相位相關(guān)計(jì)算偏移量

本文經(jīng)過(guò)反復(fù)對(duì)比，進(jìn)行相位相關(guān)法計(jì)算時(shí)，選用滑動(dòng)窗口為64×64的單一型窗口，即X和Y方向的Initial和final都選擇64；滑動(dòng)窗口在影像X和Y方向移動(dòng)的像素?cái)?shù)都采用默認(rèn)值16；迭代次數(shù)選擇2，掩模閾值選擇0.9。獲得了此次地震東西向和南北向的同震位移量，結(jié)果如圖5所示。此一步可輸出偏移量文件，共包含三個(gè)波段數(shù)據(jù)，分別為地表東西向形變值，南北向形變值及信噪比SNR值。

圖5玉樹地震同震位移圖(左為東西向形變，右為南北向形變)

圖6 去噪后玉樹地震局部同震位移分布圖

(5)形變結(jié)果優(yōu)化

通過(guò)對(duì)結(jié)果的SNR值(0效果差，1效果優(yōu))分析，得到了同震位移場(chǎng)的質(zhì)量評(píng)價(jià)。圖6中，上圖為震前SPOT5衛(wèi)星影像，紅線為同震地表破裂帶，中圖為東西向的同震位移分布圖，下圖為南北向的同震位移分布圖。從圖中可以看到，雖有一些效果比較差的地方，但主要分布在遠(yuǎn)離地表破裂帶處，在破裂帶展布區(qū)域分布很少。這是由于玉樹地處高原地區(qū)，地形起伏比較大，在一些地區(qū)得到的數(shù)據(jù)結(jié)果不是很理想，而且由于震后SPOT5衛(wèi)星影像為第二天應(yīng)急拍攝，沒有考慮云量，導(dǎo)致震后影像云量比較多，但幸運(yùn)的是云量基本沒有覆蓋地表破裂帶的范圍，而且地形起伏導(dǎo)致的不理想數(shù)據(jù)結(jié)果也很少在地表破裂帶上，這保證了得到的同震位移場(chǎng)至少在近地表破裂帶處是比較可信的。

從圖中可以明顯地看出地表破裂帶的展布以及同震位移量的分布。此次玉樹地震是左旋走滑型地震，產(chǎn)生的同震位移主要為水平位移，垂直位移分量比較小。其最大位移不超過(guò) 2 m，地表裂長(zhǎng)度為 30 km左右。

4 結(jié) 語(yǔ)

本文對(duì)光學(xué)影像匹配技術(shù)進(jìn)行了初步研究，利用光學(xué)影像空間分辨率高的特點(diǎn)，獲取了玉樹地震水平方向的形變場(chǎng)。根據(jù)野外地表破裂帶調(diào)查，此次地震較清晰的地表破裂帶由3條主破裂左階組成，總長(zhǎng)約 31 km，左旋走滑性質(zhì)，最大左旋位移量約為 1.8 m。這與本文得到的破裂帶長(zhǎng)度和最大同震位移量基本一致，這也從另一個(gè)角度說(shuō)明了實(shí)驗(yàn)得到的同震位移場(chǎng)的合理性和正確性，顯示了光學(xué)影像匹配技術(shù)在監(jiān)測(cè)地表形變方面的巨大應(yīng)用價(jià)值。