鄧 偉，宋迎春，李文娜，謝雪梅

(1.中南大學(xué) 地球科學(xué)與信息物理學(xué)院，長(zhǎng)沙 410083;2.中南林業(yè)科技大學(xué) 土木工程學(xué)院，長(zhǎng)沙 410004)

在大地測(cè)量數(shù)據(jù)處理中，隨著觀測(cè)手段的增多，觀測(cè)資料的累積，獲取先驗(yàn)信息的可能性也越來(lái)越大，通過(guò)獲取的先驗(yàn)信息建立約束條件一直是大地測(cè)量領(lǐng)域研究的熱點(diǎn)。傳統(tǒng)的最小二乘是基于殘差平方和最小的方法，但是其平差結(jié)果在一些特殊情況下并不能夠真實(shí)地表達(dá)出觀測(cè)對(duì)象的物理或者力學(xué)性質(zhì)。比如將不等式約束應(yīng)用到大壩變形監(jiān)測(cè)中，通過(guò)坡度方向和控制點(diǎn)的大致移動(dòng)方向建立不等式約束條件進(jìn)行解算，避免在無(wú)約束平差時(shí)兩期平差值可能朝向任意方向的問(wèn)題[1]。附加正確的、一定約束條件的平差模型，可以使平差結(jié)果在滿足最小二乘準(zhǔn)則的基礎(chǔ)上，得到符合觀測(cè)對(duì)象規(guī)律的更合適的解。在現(xiàn)代大地測(cè)量領(lǐng)域，運(yùn)用較為廣泛的就是等式約束[2-5]和不等式約束[1,6-8]，許多學(xué)者也對(duì)此進(jìn)行了大量的研究并提出許多可靠的方法。

但是，在實(shí)際工程應(yīng)用中，有許多用區(qū)間形式表達(dá)的參數(shù)約束信息。比如地球物理反演中的三維電阻率反演法，反演參數(shù)為電阻率，根據(jù)物理常識(shí)可以獲得電阻率范圍[9]。又如在大地測(cè)量參數(shù)反演中，反演參數(shù)巖土的密度值，由于密度是固有的正值，便能假設(shè)其大于零，且可以根據(jù)內(nèi)部巖石特征得到參數(shù)值必定在某個(gè)區(qū)間內(nèi)[10]。

針對(duì)參數(shù)帶區(qū)間約束的平差模型，許多學(xué)者進(jìn)行了各種算法研究。謝雪梅通過(guò)將平差問(wèn)題轉(zhuǎn)化為一個(gè)簡(jiǎn)單的二次規(guī)劃問(wèn)題，基于矩陣的正則分裂提出了一種迭代算法，并驗(yàn)證了該算法在系數(shù)矩陣病態(tài)情況下的有效性[11]；肖兆兵將區(qū)間約束轉(zhuǎn)化為二次規(guī)劃問(wèn)題，基于Kuhn-Tucker條件，針對(duì)參數(shù)的不確定度提出了一種平差迭代算法[12]。夏玉國(guó)基于積極集思想，將最優(yōu)化理論中的子空間截?cái)嗯ｎD法應(yīng)用到測(cè)量平差中，有效改善了區(qū)間約束下平差模型的計(jì)算效率[13]。

帶區(qū)間約束條件的最優(yōu)化問(wèn)題，有數(shù)學(xué)學(xué)者將區(qū)間約束條件轉(zhuǎn)化成橢球約束條件進(jìn)行求解[14]。即在區(qū)間約束周圍外接一個(gè)橢球，擴(kuò)大約束范圍，利用橢球約束方法進(jìn)行求解，將橢球特征矩陣的逆作為廣義嶺估計(jì)的嶺參數(shù)。這種方法雖然可以解決區(qū)間約束最小二乘問(wèn)題，遺憾的是解有時(shí)不在區(qū)間約束內(nèi)。本文借助數(shù)學(xué)上的這類方法，提出一種改進(jìn)的求解橢球約束的方法，通過(guò)對(duì)橢球特征矩陣進(jìn)行加權(quán)，減小特征矩陣的半徑來(lái)迭代求解，直至求出最優(yōu)解。

1 區(qū)間約束最小二乘平差問(wèn)題

帶有區(qū)間約束的平差模型為：

L=AX+e，權(quán)陣為P，

s.t.α≤X≤β.

(1)

式中：s.t.是subject to的縮寫，意思是“受約束”；A為m×n(m≥n)維系數(shù)矩陣;X=[x1,x2,…,xn]T為n維參數(shù)向量;L為m維線性化后觀測(cè)值與近似值之間的差值;e為n維隨機(jī)誤差向量;α=[α1,α2,…,αn]T為參數(shù)下界;β=[β1,β2,…,βn]T為參數(shù)上界。

L=AX+e,

