張可欣 田崇瑞 焦 健 曾琪明

1 北京大學(xué)地球與空間科學(xué)學(xué)院，北京市頤和園路5號(hào)，100871

InSAR的精度會(huì)受到電離層效應(yīng)的影響，需要進(jìn)行電離層校正。目前常用的校正方法為距離向頻譜分割法RSS[1]，該方法基于電離層的色散特性區(qū)分電離層延遲相位和其他非色散相位，具有清晰的物理機(jī)制，已成功應(yīng)用于GNSS雙頻觀測(cè)電離層校正中[2]。然而RSS精度受相位噪聲影響較大，在實(shí)際處理中需要通過多視和低通濾波(下文統(tǒng)稱為平滑處理)去除高頻噪聲、提高電離層估算精度。但不同平滑參數(shù)下獲取的RSS電離層相位差異較大，且RSS方法計(jì)算的電離層相位是主、從影像成像時(shí)刻電離層狀態(tài)的差值，難以通過相位條紋形態(tài)判斷其可靠性，需要處理人員依據(jù)經(jīng)驗(yàn)進(jìn)行判斷，缺乏客觀依據(jù)。

針對(duì)上述問題，本文擬選用國(guó)際參考電離層IRI模型計(jì)算的電離層分布為參考，輔助確定RSS方法估算電離層相位的平滑處理參數(shù)。重點(diǎn)探討無明顯地表形變區(qū)和緩慢形變區(qū)RSS方法校正電離層時(shí)平滑處理參數(shù)的選取策略，以期推廣到因地震或其他原因引發(fā)的強(qiáng)烈形變區(qū)RSS方法校正電離層時(shí)的平滑處理中。

1 方法與原理

1.1 RSS估算電離層相位

對(duì)于重軌InSAR，在不考慮相位噪聲的情況下，其干涉相位成分可以表示為：

ΔφInSAR=(Δφdisp+Δφgeo+Δφtropo)+Δφiono

(1)

式中，Δφdisp為地表形變相位，Δφgeo為包含平地和地形的幾何相位，Δφtropo為對(duì)流層延遲相位，上述相位統(tǒng)稱為非色散相位；Δφiono為色散的電離層延遲相位。RSS方法根據(jù)電離層的色散特性估算得到電離層相位[1]：

(2)

式中，fH和fL分別為被分割的高、低頻譜子頻帶影像的中心頻率，f0為全頻帶影像的中心頻率，ΔφH和ΔφL分別為高、低頻率的子頻帶干涉相位。最后，將電離層從全頻帶InSAR相位中減去，即可得到電離層校正后的InSAR相位：

(3)

由于現(xiàn)有的在軌衛(wèi)星帶寬較窄，頻譜分割后生成的子頻帶干涉圖中含有較多噪聲，因此通過多視和低通濾波去除相位噪聲對(duì)RSS估算精度的提升有重要作用[3]。然而多視視數(shù)和濾波窗口的增加會(huì)導(dǎo)致較小空間尺度的電離層變化無法被反演，獲取更高估算精度和保留更多電離層細(xì)節(jié)之間的矛盾導(dǎo)致平滑參數(shù)的選取困難且主觀。

1.2 IRI模型計(jì)算差分電離層相位

IRI模型是由國(guó)際空間研究委員會(huì)(Committee on Space Research，COSPAR)與國(guó)際無線電科學(xué)聯(lián)合會(huì)(International Union of Radio Science，URSI)聯(lián)合提出的經(jīng)驗(yàn)型電離層預(yù)報(bào)模型，采用最新的IRI 2016開源軟件可計(jì)算得到指定時(shí)間、地點(diǎn)的天頂方向電子總含量(certical total electron content, VTEC)[4]。由于SAR為側(cè)視成像，根據(jù)式(4)可將VTEC轉(zhuǎn)化為SAR影像中某一位置處的電離層延遲相位[1]：

(4)

式中，f為載波頻率，c為光速，K為常數(shù)40.28 m3/s2，θ為SAR側(cè)視角。IRI模型程序只能計(jì)算單點(diǎn)的電離層VTEC，全景SAR影像范圍內(nèi)的電離層延遲相位可通過循環(huán)迭代得到。

