病态问题常存在于卫星定位^[1-2]、回归分析^[3-4]、重力场延拓^[5-7]、PolInSAR植被参数反演等测绘领域。由于函数模型的病态性，常规最小二乘估计(least squares，LS)的方差较大，参数估值极不可靠。这种情况下，LS虽然仍是无偏估计，但已不是最优估计^[8]。函数模型的病态性体现在观测方程系数矩阵中出现较小甚至接近于零的奇异值，导致参数估值的方差被小的奇异值严重扩大，造成估值精度的降低。为了提高估值的稳定性和精度，学者们提出了一系列改善估计质量的有偏估计方法，如Tikhonov正则化法、岭估计法、截断奇异值法(truncated singular value decomposition，TSVD)。Tikhonov正则化法在最小二乘估计准则的基础上加入稳定泛函约束以提高估计的可靠性，并引入正则化参数调节两者的平衡^{[9, 10]}。岭估计法与Tikhonov正则化法的估计准则类似，可看作是Tikhonov正则化法的简化形式，岭估计将稳定泛函确定为待定参数的二范约束，即将稳定泛函中的正则化矩阵替换为单位矩阵^[11-12]。截断奇异值法则是基于奇异值分解技术的直接解算方法，将严重扩大方差的小奇异值截去，仅利用较大奇异值包含的观测信息进行参数估计^[13-15]。3种有偏估计方法均有效改善了病态问题参数估计的可靠性和稳定性，是应用较为广泛的解算方法。

相比于Tikhonov正则化法和岭估计法，TSVD仅需要考虑截断小奇异值对估计的影响，是一种更为简单直接的解算方法，在大地测量各领域得到了广泛应用。影响TSVD解算效果的关键因素是截断参数的选择^[16]。目前，确定奇异值截断参数的方法主要有广义交叉核实法^[17-18](generalized cross-validation，GCV)和L曲线法^[19-20]。GCV法通过交叉验证的方式选取截断参数使检验集观测值的残差平方和最小^[17]。然而对于有偏估计，观测值残差平方和最小并不能保证参数估值的均方误差最小。L曲线法通过计算TSVD参数估值二范数与观测值残差二范数变化曲线的拐点来确定截断参数^[20]，这种方式对于阶梯型变化的奇异值效果较好，而对于均匀下降型的奇异值，曲线拐点不明显，常导致L曲线法得不到可靠的截断参数。

本文针对TSVD截断参数的选取问题，首先分析了较小的奇异值对估值方差的不利影响；其次分析了截掉奇异值后参数估值方差与偏差的变化，当截掉奇异值后的方差下降量大于偏差引入量时，TSVD可提高估值的精度；进而研究提出了一种通过判断截掉奇异值后方差下降量与偏差引入量大小来确定截断参数的方法，更合理可靠地确定截断参数，提高TSVD估计精度。

1 病态性的影响及TSVD解算方法 1.1 病态性对LS的影响

测量数据处理常用的G-M模型表示为^[21]

式中，L为观测值向量；V为观测值残差向量；A为设计矩阵；X为未知参数向量。根据最小二乘准则

求得参数的估值为

其方差-协方差矩阵为

式中，P为权矩阵。对式(1) 作如下变换

式中，P为对角阵，表示对角线元素求平方根，令，式(3) 可变换为

因此，为了方便推导，可设权矩阵P为单位阵，对系数矩阵A进行奇异值分解^[22]

式中，U为左奇异向量矩阵；S为奇异值矩阵；γ₁>γ₂>…>γ_n为其奇异值；G为右奇异向量矩阵。协方差矩阵的迹是各参数估值的方差之和，可以整体上反映参数估计方差的大小，由式(4)、式(8) 得LS估值的方差之迹为

式中，tr(·)为矩阵的迹运算。如果方程病态，γ₁远大于γ_n，γ_n为接近于零的较小值。由式(9) 可以看出，较小的奇异值会对估值的方差造成严重影响，导致LS估计极不可靠。

1.2 解算病态问题的TSVD方法

TSVD法是基于奇异值分解技术的病态方程的一种直接解法，其原理是截断放大方差的较小的奇异值^[13]，保留较大的奇异值不变，来消除较小的奇异值对方差的影响，以降低估值的方差。

TSVD法通过截断参数保留较大的奇异值，截掉较小的奇异值。对设计矩阵A进行奇异值分解，得到奇异值矩阵S，如式(7)、式(8)。设由截断参数保留的奇异值为奇异值矩阵S的前k个，截掉的奇异值为S中的后n-k个，将奇异值矩阵S中保留的奇异值求逆，截掉的奇异值取零得