(2)

式(2)中，根據(jù)最小二乘平差準(zhǔn)則VTPV=min，可以等價(jià)求解如下問(wèn)題：

min(L-AX)TP(L-AX),

(3)

構(gòu)造拉格朗日函數(shù)：

f(X)=(L-AX)TP(L-AX)+

(4)

分別對(duì)X和λ求導(dǎo)：

(5)

所以，

(6)

由于最優(yōu)解Xλ都是在橢球面上取得，問(wèn)題可以轉(zhuǎn)變?yōu)榍蠼庠跈E球面上的最小二乘問(wèn)題。用傳統(tǒng)的拉格朗日法求解式(6)的過(guò)程中，關(guān)鍵在于拉格朗日乘數(shù)λ的確定。文獻(xiàn)[14]中不再將λ看作一個(gè)常數(shù)，而是看作依賴于樣本的嶺參數(shù)。通過(guò)式λ=XTATPL-XTATPAX對(duì)其進(jìn)行迭代求解。文獻(xiàn)[15]基于平衡損失函數(shù)，提出一種既考慮估計(jì)的精度，又考慮模型擬合優(yōu)良程度的方法。

雖然將區(qū)間約束轉(zhuǎn)化成橢球約束可以減少約束方程，便于加入平差模型進(jìn)行求解，尤其是在待估參數(shù)向量維數(shù)比較高時(shí)，但是這是一種不對(duì)等的轉(zhuǎn)換。為了便于從幾何上來(lái)分析，以二維平面上的區(qū)間約束轉(zhuǎn)換成橢球約束作為參考。如圖1可以看出，將區(qū)間約束轉(zhuǎn)化成橢球約束相當(dāng)于在原本區(qū)間約束周圍外接一個(gè)橢球。如果此時(shí)再使用傳統(tǒng)的拉格朗日方法對(duì)橢球約束平差模型進(jìn)行求解，可能會(huì)導(dǎo)致參數(shù)解雖然在橢球面上，但是不符合原本的區(qū)間約束條件?；谶@個(gè)問(wèn)題，文中提出一種新的求解橢球約束平差模型的方法，即將橢球約束條件作為一個(gè)懲罰項(xiàng)加入平差模型中，構(gòu)造出一個(gè)新的輔助函數(shù)。懲罰項(xiàng)不斷加權(quán)，按照一定步長(zhǎng)減小約束橢球的半徑，直至確定滿足原始區(qū)間約束條件的解為止。

圖1 區(qū)間約束區(qū)域與橢球約束區(qū)域比較

2 二次罰函數(shù)法

二次罰函數(shù)法是罰函數(shù)法的一種特殊形式[16]，懲罰項(xiàng)為約束條件的平方，優(yōu)點(diǎn)在于比較簡(jiǎn)單、直觀。對(duì)于等式約束最小二乘問(wèn)題：

minf(x) s.t.ci(x)=0.

(7)

其中f(x)為目標(biāo)函數(shù)；ci(x)=0為等式約束條件，式(7)的二次懲罰函數(shù)形式為：

(8)

其中，ρ>0為懲罰因子。二次懲罰法將違反每個(gè)約束的平方的倍數(shù)加到目標(biāo)函數(shù)上，在ρ→∞的過(guò)程中，對(duì)約束條件的“懲罰”就越大，可以求解出一系列的無(wú)約束解。

構(gòu)造罰函數(shù)：

(9)

式中：μ為加權(quán)因子;Gμ為權(quán)矩陣，Gμ=μG，0≤μ≤1。通過(guò)式(9)可以看出，在μ變化的過(guò)程中，每一個(gè)橢球面上都可以求出一個(gè)在最小二乘準(zhǔn)則下的解，但是要對(duì)該解進(jìn)行判斷，判斷其是否滿足區(qū)間約束條件。μ從1開(kāi)始減小，直到確定不等式(3)的解為止。

F(X)=

(10)

求F(X)關(guān)于X的一階導(dǎo)數(shù)并令其為0。

(11)

令二次罰函數(shù)法求得的解為：

(12)

二次罰函數(shù)法具體步驟：

1)利用先驗(yàn)信息確定參數(shù)上下界α和β。

3)系數(shù)矩陣A，觀測(cè)向量L，權(quán)矩陣P，初始懲罰項(xiàng)Gμ，μ=1作為輸入項(xiàng)。

5)若對(duì)于Xqp中的所有元素，都滿足α≤Xqp≤β(即參數(shù)值在約束區(qū)間之內(nèi))，則算法終止。

3 二次罰函數(shù)正則化

典型的正則化方法是Tikhonov正則化[17]，標(biāo)準(zhǔn)Tikhonov正則化的形式為：

(13)