IRI模型時(shí)空分辨率不高，但能客觀反映一般情況下區(qū)域內(nèi)大尺度電離層空間的變化趨勢(shì)和數(shù)量級(jí)。由于高頻空間噪聲會(huì)干擾對(duì)InSAR反演的電離層整體分布特征的判斷，因此IRI模型可為InSAR提供客觀參考，便于相互印證。

2 實(shí)驗(yàn)區(qū)選擇與結(jié)果分析

SAR干涉數(shù)據(jù)選用L波段ALOS-2 PALSAR的條帶模式數(shù)據(jù)。實(shí)驗(yàn)區(qū)優(yōu)先選擇電離層環(huán)境相對(duì)簡(jiǎn)單的中、低緯度區(qū)域以及無形變或緩慢形變區(qū)。有研究表明，地震[6]、地下核爆[7]等強(qiáng)烈地表形變會(huì)對(duì)電離層產(chǎn)生異常擾動(dòng)，而IRI模式主要反映的是平靜電離層的平均狀態(tài)，對(duì)瞬時(shí)電離層擾動(dòng)反映不足[8]。在不考慮對(duì)流層延遲及地形等其他相位時(shí)認(rèn)為，無地表形變區(qū)InSAR相位變化主要是由電離層引起的，故通過研究電離層校正前后InSAR相位標(biāo)準(zhǔn)差變化可評(píng)價(jià)校正效果。為盡量減少干擾電離層和其他InSAR 相位因素的影響，將實(shí)驗(yàn)區(qū)選擇在低緯度無顯著形變的巴西阿克里州(Acre)，以及中緯度有緩慢形變的中國(guó)甘肅省定西市。綜合考慮多視對(duì)低相干噪聲的抑制效果與高相干區(qū)域的保留效果后，選擇方位向(az)和距離向(rng)多視視數(shù)為16×8和32×16(下文均表示為az×rng)；采用高斯濾波的方法，選擇方位向和距離向尺寸均為200、400、600、800像素的濾波窗口，獲得相對(duì)平滑的電離層條紋。

2.1 無明顯形變研究區(qū)

實(shí)驗(yàn)區(qū)為巴西阿克里州東南角，位于世界上地震活動(dòng)最少的穩(wěn)定大陸區(qū)之一的南美洲中部板塊內(nèi)部[9]。主從影像成像時(shí)間為2014-08-13和2014-08-27，時(shí)間基線為14 d，研究時(shí)間段內(nèi)無地震報(bào)道，因此該區(qū)域地表形變可忽略不計(jì)。數(shù)據(jù)模式為High sensitive，數(shù)據(jù)帶寬42 mHz。

該實(shí)驗(yàn)區(qū)主要處理結(jié)果如圖1所示。圖1(a)為全頻帶InSAR相位圖，其相位主要由電離層延遲、對(duì)流層延遲和殘余地形相位等組成。圖1(c)為32×16多視下RSS方法反演出的原始相位分布，可以看出，與圖1(b)IRI計(jì)算的相位分布相差甚遠(yuǎn)，存在較多的空間高頻信號(hào)，不能明顯反映電離層的變化趨勢(shì)，需要對(duì)其進(jìn)行平滑處理。

圖1 阿克里州地區(qū)處理結(jié)果Fig.1 The processing results of Acre

不同多視視數(shù)和濾波窗口處理后的RSS電離層相位分布如圖2所示。RSS與IRI電離層相位在形態(tài)上較為相似，但不同平滑參數(shù)的RSS電離層相位在細(xì)節(jié)上有所不同。當(dāng)多視視數(shù)及濾波窗口較小時(shí)，RSS電離層存在較多細(xì)小結(jié)構(gòu)，而隨著多視視數(shù)和窗口尺寸的增加，RSS電離層條紋更加平直，與IRI電離層形態(tài)更加接近。多視和濾波共同起到了平滑電離層的作用，但二者對(duì)電離層條紋延伸趨勢(shì)的影響各有不同。對(duì)于多視16×8(az×rng)的RSS電離層相位，即使選用較大濾波窗口，仍無法獲得相對(duì)平滑的電離層條紋。因此，在平滑參數(shù)的選擇中，足夠數(shù)量的多視視數(shù)是獲取理想結(jié)果的必要條件。從目視效果上看，當(dāng)多視視數(shù)為32×16、濾波窗口達(dá)到400像素及以上時(shí)，RSS電離層條紋和IRI電離層條紋基本一致，即對(duì)于沒有更多電離層擾動(dòng)的無明顯形變區(qū)域，上述參數(shù)能取得較好的平滑效果。

