翟心儀 潘光永 陳 洋 劉國(guó)林

1 山東科技大學(xué)測(cè)繪與空間信息學(xué)院，青島市前灣港路579號(hào)，266590

短基線集干涉合成孔徑雷達(dá)(SBAS-InSAR)技術(shù)在地震、火山、滑坡、泥石流、城市沉降、河口沉降等災(zāi)害監(jiān)測(cè)中有廣泛應(yīng)用[1-4]。然而在實(shí)際使用中發(fā)現(xiàn)，SBAS-InSAR形變模型病態(tài)的可能性始終存在，這會(huì)直接影響形變反演的精度和可靠性[5-6]。

譜修正迭代法優(yōu)良的性質(zhì)及其未改變方程等量關(guān)系的特點(diǎn)，使其在病態(tài)問題的解算中得到廣泛關(guān)注[7-12]。通過譜修正迭代的收斂性證明不難看出，其充分迭代結(jié)果收斂于最小二乘估計(jì)結(jié)果，現(xiàn)有的迭代終止條件一般是根據(jù)相鄰兩次迭代估值是否小于閾值進(jìn)行判定，但這意味著不同的閾值會(huì)導(dǎo)致譜修正迭代估計(jì)結(jié)果發(fā)生變化，這種經(jīng)驗(yàn)性的閾值選取會(huì)帶來(lái)較大的不確定性。因此，本文針對(duì)譜修正迭代終止條件的優(yōu)化方法展開研究，提出一種L曲線譜修正迭代法實(shí)現(xiàn)迭代結(jié)果的優(yōu)選，并通過SBAS-InSAR形變反演實(shí)驗(yàn)對(duì)方法的有效性進(jìn)行驗(yàn)證。

1 SBAS-InSAR形變反演

假設(shè)實(shí)驗(yàn)區(qū)域內(nèi)存在Ns+1幅時(shí)間序列為t0,t1,…,tNs的單視復(fù)數(shù)(singal look complex, SLC)影像(圖1)，根據(jù)時(shí)空基線條件得到M個(gè)差分干涉對(duì)。若第j幅干涉圖由tA和tB(tA>tB)時(shí)刻的影像得到，去除地形相位可得到差分干涉圖中(r,x)處的差分干涉相位為：

δφj(r,x)=φ(tB,r,x)-φ(tA,r,x)

(1)

φatm+φnon-lin+φnoise+φtop-error

(2)

式中，d(tA,r,x)和d(tB,r,x)為L(zhǎng)OS向累積形變量，φ(tA,r,x)和φ(tB,r,x)為形變相位值，λ為雷達(dá)波長(zhǎng)，φatm為大氣延遲相位，φnon-lin為非線性形變相位，φnoise為噪聲相位，φtop-error為地形殘余相位。

圖1 SBAS-InSAR時(shí)序差分干涉對(duì)的組成方式Fig.1 Combination of time series differential interference pairs of SBAS-InSAR

不同時(shí)刻(t0,t1,…,tNs)相對(duì)于t0時(shí)刻的形變相位向量φT可表示為：

φT=[φ(t1),…,φ(tNs)]

(3)

差分干涉圖解纏后的相位觀測(cè)值δφT可表示為：

δφT=[δφ1,…,δφNs]

(4)

通過組合差分干涉圖可得：

(5)

(6)

不難看出，時(shí)序差分干涉對(duì)的質(zhì)量對(duì)形變反演精度有著重要影響。由于產(chǎn)生時(shí)序差分干涉對(duì)的方法仍以經(jīng)驗(yàn)判斷為主，因此形變反演模型中可能存在病態(tài)性，從而導(dǎo)致最小二乘估計(jì)結(jié)果難以滿足精度要求。

2 譜修正迭代法

(7)

法方程為：

(8)

式中，N=BTPB，W=BTPL，L為觀測(cè)向量，P為權(quán)陣。

(9)

