周巖，靳奉祥，梁慶華，馬德鵬

(1.山東科技大學(xué)資源與土木工程系，山東泰安，271019；2.山東建筑大學(xué)測繪地理信息學(xué)院，山東濟(jì)南，250101；3.中煤科工集團(tuán)重慶研究院有限公司，重慶，400039；4.山東農(nóng)業(yè)大學(xué)水利土木工程學(xué)院，山東泰安，271018)

半?yún)?shù)模型是20世紀(jì)80年代初發(fā)展起來的一種回歸模型，它是將系統(tǒng)誤差或者建模近似引起的模型偏差視為非隨機(jī)參數(shù)，采用補(bǔ)償最小二乘估計(jì)法，得到參數(shù)和非參數(shù)的估值[1-3]，因此，半?yún)?shù)模型廣泛地應(yīng)用于模型精化和減弱系統(tǒng)誤差方面[4-5]。半?yún)?shù)模型補(bǔ)償最小二乘估計(jì)中，關(guān)鍵的因素在于如何確定合適的正則化參數(shù)α及正則化矩陣R[6]。在測量平差中，正則化矩陣R可根據(jù)實(shí)際情況按自然樣條光滑法、時(shí)間序列的特性、先驗(yàn)方差的特性以及觀測量之間的某種距離來確定，而正則化參數(shù)α的確定方法目前主要有GCV法、L曲線法[7-13]等。GCV法在理論上能夠獲取最優(yōu)的正則化參數(shù)，但有時(shí)GCV 函數(shù)的變化平緩，定位其最小值有困難。HANSEN 等[14-15]用L曲線法確定病態(tài)問題中的嶺估計(jì)參數(shù)，并詳細(xì)介紹了L曲線法的基本思想和相關(guān)性質(zhì)；FISCHER等[16-17]在用半?yún)?shù)回歸模型處理大地測量數(shù)據(jù)時(shí)，將L曲線法引入半?yún)?shù)模型中確定正則化參數(shù)，并進(jìn)行了相關(guān)證明。現(xiàn)有文獻(xiàn)表明，L曲線法可以比較容易地獲得正則化參數(shù)，并且能夠得到比較精確的解，但L曲線法在求解過程可能不收斂。U曲線法是根據(jù)定義的U(α)函數(shù)，即以構(gòu)成L曲線法的基于正則化參數(shù)α>0的信號范數(shù)與殘差范數(shù)為基礎(chǔ)，由信號范數(shù)的倒數(shù)與殘差范數(shù)的倒數(shù)之和組成，以U(α)為縱軸，α為橫軸所表示的曲線[18]，其形態(tài)接近U 型，共包括3 個(gè)組成部分：左邊和右邊部分幾乎是“垂直”的，中間部分幾乎是“水平”的，曲線上到原點(diǎn)距離最小的點(diǎn)對應(yīng)的α或者曲線左側(cè)曲率最大點(diǎn)對應(yīng)的α即為所確定的正則化參數(shù)。通過上述分析，本文作者在簡要介紹L曲線法的基礎(chǔ)上，詳細(xì)分析U曲線法，并闡述其相關(guān)特性，同時(shí)，分別采用L曲線法和U曲線法對模擬數(shù)值算例進(jìn)行求解，通過參數(shù)、非參數(shù)估計(jì)解及單位權(quán)中誤差的相互對比來分析U曲線法和L曲線法的異同。

1 確定正則化參數(shù)的L曲線法

由半?yún)?shù)回歸模型L=BX+S+Δ，可得誤差方程式為

式中：K為拉格朗日常數(shù)；α為正則化參數(shù)；R為正則化矩陣。

得

這里正則化矩陣R正定或者半正定，則可求得參數(shù)、非參數(shù)估計(jì)解及單位權(quán)中誤差，如下式所示[20]：

在半?yún)?shù)模型補(bǔ)償最小二乘估計(jì)中，確定正則化參數(shù)和正則化矩陣尤為重要，本文在選取正則化矩陣R=GTG的基礎(chǔ)上，重點(diǎn)研究不同的正則化參數(shù)確定方法對參數(shù)估值的影響。

L 曲線法是基于正則化參數(shù)α的殘差范數(shù)和信號范數(shù)所表示的曲線，二者都是正則化參數(shù)α的函數(shù)，選擇不同的α，便得到不同的點(diǎn)(Vn(α)，Sn(α))，較為直觀的形式就是以Sn(α)為橫軸，以Vn(α)為縱軸作圖，經(jīng)擬合得到1條曲線。該曲線包括2個(gè)典型部分，分別是“水平”和“垂直”部分，越“水平”部分對應(yīng)的是正則化參數(shù)越小，主要由Sn(α)控制；越“垂直”的部分對應(yīng)的正則化參數(shù)越大，主要由Vn(α)控制；將所得到的類似字母L的曲線稱為L曲線。從圖1觀察信號范數(shù)和噪音范數(shù)之間是否取得了平衡，或兩者是否出現(xiàn)了過大或過小的情況，因此，L曲線法顧及了擬合和光滑之間的平衡，拐點(diǎn)處的值所對應(yīng)的“角”即為所求的正則化參數(shù)的最優(yōu)解。選取最后的正則化參數(shù)值有2種標(biāo)準(zhǔn)：一是選取離原點(diǎn)最近的拐點(diǎn)；二是選取曲率k最大的點(diǎn)。

