利用TSVD参数估值变化特性确定算法截断参数

林东方，姚宜斌，郑敦勇，李朝奎

1. 湖南科技大学测绘遥感信息工程湖南省重点实验室，湖南 湘潭 411201； 2. 武汉大学测绘学院，湖北 武汉 430079； 3. 湖南科技大学地理空间信息技术国家地方联合工程实验室，湖南 湘潭 411201； 4. 武汉大学地球空间环境与大地测量教育部重点实验室，湖北 武汉 430079

受观测条件限制，在地球重力场反演、GNSS空间环境测量及InSAR地表形变测量等大地测量应用领域常会出现病态问题[1-7]。病态问题导致模型解算稳定性较差，难以得到准确可靠的模型参数估值，如何处理与解算病态问题成为大地测量数据处理的重要研究内容[8-9]。

病态问题具体表现于函数模型的奇异性及模型参数估计的扰动性上[10]。病态函数模型设计矩阵往往包含较小的甚至接近于零的奇异值，该奇异值导致参数估值方差较大，观测数据中的微小误差即引起参数估值的剧烈扰动，这种情况下，常规最小二乘估计难以得到模型参数的准确估值[11]。为了提高估值的稳定性和可靠性，文献[12—15]在最小二乘估计基础上提出了正则化法、岭估计法、截断奇异值法(truncated singular value decomposition，TSVD)等有偏估计方法。在均方误差意义下，3种方法均通过增加偏差、减少方差的形式降低模型参数估值均方误差，然而如何确定偏差与方差平衡点，以最大程度降低均方误差，仍是3种方法尚未解决的难题[16-17]。TSVD通过截掉小奇异值来降低方差，是一种较为简捷高效的病态问题解算策略，在诸多领域得到了广泛应用[18-19]。影响TSVD模型参数估值均方误差的关键因素是截断参数，目前，常用的截断参数确定方法主要有L曲线法、均方误差最小法等，L曲线法通过计算模型参数估值二范数与观测值残差二范数变化曲线的拐点来确定截断参数。L曲线法确定的截断参数没有明确的理论依据，曲线拐点处的参数估值并不代表最优的参数估值[20-21]。均方误差最小法通过计算各截断参数下的模型参数估值均方误差最小值确定截断参数[22]，该方法理论依据明确，但均方误差的计算需要利用模型参数真值，而真值是未知的，以估值代替真值计算的均方误差与真实均方误差存在一定差异，限制了TSVD解算效果。

均方误差反映了模型参数的估计质量，对于有偏估计，均方误差由方差与偏差两部分组成，TSVD截掉小奇异值引起参数估值方差与偏差的变化(方差减小，偏差增大)，该变化将反映在模型参数估值的变化上。鉴于此，本文利用参数估值变化分析奇异值截掉后的偏差变化，结合奇异值截掉后的方差变化确定最优截断参数，克服参数真值未知偏差难以计算问题，并通过试验验证本文方法的可行性与有效性。

1 病态问题解算的TSVD方法

1.1 TSVD解算方法

TSVD方法是在最小二乘估计算法基础上发展的病态问题解算方法，最小二乘估计准则表示为[23-25]

Φ=VTPV=min

(1)

式中，V=AX-L，表示观测值残差向量；A表示设计矩阵；L为观测向量；X表示未知模型参数；P为权重矩阵；由最小二乘估计准则可得模型参数估值为

(2)

若函数模型病态，则设计矩阵A存在较小的奇异值，对A进行奇异值分解可得

A=USGT

(3)

(4)

式中，U表示左奇异向量矩阵；S表示奇异值矩阵；G表示右奇异向量矩阵；γ表示奇异值，γ1>γ2>…>γn>0。由此可得模型参数估值方差

(5)

由式(5)可知，小奇异值导致参数估值方差较大，严重降低了模型参数估计精度[11，22]，TSVD算法通过截掉部分小奇异值来降低方差，提高模型参数估计稳定性，具体表示为

(6)

式中，G与U均取前t列向量。

(7)

式中，t表示由截断参数确定的需保留的t个较大奇异值。

1.2 常用TSVD截断参数确定方法

影响TSVD解算效果的关键因素是截断参数，最优截断参数可最大程度提高模型参数的估计精度。目前，最为常用的截断参数确定方法有L曲线法、均方误差最小法等。