式中：α為正則化參數(shù)。確定正則化參數(shù)的方法主要有L-曲線法[18]、廣義交叉核實(shí)法(Generalized Cross-Validation，GCV)[19]和嶺跡法[20]。鑒于廣義交叉核實(shí)法適用性較弱和嶺跡法選取的主觀性太強(qiáng)的缺點(diǎn)。文中選擇使用適用性強(qiáng)，并且能獲得較為合理的正則化參數(shù)的L-曲線法來(lái)選取正則化參數(shù)。

在二次罰函數(shù)法中引入正則化項(xiàng)，即：

minf(X)=

(14)

求f(X)關(guān)于X的一階導(dǎo)數(shù)并令其為0。

(15)

(16)

對(duì)于參數(shù)帶區(qū)間約束的平差模型，由于其雙邊不等式約束的特殊形式，直接求解可能比較麻煩。在計(jì)算時(shí)通常將其轉(zhuǎn)化成橢球約束進(jìn)行求解，變成在橢球面上求解最小二乘解的問(wèn)題。針對(duì)橢球約束會(huì)弱化不等式約束，導(dǎo)致解可能不在區(qū)間約束內(nèi)的缺點(diǎn)，文中提出的二次罰函數(shù)法可以很好解決這個(gè)問(wèn)題。

當(dāng)系數(shù)矩陣病態(tài)時(shí)，雖然橢球約束會(huì)對(duì)模型有一定的改善，但是這種改善太依賴于參數(shù)的先驗(yàn)信息?；谶@樣的缺點(diǎn)，文中提出在二次罰函數(shù)的基礎(chǔ)上加入正則化項(xiàng)。既能通過(guò)約束條件使得最小二乘解滿足觀測(cè)對(duì)象的實(shí)際情況，又能更大限度解決系數(shù)矩陣病態(tài)帶來(lái)的解不穩(wěn)定的問(wèn)題。

橢球約束正則化最小二乘法步驟：

1)利用先驗(yàn)信息確定參數(shù)上下界α和β。

3)將系數(shù)矩陣A，觀測(cè)向量L，權(quán)矩陣P，初始懲罰項(xiàng)Gμ作為輸入項(xiàng)。

4)根據(jù)L-曲線法求出正則化參數(shù)α。

6)若對(duì)于Xqpr中的所有元素，都滿足α≤Xqpr≤β(解算出的參數(shù)值在約束區(qū)間之內(nèi))，則算法終止。

4 算例與分析

在地球物理領(lǐng)域常常會(huì)遇到不適定問(wèn)題，例如三維電阻率反演。其正演過(guò)程是基于有限差分或有限元一類計(jì)算方法的數(shù)值解，反演則是求數(shù)量巨大的各網(wǎng)格剖分單元電阻率值。反演需要求取的參數(shù)是各網(wǎng)格單元的電阻率值，而電阻率值是可以通過(guò)物理常識(shí)確定區(qū)間范圍的，可以將約束范圍條件加入平差模型進(jìn)行解算[9]。

將非線性的三維電阻率反演問(wèn)題線性化后，可以得到如下的反演方程：

AΔm=Δd.

式中，A為敏感度矩陣;Δm為模型參數(shù)增量向量;Δd為觀測(cè)數(shù)據(jù)dobs與正演理論值dm的殘差向量;dm是根據(jù)給定的模型參數(shù)由數(shù)值正演得到的理論觀測(cè)值。三維電阻率反演問(wèn)題往往表現(xiàn)為混定問(wèn)題，導(dǎo)致方程為病態(tài)方程。

通常在模型中加入如下的不等式約束：

ρli≤mi≤ρui.

其中，mi為第i個(gè)網(wǎng)格的電阻率;ρli和ρui分別為第i個(gè)網(wǎng)格電阻率的下限和上限。對(duì)于ρli和ρui，可以根據(jù)物理常識(shí)來(lái)得到其變化的大致范圍。

綜合反演方程和不等式約束模型，可以得到如下的三維電阻率反演的目標(biāo)函數(shù)：

F=min(Δd-AΔm)T(Δd-AΔm),

ρli≤mi≤ρui.

為了評(píng)價(jià)文中算法在大地測(cè)量反演中的效果，設(shè)計(jì)了如下算例進(jìn)行驗(yàn)證。算例中用一個(gè)希爾伯特矩陣代替三維電阻率反演中線性化后的系數(shù)矩陣，真值假設(shè)為Xzhen=[1,2,3,4,5,6]T，假設(shè)根據(jù)常識(shí)得到的電阻率反演范圍為：l=[0 0 0 0 0 0]T，u=[5 6 7 8 9 10]T，觀測(cè)值L為根據(jù)系數(shù)矩陣和真值計(jì)算出來(lái)的觀測(cè)值真值加上隨機(jī)誤差表示。A和L分別為：