圖1 典型的L曲線示意圖Fig.1 Schematic diagram of typical L curve

2 確定正則化參數(shù)的U曲線法

定義

式中：α>0；x(α)和y(α)分別與L 曲線法中Vn(α)和Sn(α)的定義相同。

從圖2可見U曲線包含3個(gè)典型部分：

1)U曲線在左邊部分和右邊部分近似垂直。

2)U曲線中間部分近似水平。

3)垂直部分對應(yīng)的正則化參數(shù)是代表了Vn(α)和Sn(α)分別各自主導(dǎo)部分。水平部分對應(yīng)的正則化參數(shù)表征Vn(α)和Sn(α)較接近。

由于x(α)和y(α)分別等同于L 曲線里定義的Vn(α)和Sn(α)，故式(6)可以表示為

對式(7)求導(dǎo)得

圖2 典型的U曲線示意圖Fig.2 Schematic diagram of typical U curve

1)在L曲線法中，隨著α增大，Vn(α)是單調(diào)遞增函數(shù)，Sn(α)是單調(diào)遞減函數(shù)，即在式(8)中，-Sn'(α)>0。由于Vn(α)是單調(diào)遞增函數(shù)，Sn(α)是單調(diào)遞減函數(shù)，根據(jù)式(7)和(8)可知，當(dāng)α較小時(shí)，占控制部分的是αSn2(α)，即Vn2(α)-αSn2(α)<0；當(dāng)α較大時(shí)，占控制部分的是Vn2(α)，即Vn2(α)-αSn2(α)>0，所以，U(α)是先減小后增大。

當(dāng)Vn(α)和Sn(α)這2 部分平衡時(shí)，通過U(α)曲線可知U(α)接近水平部分，因此，曲線呈U狀，兩端近似豎直，中間近似水平。

2)正則化參數(shù)的選擇。

方法a：

式中：αD為曲線上的點(diǎn)到原點(diǎn)的距離最小時(shí)對應(yīng)的α。

方法b：

式中：αK為曲線左側(cè)曲率最大時(shí)對應(yīng)的α。

3 算例分析

3.1 模擬算例1

設(shè)已知P1～P4這4個(gè)觀測點(diǎn)上的重力異常觀測值L及其坐標(biāo)，如表1所示。觀測誤差方差DΔ=(0.03)2Ι，試估計(jì)已測點(diǎn)重力異常Δg[19]。

物理大地測量中很早就用最小二乘配置法得到重力異常最佳估計(jì)值，下面用半?yún)?shù)模型的方法來研究重力異常。在許多地區(qū)，重力異常不僅包含隨機(jī)部分S，而且包含系統(tǒng)部分BX，故重力異常觀測方程可寫為半?yún)?shù)模型。已有的采用半?yún)?shù)模型求解的文獻(xiàn)中，把正則化參數(shù)α取為1[20]，也有的采用L 曲線法對正則化參數(shù)取值。在這里，采用本文所提出的U曲線法選取正則化參數(shù)，并同L曲線法進(jìn)行比較，同時(shí)注意到在重力異常的求解中協(xié)方差函數(shù)是距離的某種函數(shù)，因此，本文中的正則化矩陣采用簡單的距離函數(shù)，即

表1 重力異常觀測值及其位置Table 1 Observations of gravity abnormity and its locations

函數(shù)模型中，重力異常包含隨機(jī)異常和系統(tǒng)異常2部分。系統(tǒng)部分一般表示為坐標(biāo)(x,y)的線性函數(shù)：

按照半?yún)?shù)模型解算：

L曲線法與U曲線法確定正則化參數(shù)的示意圖如圖3和圖4所示，2 種方案計(jì)算的重力異常估值差異的比較見表2。

圖3所示為縱軸Vn(α)與橫軸Sn(α)的關(guān)系圖，根據(jù)其對應(yīng)的方法得到最佳α的取值，α=1.85。圖4所示為定義的U(α)函數(shù)(縱軸)與正則化參數(shù)α(橫軸)的關(guān)系圖。從圖3和圖4可以看出U 曲線法比L 曲線法在得到最佳α的取值上計(jì)算量要小，而且更直觀。

圖3 L曲線法確定的正則化參數(shù)αFig.3 Determining regularization parameterα according to L-curve method