式中，I為單位陣。則譜修正迭代法的公式為：

(10)

L=AY+Δ

(11)

(12)

(13)

(14)

由于譜修正迭代法最終收斂于最小二乘估計(jì)[7-9]，當(dāng)譜修正迭代法的迭代初值給定后，可得：

(15)

(16)

則譜修正迭代法最終估值的均方誤差可表示為：

(17)

有偏估計(jì)方法通常是引入正則化參數(shù)或正則化矩陣，通過降低解算方程的病態(tài)程度，從而降低病態(tài)性對(duì)估計(jì)結(jié)果的影響程度。而譜修正迭代法在不改變方程等量關(guān)系的情況下，通過合理確定迭代估計(jì)結(jié)果改善病態(tài)問題估計(jì)的結(jié)果。

3 L曲線譜修正迭代

譜修正迭代法的迭代終止條件直接影響估計(jì)結(jié)果的精度。通過設(shè)置閾值確定譜修正迭代法終止條件，易使得估計(jì)結(jié)果的穩(wěn)定性變差，估值精度存在不確定性。因此，有必要改進(jìn)譜修正迭代法的迭代終止條件。

(18)

為了實(shí)現(xiàn)譜修正迭代結(jié)果的優(yōu)選，提出一種L曲線譜修正迭代法，利用L曲線的最大曲率點(diǎn)，確定殘差較小且穩(wěn)定的迭代結(jié)果作為最終估值。構(gòu)造譜修正迭代法的估計(jì)準(zhǔn)則為：

(19)

(20)

(21)

(22)

選擇最大曲率點(diǎn)作為終止迭代的條件，迭代次數(shù)m對(duì)應(yīng)的估計(jì)值即為估值穩(wěn)定且殘差較小的譜修正迭代法的最終估值。

圖2 基于L曲線的譜修正迭代法Fig.2 L-curve characteristic value correction iteration method

4 SBAS-InSAR濟(jì)寧礦區(qū)沉降實(shí)驗(yàn)

濟(jì)寧市擁有豐富的煤炭資源，但隨著開采量的增加，礦區(qū)沉降問題日益嚴(yán)重[13]。采用2008-01～2010-09期間覆蓋實(shí)驗(yàn)區(qū)的13景Envisat ASAR影像(極化方式VV，波段C，波束入射角23°，方位向分辨率4 m，距離向分辨率7 m)，采用多主影像自由組合進(jìn)行干涉配對(duì)，形成46個(gè)干涉對(duì)(圖3)，選取58 751個(gè)高相干點(diǎn)，構(gòu)建形變模型系數(shù)矩陣的條件數(shù)為519，呈病態(tài)性。研究區(qū)域內(nèi)既有沉降，也存在抬升，最大年平均沉降速率大于-40 mm/a。

圖3 干涉對(duì)組合方式Fig.3 Combination of interference pairs

采用最小二乘估計(jì)(LS)、譜修正迭代法(PXZ，閾值為0.01、0.001)、L曲線譜修正迭代法(PXZ-L)進(jìn)行解算，圖4給出了不同估計(jì)方法的年平均沉降速率。表1為L(zhǎng)S、PXZ-0.01、PXZ-0.001、PXZ-L解算的高相干點(diǎn)的均方根誤差統(tǒng)計(jì)結(jié)果。將實(shí)驗(yàn)區(qū)內(nèi)5個(gè)點(diǎn)的水準(zhǔn)數(shù)據(jù)作為精度評(píng)定的依據(jù)，表2和3分別列出4種方法相對(duì)于水準(zhǔn)結(jié)果的沉降量和沉降速率中誤差。所選5個(gè)水準(zhǔn)點(diǎn)位于沉降區(qū)域的邊緣，最大年平均沉降速率均小于-6 mm/a。

表1 均方根誤差統(tǒng)計(jì)