1.2.1L曲线法

(8)

式中，ρ′、η′为一阶导数，ρ″、η″为二阶导数，以不同截断参数计算L曲线曲率，曲率最大点km对应的截断参数即为最优截断参数。L曲线法缺乏合理的理论依据，所确定的截断参数稳定性较好，但难以给出最优的截断参数。

1.2.2 均方误差最小法[22]

均方误差是模型参数估值与真值之间差值的数学期望，反映了参数估值相对于真值的离散程度。TSVD模型参数估值均方误差可表示为

(9)

由式(9)可见，均方误差的计算需要模型参数的真值，但实际应用中，参数真值是未知的。鉴于此，均方误差最小法以模型参数估值代替真值计算均方误差。第1步，利用TSVD方法截掉1个最小奇异值获得模型参数初步估计值；第2步，以初步估值代替真值计算不同截断参数下的均方误差，由最小均方误差确定最优截断参数，再次估计模型参数，得到该截断参数下的模型参数估值；第3步，将模型参数第2步估值替代真值重新计算不同截断参数下的均方误差，确定最优截断参数，如此迭代计算，直至两次模型参数估值二范数收敛到某一较小值[22]

(10)

然而，以估值代替真值计算的均方误差并非实际的均方误差，由此确定的截断参数受初值影响较大，常导致迭代过程不收敛或收敛于初值，限制了截断参数确定可靠性。

1.3 改进TSVD截断参数确定方法

1.3.1 截断参数确定方法

有偏估计均方误差由方差与偏差两部分组成，TSVD通过截掉小奇异值来降低方差，但截掉小奇异值导致模型偏离，引起模型参数估计偏差。TSVD模型参数估计均方误差可扩展为[22]

(11)

式中，Tt为TSVD模型参数估值方差总和；bt为估值偏差。通过奇异值分解，方差与偏差可表示为

(12)

(13)

式中，gi表示右奇异向量矩阵的第i列向量。

由式(12)和式(13)可知，TSVD每截掉1个小奇异值后，方差的减少量和偏差的增加量不尽相同。在截掉小奇异值后，方差的减少大于偏差的增加，则均方误差会降低，参数估值精度得到提高，估值更接近于真值。因此，判断某小奇异值是否需要截掉的关键在于该奇异值截掉后，方差减少量和偏差增加量的大小关系，即

(14)

分析方差与偏差对模型参数估值影响可知，TSVD每截掉一个小奇异值，模型参数估计方差与偏差会产生不同变化，这种变化将直接体现在参数估值变化上。由此，可计算出相邻奇异值截掉后的方差变化量与模型参数估值变化量，其中，模型参数估值变化量应包含截掉后一奇异值所引起的方差变化影响及偏差影响。在方差变化已知的情况下，参数估值变化量应可有效反映出偏差的变化。

(15)

由此，可建立奇异值截断方法

(16)

本文方法利用去除方差影响的TSVD模型参数估值变化近似描述奇异值截掉后的偏差变化，有效避免利用参数真值计算偏差，不受参数真值未知情形影响，与实际应用相符，具备较好的可行性与实用性。

1.3.2 标准差与参数估值变化确定方法

1.3.2.1 参数估值标准差变化量

标准差是方差的平方根，在无偏估计的情况下，标准差即是均方根误差，反映了参数估值与真值之间的差异。由式(12)可知，标准差的计算需要单位权方差，单位权方差反映了观测数据的观测精度，在仪器观测精度已知的情况下，可由仪器精度计算得到。在仪器观测精度未知时，可利用多余观测，通过无偏估计计算得到。

由最小二乘估计可得观测值残差向量

V=A(ATPA)-1ATPL-L

(17)

对设计矩阵A进行奇异值分解化简可得

(18)

式中，Um表示对应于奇异值的m×n阶左奇异向量矩阵，则单位权方差可通过式(19)估计为

(19)

由式(19)可知，单位权方差的估计与多余观测数有关，与奇异值大小无关，因此，在包含多余观测的情况下，可实现单位权方差的估计，利用单位权方差估值即可得到截掉奇异值后的标准差变化量

(20)

1.3.2.2 参数估值变化量