由圖2可見，電離層在空間分布上主要隨緯度變化，因此分別在IRI和RSS電離層相位中心經(jīng)線和中心緯線處繪制剖面線進(jìn)行分析。不同平滑參數(shù)下阿克里州RSS及IRI電離層相位剖面線如圖3所示(圖例中數(shù)字格式為“方位向視數(shù)與距離向視數(shù)_濾波窗口”)。由圖3(a)和(b)可見，整體上所有RSS電離層相位(藍(lán)色和橘色線)都表現(xiàn)出與IRI電離層相位(紅色粗線)較高的相似性，但I(xiàn)RI與RSS存在一個(gè)整體偏差，可能是因?yàn)镮RI差分電離層相位為絕對(duì)值而RSS電離層相位為相對(duì)值。因此比較IRI與RSS電離層相位差異的重點(diǎn)是其分布形態(tài)而非具體數(shù)值。數(shù)值變化范圍上，RSS沿經(jīng)線方向變化大，在研究范圍(約0.7°)內(nèi)約6 rad；沿緯線方向變化小，不足1 rad，因此后續(xù)主要沿經(jīng)線方向進(jìn)行討論。圖3(c)展示了放大后不同平滑參數(shù)下RSS電離層相位剖面線，整體來看，不同平滑參數(shù)RSS電離層相位較為接近。當(dāng)多視視數(shù)及濾波窗口較小(如3216_200)時(shí)，曲線仍有明顯波動(dòng)，但值域范圍較小。

圖3 不同平滑參數(shù)下阿克里州RSS及IRI電離層相位剖面線Fig.3 RSS and IRI ionospheric phase profile lines with different smoothing parameters in Acre

為定量比較不同參數(shù)RSS電離層與IRI的相似程度，計(jì)算電離層相位變化范圍較大的中心經(jīng)線上IRI與RSS差值的均方根偏差RMSD。選取中心經(jīng)線上相干系數(shù)最高時(shí)(上標(biāo)h表示)RSS與IRI電離層相位的差Δφconst為參考值，以減少該參考點(diǎn)自身誤差及解纏錯(cuò)誤的影響：

(5)

RMSD表示為：

(6)

綜合考慮目視效果、RMSD統(tǒng)計(jì)結(jié)果和保留小尺度電離層空間結(jié)構(gòu)的需求，最終選擇多視視數(shù)為32×16、濾波窗口為600像素的RSS電離層，電離層校正后的InSAR相位如圖1(d)所示。由圖1(d)可見，電離層校正后影像中仍存在部分小區(qū)域殘存相位，由于此處未完成對(duì)流層延遲校正，其結(jié)果仍存在對(duì)流層延遲相位或殘留地形誤差的影響，但原始InSAR相位中電離層基本得到去除，從而在一定程度上說明了RSS電離層校正的有效性。

由表1可見，隨著多視視數(shù)和濾波窗口的增加，電離層校正后InSAR相位標(biāo)準(zhǔn)差逐漸減小，由校正前的1.499rad減小至0.637rad，并在多視視數(shù)為32×16、濾波窗口為600像素時(shí)達(dá)到最小，這與基于IRI選取的RSS參數(shù)結(jié)果一致，初步驗(yàn)證了上述方法選取平滑參數(shù)的可行性。

表1 阿克里州中心經(jīng)線處RSS與IRI電離層差值的均方根偏差及電離層校正前后InSAR相位標(biāo)準(zhǔn)差

2.2 緩慢形變研究區(qū)

實(shí)驗(yàn)區(qū)位于甘肅省東南部的定西市北部安定區(qū)，主從影像成像時(shí)間為2018-05-22和2018-07-17。由于受到青藏高原東北部東昆侖斷裂的擠壓和西秦嶺北緣斷裂走滑等動(dòng)力因素的影響，該區(qū)域地震頻發(fā)[10]。主從影像成像時(shí)間段內(nèi)，距離影像西南部邊緣約50 km處的甘肅甘南州卓尼縣于2018-06-30發(fā)生3.1級(jí)地震，但該地震對(duì)研究區(qū)影響不大，因此可認(rèn)為該研究區(qū)為緩慢地表形變區(qū)。