圖4 U曲線法確定的正則化參數(shù)αFig.4 Determining regularization parameterα according to U-curve method

由表2可以看出：U曲線法確定的正則化參數(shù)使參數(shù)估計(jì)解的精度要比L曲線法的高，U曲線法計(jì)算的重力異常估計(jì)結(jié)果較接近于文獻(xiàn)[19]中的結(jié)果，這說明在此算例中，U曲線法求得的正則化參數(shù)的可靠性要比L曲線法的高。

3.2 模擬算例2（周期性系統(tǒng)誤差）

設(shè)L=BX+S(ti)+Δ，B=(Bi,j)100×2，Bi,1=ti，Bi,2=(ti)2，S(ti)=3sin(ti)sin(3ti)，ti=2π(i-1)/100，X=[1 0.25]T，Δ～(0,1)[20]；i=1,2,…,100。

表2 2種方案計(jì)算的重力異常估值差異的比較Table 2 Comparison of value difference of gravity anomaly between two calculate solutions

該算例模擬的是周期性系統(tǒng)誤差，采用基于正則化矩陣的補(bǔ)償最小二乘法求解，R矩陣由相鄰兩觀測點(diǎn)模型誤差之差的平方和形式確定，設(shè)計(jì)以下2種方案：方案1，采用L曲線法確定正則化參數(shù)；方案2，采用U曲線法確定正則化參數(shù)。

這2種方案計(jì)算的模擬系統(tǒng)誤差的估值與殘差估值分別見圖5和圖6(模擬的殘差為模擬參數(shù)的真值與模擬系統(tǒng)誤差的值代入誤差方程得到，半?yún)?shù)估計(jì)殘差為模擬參數(shù)的估值與模擬系統(tǒng)誤差的估值代入誤差方程得到，下同)。計(jì)算所得的X^ 估值與模擬X真值差異的比較如表3所示。

從圖5(a)和圖6(a)可以看出：當(dāng)模擬的系統(tǒng)誤差呈周期性時(shí)，L曲線法和U曲線法確定的正則化參數(shù)都使得模擬的系統(tǒng)誤差與模擬系統(tǒng)誤差的估值變化的趨勢基本一致，但U曲線法確定的正則化參數(shù)解算得到的系統(tǒng)誤差的估值與其真值的吻合程度要比L曲線法的高。從圖5(b)和圖6(b)可以看出：L 曲線法和U曲線法確定的正則化參數(shù)同樣也使得半?yún)?shù)估計(jì)殘差與模擬殘差變化的趨勢基本一致，并且半?yún)?shù)模型解算觀測值的殘差經(jīng)t檢驗(yàn)，也都以大于95%的概率服從正態(tài)分布。這是由于半?yún)?shù)模型能夠?qū)⒂^測值中的系統(tǒng)誤差分離出來，減弱其對觀測值的影響，從而得到較可靠的參數(shù)估值；同時(shí)，U曲線法確定的正則化參數(shù)解算得到的半?yún)?shù)估計(jì)殘差與模擬殘差的吻合程度也比L曲線法的高。

圖5 L曲線法確定正則化參數(shù)α(周期性系統(tǒng)誤差)Fig.5 Determination of regularization parameterα according to L-curve method(system error of periodicity)

圖6 U曲線法確定正則化參數(shù)α(周期性系統(tǒng)誤差)Fig.6 Determination of regularization parameterα according to U-curve method(system error of periodicity)

從表3可知：L 曲線法和U 曲線法確定的正則化參數(shù)差別不大，都能使所求結(jié)果與真值較接近，均提高了半?yún)?shù)模型解的精度。但通過進(jìn)一步對比分析發(fā)現(xiàn)：采用U 曲線法得參數(shù)估值為[1.011 3 0.235 1]T，與其模擬X真值差異=3.497 0×10-4；采用L 曲線法得到的參數(shù)估值為[1.018 0 0238 2]T，與其模擬X真值差異=4.632 4×10-4。與后者相比，U 曲線法確定的正則化參數(shù)求得參數(shù)估值更接近于真值，表明U 曲線法對正則化參數(shù)的確定具有一定的有效性，所得結(jié)果具有可信度。

3.3 模擬算例3（線性、周期性系統(tǒng)誤差）

為了更接近于實(shí)際情況，該算例模擬的是線性、周期性系統(tǒng)誤差，采用基于正則矩陣的補(bǔ)償最小二乘法求解，R矩陣由相鄰兩觀測點(diǎn)模型誤差之差的平方和的形式確定，仍設(shè)計(jì)以下2 種方案：方案1，采用L 曲線法確定正則化參數(shù)；方案2，采用U 曲線法確定正則化參數(shù)。