则病态方程的截断奇异值解法为^[13]

式中，X_α为TSVD参数估值，由式(10) 可以看出，截断奇异值法认为设计矩阵A较大的奇异值不会对估计方差造成不良影响，而较小的奇异值则会扩大方差。截掉较小的奇异值，可消除其对估值方差的放大影响，但改变了原设计矩阵的结构，导致估值有偏。因此，截断奇异值法是一种有偏估计方法，在改善矩阵病态性的同时引入了偏差。偏差的大小与截断参数息息相关，截掉的小奇异值越多，方差下降越大，而模型偏离也越大，引起的偏差越大，因而截断参数起到平衡估计方差与偏差的作用，是影响TSVD解算效果的关键因素。

1.3 常用截断参数确定方法

1.3.1 广义交叉核实法(GCV)

GCV法是由文献[17]提出的截断参数确定方法。从统计检验的角度确定截断参数，文献[17]给出了GCV函数，通过计算GCV函数最小值确定截断参数^[17]。

对于函数模型

GCV函数表示为

式中, α为截断参数；m为观测值的个数；H_α=AGS_I^TU^T为影响矩阵，满足AX_α=H_αL；I为单位阵。GCV法就是求解GCV函数的最小值所对应的α作为截断参数。GCV法统计理论较为完善，通过交叉验证的方式确定截断参数，但是GCV法以观测值残差作为检验条件并不能保证TSVD改善了参数估计的均方误差。

1.3.2 L曲线法

文献[14]提出用L曲线法确定TSVD解算方法的截断参数，其原理是以log‖X_α‖²为纵坐标，log‖AX_α-L‖²为横坐标绘制曲线，计算曲线上曲率最大的点，该点对应的α确定为截断参数。该曲线形似L，因此称为L曲线法。L曲线法的核心是认为L曲线上曲率最大点对应的截断参数即是最优的截断参数^[20-21]。

令η=2log‖X_α‖, ρ=2log‖AX_α-L‖，则L曲线上曲率k的计算公式可表示为

式中，ρ′、η′表示一阶导数；ρ″、η″表示二阶导数；计算式(14) 的最大值即可得到曲线上的最大曲率k_max、k_max对应的α即为确定的截断参数。L曲线法在病态方程解算中得到了广泛应用，针对阶梯型变化的奇异值效果显著，然而对于非阶梯型奇异值，容易发散或陷于局部最优。

目前，常用的截断参数确定方法GCV法和L曲线法未从TSVD改善模型参数估值质量的角度出发，稳定性和可靠性较差。文献[13]提出了基于均方误差最小的截断参数确定方法，理论依据更为充分，有效确定了截断参数，然而由于待定参数真值未知，均方误差估值的可靠性受限于参数的LS估计或岭估计的可靠性。本文基于文献[13]提出的均方误差最小确定截断参数的思想，通过分析截掉奇异值后方差与偏差的变化，尝试提出一种可靠的截断参数确定方法。

2 TSVD参数估计的方差与偏差分析及截断参数确定新方法 2.1 TSVD法的方差与偏差分析

TSVD是均方误差(MSE)优于LS估计的有偏估计方法，其均方误差反映了参数估值与参数真值之间的差异，其均方误差表示为^{[13, 23]}

式中，表示参数真值。有偏估计的均方误差由方差和偏差组成，因此，对TSVD的方差与偏差进行分析，可得到截断参数对TSVD解算效果的影响。

由式(11) 结合协方差传播律可得TSVD协方差矩阵迹为

由式(11) 可得TSVD的估计偏差与偏差平方和为

式中，g_i为奇异值γ_i的右奇异向量，是矩阵G的第i组基向量。比较式(9) 与式(17) 可以得出，TSVD截掉小奇异值可消除小奇异值对方差的扩大，有效减少估值的方差。由式(19) 可以看出，估值的偏差是由于截掉小奇异值所引起，截掉的小奇异值越多引起的偏差越大。分析式(17) 与式(19) 可以看出，TSVD截掉一个小奇异值降低的方差量为σ₀²/γ_i²，同时引入的偏差量为。因此，小奇异值的截掉与否与方差下降量和偏差引入量大小有关，在σ₀²/γ_i²大于，小奇异值需截掉，而σ₀²/γ_i²小于时，小奇异值应保留。鉴于此，可通过比较截掉奇异值的方差下降量与偏差引入量来确定截断参数。