病态性对最小二乘估计的影响主要体现在小奇异值对参数估值方差的严重扩大。TSVD通过截掉小奇异值来消除其对方差的影响。截掉小奇异值降低了参数估值方差，但同时向参数估值中引入了偏差。因此，截掉某小奇异值后，方差和偏差的变化共同引起参数估值的变化。依次截掉相邻两个奇异值后的参数估值变化可计算为

(21)

从参数估值变化量中除去标准差变化影响即可估计出偏差变化影响

(22)

由式(20)和式(22)可知，截掉奇异值后的标准差变化量和偏差影响量均可通过式(22)计算得到。因而，该截断参数确定方法的理论依据更为明确客观，容易实现。但是，该方法依然存在缺陷，即在观测精度未知情况下，单位权方差的估计需利用多余观测来实现，多余观测的质量和数量决定了单位权方差估计的可靠性，进而影响到标准差计算的准确性。因此，该方法更适用于包含丰富多余观测的情形。

2 试验分析

2.1 空间测量试验

采用空间测边网算例进行试验分析，算例中包含2个未知点，9个已知点，通过19个等精度观测确定未知点坐标，观测中误差0.01 m。未知点坐标真值为A(0，0，0)与B(7，10，-5)，两未知点之间的约束观测值为13.185 9 m，已知点坐标及其他观测值情况见表1。

表1 空间测边网观测信息Tab.1 Spatial ranging information m

为了验证本文算法的有效性，针对上述测量问题，制定了两种解算方案：方案1利用18个观测值对两个未知点进行联合解算；方案2引入A、B约束观测，利用19个观测值联合解算两个未知点。

(1) 方案1。在A、B点联合解算时，设计矩阵的奇异值情况见表2，奇异值中包含较小接近于0的奇异值，设计矩阵存在病态问题。由图1可见，在截掉奇异值数为1时，参数估值相对于最小二乘估计变化量为2.55 m，标准差减少量为1.84 m。可见参数估值变化近70%是由标准差变化引起的，剩余30%则可认为是由偏差引起。因此，截掉最后1个奇异值，有利于降低参数估值均方误差。在截掉奇异值数为2时，参数估值变化量为0.22 m，标准差变化量为0.04 m，可见参数估值变化20%是由标准差变化引起，而剩余80%可认为由偏差引起，截掉该奇异值不利于降低均方误差。由此，本文方法确定的截断参数为1。L曲线法、均方误差最小法确定的截断参数见表3。

表2 A、B点联合解算设计矩阵奇异值Tab.2 Design matrix singular value of joint estimation

表3 不同方法模型参数估计结果Tab.3 Model parameter estimation results of different methods m

由表3可知，病态问题影响了最小二乘估计的参数估值精度。TSVD通过截掉小奇异值可有效改善最小二乘方法参数估计结果，降低参数估计误差。但L曲线法截掉3个小奇异值的坐标参数估计误差要大于均方误差最小法与本文方法截掉一个最小奇异值，可见均方误差最小法与本文方法确定的截断参数优于L曲线法，由此表明，两种方法均可得到较优的截断参数。

(2) 方案2。引入联测约束的设计矩阵奇异值见表4。由表4可知，引入联测约束后，观测方程设计矩阵病态性得到改善，对比表2，各奇异值均有所增大，病态性影响减弱。继续采用最小二乘与TSVD方法解算参数进行对比分析。本文截断参数确定方法参数估值及标准差变化情况如图2所示，截掉任一奇异值后的参数估值变化量均要远大于标准差变化量，即截掉奇异值后的偏差影响要大于方差影响，因此，不应截掉奇异值，本文方法确定截断参数为0。L曲线与均方误差最小法确定截断参数见表5。

表4 引入联测约束的设计矩阵奇异值Tab.4 Singular values of design matrix with constraint

图2 截掉奇异值后的参数估值及标准差变化Fig.2 Parameter estimates and standard deviation changes after truncating singular values

表5 不同方法模型参数估计结果Tab.5 Model parameter estimation results of different methods m

由表5可知，引入联测约束后，病态问题对最小二乘估计的影响减弱，最小二乘方法可得到参数的可靠估值。采用TSVD方法进行解算，通过L曲线法、均方误差最小法确定的截断参数，TSVD分别需要截掉3个和1个奇异值，但截掉奇异值引入偏差对参数估值的影响要大于降低方差，从而降低了参数估值精度。采用本文方法确定的截断参数为0，无须截掉奇异值，则TSVD与最小二乘估计结果相同，参数估值精度最优。