這2種方案計(jì)算的模擬系統(tǒng)誤差的估值與殘差估值分別見圖7和圖8，計(jì)算所得的X^ 估值與模擬X真值差異的比較如表4所示。

從圖7(a)和圖8(a)可以看出：當(dāng)模擬的系統(tǒng)誤差呈線性周期性時(shí)，L曲線法和U曲線法確定的正則化參數(shù)都能使模擬的系統(tǒng)誤差與模擬系統(tǒng)誤差的估值較吻合，但U曲線法確定的正則化參數(shù)解算得到的系統(tǒng)誤差的估值與其真值的吻合程度也明顯比L曲線法的高。從圖7(b)和圖8(b)可以看出：L 曲線法和U 曲線法確定的正則化參數(shù)同樣都能使半?yún)?shù)估計(jì)殘差與模擬殘差變化的趨勢基本一致，但是吻合程度不高。這是因?yàn)殡S著加入系統(tǒng)誤差的復(fù)雜，半?yún)?shù)模型解算難度增加。但兩者相比而言，采用U曲線法確定的正則化參數(shù)吻合程度要高。

從表4可以看出：L 曲線法和U 曲線法所確定的正則化參數(shù)差異較大，但都能使所求結(jié)果與真值較接近，均提高了半?yún)?shù)模型解的精度。通過進(jìn)一步分析發(fā)現(xiàn)：采用U曲線法得到的α為0.095 5，X的估值為=[4.995 7 2.001 0]T，與其模擬X真值差異為=4×10-4；采用L曲線法得到的α為0.715 4，與其模擬X真值差異為=7×10-4。與L 曲線法相比，U曲線法確定的正則化參數(shù)求得參數(shù)估值更接近于真值。

實(shí)驗(yàn)結(jié)果表明，L曲線法和U曲線法都可以較好地分離觀測值中的系統(tǒng)誤差，求得參數(shù)的估值；當(dāng)模擬的系統(tǒng)誤差呈周期性或者線性周期性時(shí)，U曲線法確定的正則化參數(shù)求得參數(shù)估值的精度要比L曲線法的高。

表3 2種方案計(jì)算的周期性誤差下估值與模擬X真值差異的比較Table 3 Comparison of calculatedand simulated valuesX of 2 schemes under periodic error

表3 2種方案計(jì)算的周期性誤差下估值與模擬X真值差異的比較Table 3 Comparison of calculatedand simulated valuesX of 2 schemes under periodic error

方法L曲線法U曲線法αX^σ0 ΔX^=‖ ‖1.438 4 1.398 1[1.018 0 0.238 2]T[1.011 3 0.235 1]T 1.021 4 0.967 7X^- X 4.632 4×10-4 3.497 0×10-4

圖7 L曲線法確定正則化參數(shù)α(線性周期系統(tǒng)誤差)Fig.7 Determining regularization parameterα according to L-curve method(system error of linear periodic)

圖8 U曲線法確定正則化參數(shù)α(線性周期系統(tǒng)誤差)Fig.8 Determination of regularization parameterα according to U-curve method(system error of linear periodic)

表4 2種方案計(jì)算的線性周期性誤差下估值與模擬X真值差異的比較Table 4 Comparison of calculated and simulated valuesX of 2 schemes under linear periodic error

表4 2種方案計(jì)算的線性周期性誤差下估值與模擬X真值差異的比較Table 4 Comparison of calculated and simulated valuesX of 2 schemes under linear periodic error

方法L曲線法U曲線法αX^σ0 ΔX^=‖ ‖0.715 4 0.095 5[4.973 8 1.968 1]T[5.006 5 2.019 8]T 0.953 8 0.509 8X^- X 7×10-4 4×10-4

4 結(jié)論

1) 鑒于半?yún)?shù)模型補(bǔ)償最小二乘估計(jì)中正則化參數(shù)合理確定的問題，研究了一種新的正則化參數(shù)確定方法即U曲線法。模擬算例分析表明采用U曲線法確定合適的正則化參數(shù)，能夠有效地控制VTPV與STRS之間的平衡，得到了較準(zhǔn)確的參數(shù)估值。

2) 當(dāng)模擬的系統(tǒng)誤差分別呈周期性、線性周期性時(shí)，采用U曲線法確定的正則化參數(shù)所求得參數(shù)估值和其真值差值向量的范數(shù)與L曲線法確定的范數(shù)相比分別提高了1.135 4×10-4和3.000 0×10-4。

3) 由于在實(shí)際問題中，非參數(shù)部分即系統(tǒng)誤差或者模型誤差具有一定的復(fù)雜性，會(huì)影響正則化參數(shù)的選擇。本文提出的U曲線法在一定條件下求得正則化參數(shù)使參數(shù)估計(jì)解的精度較高，在不同條件下，正則化參數(shù)的確定要根據(jù)相應(yīng)的情況具體選取合適的方法。