2.2 病态方程奇异值截断新方法

由最小二乘估计可得观测值的残差为

对式(20) 中的系数矩阵A进行奇异值分解，化简得观测值残差与单位权方差估计量为

式中，W表示对应于奇异值的前n列左奇异向量；m为观测值的个数；n为待定参数的个数。由式(22) 可以看出，在观测值个数大于待定参数时，单位权方差可被估计出来且与奇异值无关，因此，由最小二乘估计得到的单位权方差估计值是可靠的，截掉奇异值的方差下降量可通过准确计算出。

对于偏差引入量的计算，由于参数X的真值是未知的，也无法准确求出。文献[13]提出采用X的最小二乘估值或其他估值替代X的真值。由式(15) 可以看出，X真值由估值替代时，均方误差实际上反映的是TSVD估值与真值替代值之间的离散程度，因此，真值替代值须具备较高的准确性才能得到可靠的均方误差值。然而，由于病态性的影响，X的准确估值往往难以得到。考虑替代值由最小二乘估计获得，对估计式(3) 进行奇异值分解化简得

式中，D为n×n奇异值对角矩阵包含S中的n个奇异值；W为m×n矩阵表示对应于奇异值的前n列左奇异向量。由式(8) 可以看出，原S矩阵中存在部分零矩阵，在运算过程中这部分零矩阵被化简掉，从而得到式(23)。式(23) 两边同乘以可得

式中

因此

式中，u_i表示第i列左奇异向量。由式(26) 可知，采用估值替代不能可靠地计算偏差，偏差量会受到小奇异值的影响，奇异值越小其值越大。此外，由式(26) 结合协方差传播律可得估计量的方差为

由式(27) 可以看出，的估计方差与截掉相应奇异值的方差下降量数值相同，在奇异值γ_i为较小值时，估值的方差较大，估值准确性较低，因此，采用X的最小二乘估值替代真值难以保证均方误差计算的可靠性。然而，在奇异值γ_i为较大值时，的估计方差较小，其估值是准确的，结合式(26) 可知，由于奇异值为大于零的正数，的估计方差随奇异值的变小不断增大，因此则随奇异值的变小而不断增大。考虑到特征向量矩阵G为一组标准正交基，是参数向量X在g_i方向的投影值，在参数向量X取真值时，其在正交向量组G各向量上的投影值虽大小不同，但必然在一定范围内波动，因此，也在一定范围内波动，而在参数向量X取最小二乘估值时，其各投影值波动范围会随着方差的增大不断增大，则呈现波动上升。因此，通过较大奇异值对应的投影估值波动范围可以一定程度上反映出整个向量组G上的投影值的波动范围。

较小的奇异值是影响参数估计可靠性的根源，当奇异值减小到一定程度时，由于其对方差的严重影响，无需考虑偏差引入多少，都是需要截掉的。因此，截掉过小的奇异值的偏差引入量也不需要准确估计，截断参数选择的关键在于判定中度大小的奇异值截掉后的方差下降量与偏差引入量大小。通过确定较大奇异值的偏差引入量及其波动范围，可有效反映出中度大小奇异值的偏差引入量波动范围，考虑到方差下降量σ₀²/γ₁² < σ₀²/γ₂² < … < σ₀²/γ_n²，随奇异值变小持续增大。当方差下降量超出此波动范围时，可认为截掉后续的小奇异值的方差下降量已大于偏差引入量，有利于减少估值的均方误差。由式(26) 可以看出，偏差引入量估值会受小的奇异值影响而被扩大，偏差引入量估值曲线随奇异值变小则呈现波动上升。然而由于较大奇异值的偏差引入量估值未被明显扩大，会呈现平缓波动，在实际计算中可容易判定。

设前j个奇异值为较大的奇异值，其偏差引入量估值是可靠的，可看作向量J= ，向量J中的元素大小在一定范围内波动，反映了偏差引入量数值的可信区间。由于J向量中的部分引入量也被一定程度的高估，实际应用中可设定方差下降量增大到大于J向量中50%以上元素时，截掉后续的小奇异值。

3 试验分析 3.1 模拟试验分析

Fredholm第一类积分方程是一类典型的病态问题，大地测量中的重力场延拓问题本质上就是对该方程的解算^[24-25]，该积分方程的表达式为

本文取其核函数为

精确解函数为

式中， [0,1]；y∈[-2,2]。

在采样间隔为Δx=Δy=0.02，对式(28) 积分方程进行离散化可得