TSVD是在最小二乘估计基础上建立的病态问题解算方法，其解算效果取决于截断参数的选择。综合分析两种方案的解算结果可知，L曲线法确定截断参数依赖于奇异值大小差异，差异较大时，则容易出现拐点，但该拐点不能保证模型参数估计质量；均方误差最小法易受最小奇异值截掉后模型参数估计质量的影响，初步估计结果往往决定了截断参数的选择；本文方法综合考虑了奇异值截掉后的方差变化及参数估值变化(即方差与偏差影响)，截断参数确定依据充分，两种方案下均给出了合理的截断参数，有效提高了TSVD模型参数估计精度。

2.2 PolInSAR植被高反演试验

PolInSAR是大地测量中具备穿透测量能力的新兴热门测量技术，已在大范围植被高度与林下地形测量中得到了广泛应用。PolInSAR通过多极化穿透观测，能够有效地获取植被覆盖区地表及植被体散射信息，为植被高度及林下地形测绘提供了可能。然而，PolInSAR多极化观测模型参数之间存在一定的相关性，导致利用多极化数据进行植被高反演时常出现病态问题[27-29]，限制了植被高度的反演精度。为了验证本文方法的可行性与有效性，选取了德国宇航局BioSAR2008项目的E-SAR P波段多极化数据进行植被高反演试验，观测数据信息见表6。此外，试验区拥有高精度LiDAR植被高测量数据，可用于对比分析PolInSAR植被高反演结果。

表6 观测数据参数信息Tab.6 Information of multi-polarization observation

散射模型是刻画雷达波在植被覆盖区穿透传播过程的物理模型，是利用PolInSAR多极化观测信息反演植被高度的基础，随机地体二层散射(RVoG)模型是目前应用最为广泛的散射模型，该模型有效建立了极化观测量与植被参数之间的函数关系。具体表达为

(23)

式中，ω表示极化散射状态；γ(ω)为对应于极化状态ω的复相干系数，为已知观测值；φ0表示地表相位，为未知参数；μ(ω)表示地体幅度比，为未知参数；γv表示纯体相干性，与植被高参数相关联。γv具体表达为

(24)

式中，σ表示消光系数，为未知参数；θ表示雷达入射角，为已知值；hv表示植被高参数，为未知参数；kz表示垂直向有效波数，为已知值。极化观测γ(ω)及纯体相干性γv为复数值。由散射模型可构建函数模型

γ(ω)=f(φ0,σ,hv,μ(ω))

(25)

由式(25)结合复数最小二乘估计准则，可构建高斯-马尔可夫(Gauss-Markov，G-M)模型[28-29]

Vγ=AγXγ-Lγ

(26)

式中，Vγ表示残差向量；Aγ为系数矩阵；Xγ为模型参数改正数向量；Lγ表示实部与虚部残余常数向量。利用式(26)实现植被高参数平差估计。

本次试验利用HH、HV、VV、HHpVV、HHmVV、opt1、opt2、opt3、PDHigh、PDLow这10种极化方式观测数据，进行植被高参数估计。误差方程共包含植被高、地体幅度比等13个未知模型参数，设计矩阵奇异值情况见表7。

由表7中的奇异值情况可知，设计矩阵存在严重的病态性。图3为奇异值截掉后的TSVD参数估值及标准差变化情况。由于最小奇异值过小，截掉后标准差变化过大，参数估值及标准差变化趋势拆分为图3(a)、图3(b)两部分展示。由图3可知，前2次截掉小奇异值，方差减少影响均要大于偏差增加影响，有利于降低模型参数估值均方误差。在第3次截掉奇异值后(截掉奇异值数为3)，方差减少对参数估值的影响要小于偏差增加对其影响，第3次截掉奇异值不利于均方误差的降低。因此，本文方法确定截断参数为2。为了对比分析，分别采用L曲线法、均方误差最小法确定截断参数，解算观测方程，获得植被高、地体幅度比等模型参数的TSVD估值。表8为不同方法模型参数估计的结果。

表7 观测方程设计矩阵奇异值Tab.7 Design matrix singular values of observation equation