采用與巴西阿克里州相同的實(shí)驗(yàn)流程及參數(shù)，得到結(jié)果如圖4和圖5所示。對(duì)比圖4(b)和圖5可見， RSS電離層條紋相比于IRI電離層條紋較為彎曲，但二者在整體趨勢(shì)上仍具有一定的相似性。由此推斷，緩慢地表形變區(qū)的電離層雖然受到一定擾動(dòng)因素的影響，但仍能將IRI模型作為參照。當(dāng)多視視數(shù)為16×8、濾波窗口達(dá)600像素及以上，以及多視視數(shù)為32×16、濾波窗口達(dá)400像素及以上時(shí)，增大濾波窗口對(duì)電離層條紋的改變不再明顯，因此需要在上述范圍中進(jìn)一步選擇平滑參數(shù)。

圖4 定西研究區(qū)處理結(jié)果Fig.4 The processing results of Dingxi city

圖5 定西研究區(qū)不同平滑參數(shù)下RSS法估算的電離層相位Fig.5 Ionospheric phase estimated by RSS with different smoothing parameters in Dingxi city

定西研究區(qū)中心經(jīng)線上RSS與IRI差值的RMSD計(jì)算結(jié)果如表2所示，總體上看，RMSD隨多視視數(shù)和濾波窗口的增加而減小。當(dāng)多視視數(shù)為16×8時(shí)，RMSD隨濾波窗口的增加緩慢減?。划?dāng)多視視數(shù)為32×16時(shí)，RMSD隨濾波窗口的增大迅速減小，但平滑參數(shù)選取過大時(shí)，小于濾波窗口的電離層結(jié)構(gòu)被平滑掉。為獲得較小的RMSD同時(shí)保留更多電離層空間結(jié)構(gòu)細(xì)節(jié)，在定西實(shí)驗(yàn)區(qū)仍選用32×16的多視視數(shù)和600像素尺寸的濾波窗口。

表2 定西地區(qū)中心經(jīng)線處RSS與IRI電離層差值的均方根偏差

由電離層校正后的InSAR結(jié)果(圖4(c))可以看出，在電離層校正后干涉相位中仍存在較為復(fù)雜的變化。不同于阿克里州實(shí)驗(yàn)區(qū)，定西地區(qū)在圖4(a)的InSAR相位條紋中表現(xiàn)出一些地形特征，因此在電離層校正后同樣殘存了對(duì)流層延遲和地形相位殘差，且存在一定的地表形變，這些因素導(dǎo)致校正后的InSAR相位仍存在一定的空間變化。

3 結(jié) 語

本文研究在InSAR電離層校正中以IRI電離層相位為參考選取合適的RSS平滑參數(shù)。通過目視對(duì)比可知，在無顯著地表形變的巴西阿克里州實(shí)驗(yàn)區(qū)和存在緩慢形變的中國(guó)甘肅省定西市實(shí)驗(yàn)區(qū)，IRI和RSS電離層相位在空間模式上均具有良好的相似性，從而進(jìn)一步驗(yàn)證了RSS方法的可行性。選用ALOS-2 PALSAR條帶模式數(shù)據(jù)計(jì)算RSS與IRI電離層相位偏差，確定一組多視視數(shù)為32×16(az×rng)、濾波窗口尺寸為600像素的平滑處理參數(shù)。根據(jù)此參數(shù)估算的電離層相位在兩個(gè)實(shí)驗(yàn)區(qū)均取得較好的校正結(jié)果，從而證明參考IRI模型選取RSS平滑參數(shù)在無形變區(qū)域和緩慢形變區(qū)的可行性。

RSS方法平滑參數(shù)的選取主要與數(shù)據(jù)帶寬和相干性有關(guān)，當(dāng)使用ALOS-2 PALSAR條帶模式數(shù)據(jù)在類似區(qū)域選取平滑參數(shù)時(shí)，可以簡(jiǎn)化處理環(huán)節(jié)、提高效率、取得更為可靠的結(jié)果，對(duì)于其他有顯著地表形變區(qū)域的RSS電離層平滑處理也具有一定的參考意義。

致謝:感謝日本宇航探測(cè)研究機(jī)構(gòu)JAXA提供PALSAR-2數(shù)據(jù)(PI No. 3321)及北京大學(xué)高性能計(jì)算校級(jí)公共平臺(tái)提供計(jì)算環(huán)境。