式中，j=1, 2, …，201。由式(31) 可知，K(x_i, y_j)和Δx组成了维数为201×51的设计矩阵A，观测向量和参数向量可分别表示为：，设计矩阵A的条件数为2.5e+06，病态性严重。为了验证新方法的有效性，在观测向量中加入正态分布为N(0, 1.0e-04) 的随机误差，采用LS和TSVD方法进行解算，模拟试验500次得到解算结果的统计特性。

图 1反映了设计矩阵的奇异值大小变化情况，由图中可以看出，奇异值由大到小呈均匀下降，未出现明显断层，需应用合理的截断参数确定方法获取可靠的截断参数。

图 2给出了截掉奇异值后方差下降量和偏差引入量估值情况，由图中可以看出，随着奇异值变小，方差下降量持续增大，偏差引入量估值波动增大且呈上升趋势，表明小奇异值不仅对待定参数估值的方差影响较大，同时也影响了偏差引入量的估计，其估值随奇异值变小不断增大。然而在奇异值较大时，方差下降量较小，其对偏差引入量估值的影响也较小，表明较大的奇异值对应的偏差引入量估值是可靠的，因此图中矩形区域的偏差引入量估值波动平缓，未出现上升趋势。选取矩形区域奇异值的偏差引入量估值建立向量J，作为与方差下降量进行比较的数据集。实际计算中，可设定矩形末端偏差引入量估值方差为上限，选取估值方差小于方差上限的偏差引入量估值建立向量J。

图 3反映了偏差引入量估值与偏差引入量真值的对比情况。由图 3可以看出，各奇异值对应的偏差引入量真值在一定区间内波动，并非随着奇异值的增大而增大。而偏差引入量的估值在前端波动平缓，后端随着奇异值的变小，估值不断增大，可见小奇异值对应的偏差估值会受小奇异值的影响而被高估。前端波动平缓处的偏差估值与偏差真值基本在同一波动区间内，在此波动区间前端较大的奇异值，其偏差估值与真值相近，后端的奇异值偏差估值略有高估。因此，在实际计算中，可设定方差下降量大于平缓处偏差估值中50%以上数值时，确定截断参数，截掉后续的小奇异值。

为了比较分析新方法确定截断参数的可靠性和有效性，分别应用GCV法、L曲线法、LS估计偏差的最小MSE法和新方法确定截断参数，采用LS估计和TSVD方法进行解算的统计结果见表 1。

表中AD(absolute difference)表示绝对误差；RMSE(root mean square error)表示均方根误差。由表 1可以看出，由于病态性，最小二乘估计已不是参数的有效估计方法，参数估值严重偏离真值。采用TSVD方法进行解算，估计结果得到有效改善，参数估值接近真值，但截断参数的选取严重影响TSVD的解算效果。由表 1可以看出，采用GCV法可得到有效的截断参数，但估计误差的最大值较大，表明GCV法稳定性较差，常得到不合理的截断参数，导致误差较大，影响估计的可靠性。采用L曲线法确定的截断参数稳定性较好，估计误差的最大值与最小值相差较小，但估计误差的最小值与平均值明显大于其他方法，表明L曲线法容易发散，得到的截断参数非最优的截断参数。LS估计偏差的最小MSE法在本试验中未能有效的确定截断参数，估值误差较大。采用新方法确定的截断参数在估计误差的最大值、最小值和均值上均优于GCV法和L曲线法，且误差最大值与最小值相差较小，表明新方法具有较高的稳定性和可靠性，是一种有效的截断参数确定方法。3种截断参数确定方法的各次试验均方根误差情况绘于图 4。

由图 4(a)可以看出，GCV法确定截断参数稳定性较差，常给出不合理的截断参数，导致TSVD估计误差曲线波动较大，参数估值不可靠。L曲线法确定的截断参数具有较高的稳定性，给出的截断参数使TSVD估计误差曲线变化平缓，然而L曲线法确定的截断参数明显非最优的截断参数，表明L曲线法容易发散，尤其针对均匀下降型奇异值，难以得到合理的截断参数。相比于GCV法和L曲线法，新方法确定的截断参数具有较高的稳定性和可靠性，TSVD估计误差曲线变化平缓，波动幅度明显低于GCV法，且估计误差显著小于L曲线法。由图 4(b)可知，LS估计偏差的最小MSE法在本例中未能有效改善TSVD估计。因此，新方法是一种可行且合理的截断参数确定方法，可有效提高TSVD解算病态问题的可靠性和稳定性。