图3 截掉奇异值后的参数估值及标准差变化Fig.3 Parameter estimates and standard deviation changes after truncating singular values

表8 不同方法模型参数估计结果Tab.8 Model parameter estimation results of different methods

由于函数模型中包含一些过程参数无参数真值进行对比分析，因此表8仅给出了植被高及影响植被高反演的地体幅度比参数的估计结果。由式(23)可知，地体幅度比参数影响了纯体相干性的相位高度，相位高度越接近于植被冠层，则植被高反演越准确。由表8中两类参数的反演结果可知，3种TSVD解算方法有效改善了最小二乘方法参数估计结果。相较于LiDAR植被高测量结果(20.02 m)，L曲线法确定截断参数时，TSVD植被高估值偏低，这与地体幅度比的估值结果相符，因为地体幅度比估值相较于其他两种方法也偏低，则纯体相干性相位高度位于植被冠层以下，从而造成植被高参数低估。由均方误差最小法确定截断参数时，植被高与地体幅度比估值均偏高，这与L曲线法相反，过高的地体幅度比导致纯体相干性相位高度位于植被冠层以上，造成植被高参数高估。采用本文方法确定截断参数时，尽管植被高与地体幅度比估值仍偏高，但相较于其他两种方法有明显改善，植被高估计结果最接近于LiDAR结果，精度最高。便于直观对比分析，整幅数据的植被高估计结果如图4所示。

由图4可以看出，受模型病态性影响，最小二乘估计已无法获得植被高参数的可靠估值。采用TSVD方法进行解算，不同截断参数下的植被高反演结果存在较大差异。由L曲线法确定截断参数，TSVD方法反演的植被高结果相较于三阶段初值未有明显改善，这主要由于L曲线法截掉的奇异值过多，出现过度平滑，解算结果未能得到改善。均方误差最小法确定截断参数，植被高反演结果相较于初值有明显改善，但与LiDAR结果相比，存在一定高估。这是由于均方误差最小法仅截掉一个最小奇异值，虽然病态性得到了较大改善，但剩余小奇异值病态性仍较为严重，植被高反演均方误差仍较大，导致存在高估。由本文方法确定截断参数，TSVD表现最优，植被高反演结果最接近于LiDAR植被高结果，表明新截断参数确定方法可有效改善TSVD解算效果，提高模型参数估计精度。为了量化分析，由图中均匀选取1377块样地，统计分析植被高反演误差情况，统计结果如图5所示。

图4 各方法植被高反演结果Fig.4 Vegetation height inversion results of each method

由图5可知，最小二乘估计均方根误差较大且离散程度较高，估计精度较差，无法获得植被高可靠估值。由L曲线法确定截断参数的TSVD方法，植被高估计结果离散度有较大改善，均方根误差相较于最小二乘估计有显著降低；但相较于其他方法仍较大。均方误差最小法确定截断参数，TSVD植被高估值均方根误差优于L曲线法，但差于本文方法。由本文方法确定截断参数，TSVD植被高估值均方根误差均低于其他方法，与LiDAR样地植被高符合程度最高，表明本文方法确定截断参数的TSVD植被高估计结果最佳，验证了本文方法在改善TSVD解算效果上的可行性与有效性。

图5 各方法样地植被高反演误差Fig.5 Standards vegetation height inversion errors of each method

3 结 论

TSVD方法是处理病态问题的常用有效方法，在大地测量各领域得到了广泛应用。截断参数的选择决定了TSVD方法的解算效果，常用的截断参数确定方法从不同角度给出了有效的截断参数，但仍难以确定最优截断参数。均方误差最小法考虑了TSVD模型参数估值均方误差的变化，可在均方误差理论下保证模型参数估值精度，是一种理论依据较为完备的截断参数确定方法，但均方误差的计算需要模型参数真值，在实际应用中无法满足，导致该方法难以确定出理论上的最优截断参数。鉴于此，本文提出了考虑均方误差视角下参数估值变化差异的TSVD截断参数确定方法，利用奇异值截掉后的方差与参数估值变化关系分析确定偏差影响，综合方差与偏差影响，最终实现基于均方误差最小准则的截断参数确定。模拟与实际应用的试验结果表明，本文方法确定的截断参数可有效改善TSVD解算效果，提高模型参数估值的精度与可靠性。