由圖4可知，最小二乘估計(jì)的年平均沉降速率范圍為-43.91～14.94 mm/a；譜修正迭代法(閾值分別為0.01、0.001)的年平均沉降速率范圍為-43.70～14.70mm/a、-43.88～14.87mm/a；L曲線譜修正迭代法的年平均沉降速率范圍為-48.90～14.89 mm/a?？梢姡懿B(tài)影響，各種方法的估計(jì)結(jié)果存在差異。

表2 不同估計(jì)方法解算的沉降量與水準(zhǔn)點(diǎn)監(jiān)測(cè)結(jié)果的差值

表3 不同估計(jì)方法解算的年平均沉降速率與水準(zhǔn)點(diǎn)監(jiān)測(cè)結(jié)果的差值

圖4 不同估計(jì)方法解算得到的年平均沉降速率Fig.4 Annual average subsidence rate calculated by different estimation methods

由表1可以看出，形變反演結(jié)果的均方根誤差集中分布于0～15 mm之間，最小二乘估計(jì)、譜修正迭代法(閾值為0.01、0.001)、L曲線譜修正迭代法的較小均方根誤差頻數(shù)占比分別為97.82%、97.85%、97.82%和97.88%。因此，不同閾值條件下的譜修正迭代法的估計(jì)結(jié)果存在差異，使得對(duì)病態(tài)問題的改善程度存在較大不確定性；而L曲線譜修正迭代法通過選擇估值穩(wěn)定且殘差較小的估計(jì)結(jié)果，實(shí)現(xiàn)了估計(jì)結(jié)果的優(yōu)選，其較小均方根誤差頻數(shù)占比也最高。

比較表2中水準(zhǔn)監(jiān)測(cè)結(jié)果與各估計(jì)方法反演得到的最鄰近高相干點(diǎn)的沉降量差值可以看出，譜修正迭代法(閾值為0.01、0.001)以及L曲線譜修正迭代法反演后的沉降量結(jié)果依次為σPXZ-0.01=2.908 mm，σPXZ-0.001=2.977 mm及σPXZ-L=2.867 mm，小于最小二乘估計(jì)方法反演得到的沉降量結(jié)果σLS=3.035 mm。表3中各方法估計(jì)的年平均沉降速率差值結(jié)果顯示，3種方法估計(jì)的沉降速率結(jié)果依次為σPXZ-0.01=1.077 mm/a，σPXZ-0.001=1.103 mm/a及σPXZ-L=1.062 mm/a，同樣小于最小二乘方法估計(jì)的沉降速率結(jié)果σLS=1.124 mm/a。沉降量與沉降速率的結(jié)果均表明，譜修正迭代法可減弱病態(tài)性對(duì)估計(jì)結(jié)果的不良影響，其中，基于L曲線的譜修正迭代法的差值最小，證明L曲線譜修正迭代法是可行的。

5 結(jié) 語(yǔ)

譜修正迭代法通過在法方程系數(shù)矩陣的兩端同時(shí)增加未知參數(shù)，在不改變方程等量關(guān)系的情況下，實(shí)現(xiàn)對(duì)系數(shù)矩陣病態(tài)性的改善，在整個(gè)計(jì)算過程中不需要額外計(jì)算正則化參數(shù)或正則化矩陣，因此在病態(tài)問題的解算中具有獨(dú)特的優(yōu)勢(shì)。但譜修正迭代終止條件通常是依據(jù)相鄰迭代結(jié)果小于閾值，而閾值的選取越小，迭代結(jié)果越接近于最小二乘估計(jì)，因此基于閾值的迭代終止條件無(wú)益于估計(jì)結(jié)果的優(yōu)選。為此，本文提出一種利用L曲線進(jìn)行譜修正迭代估值優(yōu)選的方法，并通過SBAS-InSAR形變模型反演實(shí)驗(yàn)進(jìn)行驗(yàn)證分析。結(jié)果表明，基于L曲線的譜修正迭代法能夠提高解的精度和可靠性。