3.2 实测试验分析

利用PolInSAR数据进行植被高度反演时常会出现病态问题，严重影响了植被高反演的精度和可靠性^[26-27]。为了验证本文方法的可行性和有效性，选取了BioSAR2008项目的E-SAR全极化数据进行植被高度反演，应用TSVD方法解算反演中的病态问题，并分别采用GCV法、L曲线法、LS估计偏差的最小MSE法和新方法确定截断参数进行对比分析。

本试验采用文献[26]提出的观测值选取方法，选取HH、HV、VV、HHpVV、HHmVV、MCD(相干最大分离算法：opt1、opt2、opt3) 与PD(相位分离算法：PDhigh、PDlow)10种极化方式对应的复相干性作为观测值。数据解算中每个像素需要估计的参数为13个，观测方程20个，则观测方程系数矩阵的奇异值为13个。各奇异值对应的方差下降量与偏差引入量估值情况见表 2。

由表 2可以看出，第13个奇异值为过小的奇异值，由于其对方差的严重影响，无论偏差引入多少都需要截掉。前10个奇异值的偏差引入量估值在同一区间内波动，未出现明显增大，因此可设定方差上限为0.2，选择方差小于0.2的偏差估值建立向量J，当方差下降量增大到大于J中50%以上的元素时，确定截断参数，截掉后续的小奇异值。由于同一数据源的观测精度大致相同，各像素可设定同一方差上限，对解算效果的影响较小。

图 5中给出了不同方法反演的植被高结果。为了便于分析，图 5(a)给出了PolInSAR植被高反演常用的三阶段算法的反演结果，图 5(f)给出了LiDAR植被高测量结果作为植被高真值进行对比分析。由图 5(b)可以看出，截断参数由L曲线法确定时，TSVD反演结果与三阶段法的反演结果相似，有效改善了病态性的影响，但未能提高反演精度。而图 5(c)和图 5(d)的植被高反演结果均在右端出现了过度高估。相比于其他结果，图 5(e)中植被高反演结果更接近于LiDAR结果，表明采用新方法确定截断参数可有效改善TSVD的解算效果，提高参数估计精度。

由于LiDAR结果为平滑后的植被高度，笔者选取了800个30×30像素的样地植被高反演结果与LiDAR结果进行比较，统计结果见表 3。

由表 3可以看出4种TSVD解算结果均要优于传统三阶段方法结果，表明TSVD方法有效解决了文献[26]方法中出现的病态问题。其中由新方法确定截断参数的TSVD植被高反演结果最优，表明新方法有效改善了TSVD的解算效果，提高了植被高的反演精度。

由植被高反演误差分布统计图 6可以看出，三阶段法以及L曲线确定截断参数的TSVD法相比于其他3种方法，误差分布较为分散，多分布于-8~3 m，表明两种方法反演植被高的稳定性较差。而截断参数由GCV法和LS估计偏差的最小MSE法确定时，估计误差集中分布于1~5 m，表明植被高被过度高估。在截断参数由新方法确定时，估计误差集中分布于-1~3 m，可见新方法有效提高了植被高反演的可靠性和稳定性。因此，新方法是选取TSVD截断参数的一种可行有效方法。

4 结语

TSVD方法是一种简单有效的病态问题解算方法，相比于正则化方法，仅需考虑截断参数的影响，解算病态问题更直观高效，在大地测量各领域得到了广泛应用。然而，常用的截断参数确定方法GCV法和L曲线法均未从TSVD改善参数估值质量的角度出发，所确定的截断参数稳定性和可靠性较差，影响了TSVD解算效果，文献[13]提出的最小MSE法理论依据较为充分，但受限于参数真值未知，MSE难以准确计算。本文分析了TSVD由小到大依次截掉奇异值后估计方差与偏差的变化情况，提出了依据降低方差量小于引入偏差量确定截断参数的思想，理论依据充分。针对参数真值未知截掉奇异值所引入的偏差难以准确计算的问题，提出通过估计出较大奇异值截掉后的偏差引入量建立偏差估值可信区间，利用可信区间内偏差引入量与各奇异值截掉后的方差下降量进行比较，从而避免了较小奇异值截掉后的偏差引入量的计算，实现了任意奇异值截掉后的方差下降量与偏差引入量的合理比较。试验表明，相比于GCV法、L曲线法和最小MSE法，新方法确定的截断参数具有更高的稳定性和可靠性，可有效改善TSVD解算病态问题的效果。