L=AXzhen+e.

e為服從正態(tài)分布的隨機(jī)數(shù)組成的向量，用其來(lái)模擬測(cè)繪實(shí)際工程中的觀測(cè)誤差。為了避免數(shù)據(jù)的特殊性，采用50次數(shù)值模擬實(shí)驗(yàn)，分別采取以下5種方案進(jìn)行計(jì)算：

1)最小二乘法。

2)嶺估計(jì)法(嶺參數(shù)采用L-曲線法確定)。

3)奇異值分解法。

4)傳統(tǒng)橢球約束法。

5)二次罰函數(shù)正則化法。

令最小二乘解為Xls，嶺估計(jì)解Xre，截?cái)嗥娈愔捣纸獾慕鉃閄svd，傳統(tǒng)橢球約束的解為Xec，本文算法的解為Xqpr。

m1=(Xls-Xzhen)T(Xls-Xzhen),

m2=(Xre-Xzhen)T(Xre-Xzhen),

m3=(Xsvd-Xzhen)T(Xsvd-Xzhen),

m4=(Xec-Xzhen)T(Xec-Xzhen),

m5=(Xqpr-Xzhen)T(Xqpr-Xzhen).

m1，m2，m3，m4，m5分別代表每種方法求解出的參數(shù)偏離模擬真值的程度，結(jié)果如表1所示。

圖2 最小二乘法計(jì)算結(jié)果

圖3 嶺估計(jì)法計(jì)算結(jié)果

圖4 截?cái)嗥娈愔捣纸夥ㄓ?jì)算結(jié)果

圖5 傳統(tǒng)橢球約束計(jì)算結(jié)果

圖6 文中算法計(jì)算結(jié)果

對(duì)模擬算例結(jié)果進(jìn)行分析：

1)從表1可以看出，解偏離模擬真值的程度從大到小依次為：Xls>Xre>Xsvd>Xec>Xqpr。

2)設(shè)計(jì)矩陣A的條件數(shù)為1.495 1×108，屬于嚴(yán)重病態(tài)矩陣。從圖2可以看出，最小二乘解極不穩(wěn)定，此時(shí)最小二乘解已經(jīng)嚴(yán)重失真。

3) 嶺估計(jì)法對(duì)于參數(shù)的估計(jì)精度太依賴于嶺參數(shù)的選取，對(duì)模型的病態(tài)程度改善不夠理想。解的穩(wěn)定性相較于最小二乘法雖然有所改善，但是其只注重于降低病態(tài)性，沒(méi)有考慮被觀測(cè)對(duì)象的實(shí)際情況。

4) 截?cái)嗥娈愔捣纸夥ㄊ峭ㄟ^(guò)去除設(shè)計(jì)矩陣接近于0的奇異值，通過(guò)損失待求參數(shù)的無(wú)偏性為代價(jià)來(lái)?yè)Q取均方根誤差的降低。解的偏離真值程度相較于最小二乘雖有所改善，但是依然很大，且解不穩(wěn)定。

5) 利用參數(shù)先驗(yàn)信息，將其作為約束條件加入平差模型，可以極大提高解的穩(wěn)定性。但是傳統(tǒng)的橢球約束算法會(huì)出現(xiàn)解不在約束區(qū)間里面的情況，從圖5可以看出，部分x1出現(xiàn)了小于0的情況。

6) 利用文中算法求出的解穩(wěn)定性較好，不僅可以解決傳統(tǒng)橢球約束算法的解可能不在區(qū)間內(nèi)的缺點(diǎn)，而且通過(guò)加入正則化項(xiàng)，可以求得正則化參數(shù)與橢球特征矩陣逆的最優(yōu)組合，更好地降低模型的病態(tài)性。

5 結(jié)束語(yǔ)

在大地測(cè)量數(shù)據(jù)處理領(lǐng)域，通過(guò)在平差模型中加入有效的先驗(yàn)信息作為約束條件，可以降低模型的不穩(wěn)定性，提高解的精度和可靠性。文中在傳統(tǒng)的橢球約束解算方法上，通過(guò)引入數(shù)學(xué)領(lǐng)域的二次罰函數(shù)法，對(duì)橢球特征矩陣進(jìn)行加權(quán)，解決了傳統(tǒng)橢球約束算法的解可能不在約束區(qū)間內(nèi)的問(wèn)題。傳統(tǒng)橢球約束是以橢球特征矩陣的逆為嶺參數(shù)的廣義嶺估計(jì)，文中在二次罰函數(shù)的目標(biāo)函數(shù)中引入正則化項(xiàng)，選出正則化參數(shù)與橢球特征矩陣逆的最優(yōu)組合。算例表明，文中算法在求解區(qū)間約束病態(tài)模型上有優(yōu)勢(shì)。

表1 解偏離真值的程度