CN107613458A - 一种tdoa条件下最优联合时间同步与定位的定位方法 - Google Patents

Abstract

本发明提出了一种TDOA条件下最优联合时间同步与定位的定位方法，用于解决现有移动端定位技术中存在的定位精度低的技术问题。实现步骤为：设定系统参数；建立联合时间同步与定位到达时间差TDOA的观测模型；基于联合时间同步与定位到达时间值TDOA的观测模型，建立双约束的移动端位置估计目标函数；求解双约束的移动端位置估计目标函数，得到待估向量y的全局最优解，并从中读出移动端位置向量p的全局最优解。本发明避免了定位误差的传递，获得了定位非凸问题的全局最优解，提高了定位精度，可用于视距场景下依靠多个蜂窝基站对移动端定位估计。

Description

一种TDOA条件下最优联合时间同步与定位的定位方法

技术领域

本发明属于无线通信技术领域，具体涉及一种TDOA条件下最优联合时间同步与定位的定位方法，可用于视距场景下依靠多个蜂窝基站对移动端定位估计。

背景技术

目前定位技术在无线通信、地质勘测、智能交通、抢险救援等方面扮演着举足轻重的角色，研究精确的定位方法具有实际的意义。定位按照是否采用测距方式分为基于测距的定位方法和非基于测距的定位方法两大类。基于测距的定位方法主要使用到达时间差TDOA测距法，到达时间测距法，接收信号强度测距法，以及到达角度测距法实现测距。目前基于TDOA的测距方法都是假设移动端与基站是在精确同步的情况下进行的，但是在实际环境中，往往存在时钟偏差，因此很难做到精确同步，并且这种方法会得到一系列的关于移动端位置的定位方程，这些方程具有高度的非线性和非凸性，很难求解出全局最优解，只能做到局部收敛，从而对定位精度产生严重的影响。例如：授权公告号为CN 102811419 B，名称为“一种基于迭代的最小二乘定位方法”的中国专利，公开了一种基于迭代的最小二乘定位方法，实现步骤为：首先对蜂窝基站数据进行分组；再利用球面相交(SSI)技术进行单次SSI-LS估计，得到中间估计值；利用中间估计值计算残差，得到相应的权值，并归一化加权，得到移动站的位置初始估计值；将初始估计值代入残差Taylor法，通过迭代得到最终位置估计值。虽然该技术能提升抗NLOS能力，能够实现在蜂窝通信系统中可靠的移动站定位，但是由于忽略了中间变量与待估参数的相关性，使得局部收敛，未获得全局最优解，定位精度受到影响。

发明内容

本发明的目的在于克服上述现有技术存在的缺陷，提出了一种TDOA条件下最优联合时间同步与定位的定位方法，用于解决现有基于测距的定位方法中存在的定位精度低的技术问题。

本发明的技术思路是：针对现有基于测距的定位方法无法获得全局最优解的情况，利用时间观测差值建立联合时间同步与定位的TDOA模型，利用中间变量将非线性模型方程线性化，建立双约束的移动端位置估计目标函数，通过解目标函数，得到非凸定位问题下移动端位置向量的全局最优解。

根据上述技术思路，实现本发明目的采用的技术方案包括如下步骤：

(1)设定系统参数：

参与移动端定位且位置已知的蜂窝基站数量为m，移动端数量为1，基站与移动端均位于n维坐标系；移动端在外部时钟为0时发送信号给m个基站，选定坐标零点处的第1个基站作为参考基站，相对于外部参考时钟的时钟相偏为τ₁，接收到移动端信号的外部参考时钟为T₁，观测噪声为ΔT₁；其余各个蜂窝基站的位置向量为p_i，各个基站相对于外部参考时钟的时钟相偏为τ_i，各个基站接收到移动端信号的外部参考时钟为T_i，观测噪声为ΔT_i；待定位的移动端的位置向量为p，移动端相对于外部参考时钟的时钟相偏为τ，其中，m≥4，i＝2,3…m，p_i∈Rⁿ，p∈Rⁿ，τ₁∈R，τ_i∈R，τ∈R，T₁∈R，T_i∈R，R为实数域，位置向量p的单位为米，时钟相偏τ的单位为秒；

(2)建立联合时间同步与定位到达时间差TDOA的观测模型：

根据设定的系统参数，将联合时间同步与定位到达时间值TOA的观测模型中所有T_i分别与T₁作差，得到联合时间同步与定位到达时间差TDOA的观测模型：

其中，c为信号的传播速度，||p-p_i||为移动端与第i个基站的欧几里得距离，||p||为移动端与参考基站的欧几里得距离；

(3)通过联合时间同步与定位到达时间值TDOA的观测模型，建立移动端位置估计的有约束目标函数：

(3a)对公式(1)两端同乘c，得到关于距离差的观测模型：

r_i,1＝||p-p_i||-||p||+τ_1,i+Δr_i，1 (2)

其中，r_i,1＝cT_i-cT₁，τ_1,i＝cτ₁-cτ_i，Δr_i,1＝cΔT_i-cΔT₁，Δr_i,1为测量误差，其服从均值为零且方差为σ²的高斯分布，且σ²＞0；

(3b)将公式(2)转化为关于移动端位置向量p的线性等式：

对公式(2)两边分别求平方：

2r_i,1||p||+2p_i ^Tp-(p_i ^Tp_i-r_i,1 ²)＝2||p-p_i||Δr_i,1 (3)

定义一个与移动端位置向量p相关的中间变量β，β＝||p||，β≥0，并将其带入(3)式中，得到m-1组关于移动端位置向量p的线性等式：

(3c)将m-1组关于移动端位置向量p的线性等式转化为矩阵形式：

By-g＝ε (5)

其中，ε为m-1行的列向量，且A为m-1行n+2列的矩阵，且g为m-1行的列向量，且待估向量y＝[p β]^T，y∈Rⁿ⁺¹；

(3d)根据最小二乘算法，将(5)式转化为无约束目标函数：

(3e)为公式(6)添加两个约束条件得到移动端位置估计的双约束目标函数：

其中，I_n为n行n列的单位矩阵，0_n×1为n行1列的零矩阵，0_1×n为1行n列的行向量，y_n+1为y向量的第n+1个元素；

(4)求解移动端定位估计的双约束目标函数，得到待估向量y的全局最优解：

(4a)获取待估向量y的解算条件式：

设c(y)，q(y)均为y的函数，且c(y)＝y^TCy，q(y)＝||By-g||²，并将(7)式转化为：

令λ满足以下三个条件，λ∈R：

c(y)＝0 (10)

其中符号表示求一阶导数，符号表示求两阶导数；

将c(y)和q(y)带入公式(9)、(10)和(11)中，得到待估向量y的解算条件式：

(B^TB+λC)y＝B^Tg (12)

y^TCy＝0 (13)

其中，表示B^TB+λC为正定矩阵；

(4b)根据待估向量y的解算条件式，将移动端位置估计的双约束目标函数的求解问题转换为多项式零点的求解问题：

根据式(12)，计算待估向量y：y＝(B^TB+λC)^-1B^Tg (15)

设是关于λ的函数，令：

将(15)式带入公式(16)得：

将公式(16)代入(13)，得方程实现将移动端位置估计的双约束目标函数的求解问题转化为多项式零点的求解问题；

(4c)利用式(14)计算的定义域I₁，并证明在定义域I₁内严格单调递减；

(4d)根据函数在定义域I₁上严格单调递减，利用二分法在定义域I₁内对λ进行有限次的迭代，得到的零点值；

(4e)将的零点值代入式(15)，得到待估向量y的估计值,并从中读出y_n+1；

(4f)判断y_n+1≥0是否成立，若是，则y就为全局最优解，从中读出移动端的位置向量p的全局最优解；否则执行步骤(4g)；

(4g)将待估向量y的解算条件式(14)更换为：B^TB+λC最多有一个负特征值，并在该条件下计算定义域I₂；

(4h)在定义域I₂内求出方程的所有实数根λ₁,λ₂…λ_k，其中1＜k≤n；

(4i)将λ₁,λ₂…λ_k代入y＝(B^TB+λC)^-1B^Tg，解出待估向量y的多个估计值y₁,y₂…y_k，并将y₁,y₂…y_k代入q(y)＝||By-g||²，将使q(y)取得最小值的y_j作为全局最优解，从中读出移动端的位置向量p的全局最优解，其中1≤j≤k。

本发明与现有技术相比，具有如下优点：

本发明建立基于TDOA的联合时间同步与定位模型，避免了现有基于测距的定位方法中定位误差的传递，利用对距离差平方处理，通过添加中间变量使模型方程线性化，利用最小二乘，并构造关于待估变量的两个约束条件，建立双约束的移动端位置估计目标函数，并求解该目标函数，获取待估变量的全局最优解，避免了现有技术中待估变量的局部收敛，实现了定位非凸问题的全局最优解的获得，提高了定位精度。

附图说明

图1为本发明的实施流程图；

图2为本发明与现有基于迭代的最小二乘定位方法的定位精度仿真对比图。

具体实施方式

以下参照附图和具体实施例，对本发明作进一步详细描述。

参照图1.TDOA条件下最优联合时间同步与定位的定位方法，包括如下步骤：

步骤1)设定系统参数：

参与移动端定位且位置已知的蜂窝基站数量为m，移动端数量为1，基站与移动端均位于n维坐标系；移动端在外部时钟为0时发送信号给m个基站，选定坐标零点处的第1个基站作为参考基站，相对于外部参考时钟的时钟相偏为τ₁，接收到移动端信号的外部参考时钟为T₁，观测噪声为ΔT₁；其余各个蜂窝基站的位置向量为p_i，各个基站相对于外部参考时钟的时钟相偏为τ_i，各个基站接收到移动端信号的外部参考时钟为T_i，观测噪声为ΔT_i；待定位的移动端的位置向量为p，移动端相对于外部参考时钟的时钟相偏为τ，其中，m≥4，i＝2,3…m，p_i∈Rⁿ，p∈Rⁿ，τ₁∈R，τ_i∈R，τ∈R，T₁∈R，T_i∈R，R为实数域，位置向量p的单位为米，时钟相偏τ的单位为秒；

步骤2)建立联合时间同步与定位到达时间差TDOA的观测模型：

根据设定的系统参数，将联合时间同步与定位到达时间值TOA的观测模型中所有T_i分别与T₁作差，得到联合时间同步与定位到达时间差TDOA的观测模型：

其中各个基站接收到移动端信号的外部参考时钟的表达式为：

c为信号的传播速度，||p-p_i||为移动端与第i个基站的欧几里得距离，||p||为移动端与参考基站的欧几里得距离；

步骤3)通过联合时间同步与定位到达时间值TDOA的观测模型，建立移动端位置估计的有约束目标函数：

步骤3a)对公式(1)两端同乘c，得到关于距离差的观测模型：

r_i,1＝||p-p_i||-||p||+τ_1,i+Δr_i，1 (2)

其中，r_i,1＝cT_i-cT₁，τ_1,i＝cτ₁-cτ_i，Δr_i,1＝cΔT_i-cΔT₁，Δr_i,1为测量误差，其服从均值为零且方差为σ²的高斯分布，且σ²＞0；

步骤3b)将公式(2)转化为关于移动端位置向量p的线性等式：

对公式(2)两边分别求平方：

2r_i,1||p||+2p_i ^Tp-(p_i ^Tp_i-r_i,1 ²)＝2||p-p_i||Δr_i,1 (3)

定义一个与移动端位置向量p相关的中间变量β，β＝||p||，β≥0，并将其带入(3)式中，得到m-1组关于移动端位置向量p的线性等式：

步骤3c)将m-1组关于移动端位置向量p的线性等式转化为矩阵形式：

By-g＝ε (5)

其中，ε为m-1行的列向量，且A为m-1行n+2列的矩阵，且g为m-1行的列向量，且待估向量y＝[p β]^T，y∈Rⁿ⁺¹；

步骤3d)根据最小二乘算法，将(5)式转化为无约束目标函数：

步骤3e)为公式(6)添加两个约束条件得到移动端位置估计的双约束目标函数：

其中，I_n为n行n列的单位矩阵，0_n×1为n行1列的零矩阵，0_1×n为1行n列的行向量，y_n+1为y向量的第n+1个元素；

步骤4)求解移动端定位估计的双约束目标函数，得到待估向量y的全局最优解：

步骤4a)获取待估向量y的解算条件式：

设c(y)，q(y)均为y的函数，且c(y)＝y^TCy，q(y)＝||By-g||²，并将(7)式转化为：

令λ满足以下三个条件，λ∈R：

c(y)＝0 (10)

其中符号表示求一阶导数，符号表示求两阶导数；

将c(y)和q(y)带入公式(9)、(10)和(11)中，得到待估向量y的解算条件式：

(B^TB+λC)y＝B^Tg (12)

y^TCy＝0 (13)

其中，表示B^TB+λC为正定矩阵；

步骤4b)根据待估向量y的解算条件式，将移动端位置估计的双约束目标函数的求解问题转换为多项式零点的求解问题：

根据式(12)，计算待估向量y：y＝(B^TB+λC)^-1B^Tg (15)

设是关于λ的函数，令：

将(15)式带入公式(16)得：

将公式(16)代入(13)，得方程实现将移动端位置估计的双约束目标函数的求解问题转化为多项式零点的求解问题；

步骤4c)利用式(14)计算的定义域I₁，并证明在定义域I₁内严格单调递减：

计算的定义域I₁：

步骤4c1)令矩阵得到矩阵B^TB+λC与矩阵合同，其中I为n+1行n+1列的单位矩阵；

步骤4c2)将矩阵表示为(C,B^TB)；

步骤4c3)令得到矩阵与矩阵C合同；

步骤4c4)利用线性代数中矩阵与矩阵C合同定义，矩阵与矩阵C拥有相同的正负惯性，即矩阵与矩阵C拥有相同的正特征值个数，负特征值个数和零特征值个数，根据得到矩阵有1个负特征值，n个正特征值，即：

λ₁(C,B^TB)＞λ₂(C,B^TB)＞…＞λ_n(C,B^TB)＞0＞λ_n+1(C,B^TB)

其中，λ_j(C,B^TB)表示矩阵的第j个特征值，j＝1,2…n+2，λ_j∈R，,且λ_j降序排列；

步骤4c5)再次利用线性代数中矩阵B^TB+λC与矩阵I+λ(C，B^TB)合同定义，即矩阵B^TB+λC与矩阵I+λ(C，B^TB)拥有相同的正特征值个数，负特征值个数和零特征值个数，根据公式(14)，得到矩阵I+λ(C,B^TB)的特征值均为正数，即：

1+λλ_j(C,B^TB)＞0

其中，λ_j(C,B^TB)表示矩阵的第j个特征值，j＝1,2…n+2，λ_j∈R，,且λ_j降序排列；

步骤4c6)根据λ1₁(C,B^TB)＞λ₂(C,B^TB)＞…＞λ_n(C,B^TB)＞0＞λn_n+1(C,B^TB)，j＝1,2…n+1，解不等式1+λλ_j(C,B^TB)＞0的步骤如下：

当j∈[1,n]，

当j＝n+1，

取不等式(17)和(18)的交集得λ的定义域I₁：

证明在定义域I₁内严格单调递减：

步骤4c7)对公式(12)两端同时对λ求一阶导得：

步骤4c8)根据公式(9)和(12)，得：

对上式进行移项得：

步骤4c9)根据公式(16)，将的两端同时对λ求一阶导：

步骤4c10)根据矩阵B^TB+λC为正定矩阵，利用正定矩阵的定义得即在定义域I₁内严格单调递减。

步骤4d)根据函数在定义域I₁上严格单调递减，利用二分法在定义域I₁内对λ进行有限次的迭代，得到的零点值；

步骤4e)将的零点值代入式(15)，得到待估向量y的估计值,并从中读出y_n+1；

步骤4f)判断y_n+1≥0是否成立，若是，则y就为全局最优解，从中读出移动端的位置向量p的全局最优解；否则执行步骤(4g)；

步骤4g)将待估向量y的解算条件式(14)更换为B^TB+λC最多有一个负特征值，并在该条件下计算定义域I₂；

步骤4g1)计算定义域I₂的前四步与计算定义域I₁的前四个步骤完全一致；

步骤4g2)根据矩阵B^TB+λC最多只有一个负特征值，j＝1,2…n+1，且

λ₁(C,B^TB)＞λ₂(C,B^TB)>…>λ_n(C,B^TB)＞0＞λ_n+1(C,B^TB)得：

当λ≥0，j∈[1,n]时，1+λλ_j(C,B^TB)＞0恒成立，要使B^TB+λC最多有一个负特征值，需要满足不等式：

1+λλ_n+1(C,B^TB)＞0 (19)

解不等式(19)得：

当λ＜0，j＝n+1时，1+λλ_n+1(C,B^TB)＞0恒成立，当j∈[1,n]时，1+λλ₁(C,B^TB)＜1+λλ₂(C,B^TB)＜…1+λλ_n(C,B^TB)，要使B^TB+λC最多有一个负特征值，需满足1+λλ_j(C,B^TB)中有n-1个正数和一个负数，即需要满足不等式：

1+λλ₁(C,B^TB)＜0＜1+λλ₂(C,B^TB) (21)

解不等式(21)得：

取不等式(20)和(22)的交集得λ的定义域I₂：

步骤4h)在定义域I₂内求出方程的所有实数根λ₁,λ₂…λ_k，其中1＜k≤n：

步骤4h1)利用矩阵的相似对角化规则求解出一个正交矩阵P，该矩阵能够将矩阵B^TB和矩阵C同时对角化：

P^T(B^TB)P＝diag(γ₁,γ₂…γ_n+1)

P^TCP＝diag(δ₁,δ₂…δ_n+1)

其中P为n+1行n+1列的正交矩阵，diag(…)表示n+1行n+1列的对角矩阵，γ₁,γ₂…γ_n+1，δ₁,δ₂…δ_n+1为对角矩阵的对角线上的元素；

步骤4h2)将矩阵方程式转换为代数方程式：

其中f为n+1行的列向量，f＝P^TB^Tg，f_z为列向量f的第z个元素，1≤z≤n+1；

步骤4h3)根据公式(23)，利用高阶多项式方程求根法求出方程的所有实数根λ₁,λ₂…λ_k，其中1＜k≤n；

步骤4i)将λ₁,λ₂…λ_k代入y＝(B^TB+λC)^-1B^Tg，解出待估向量y的多个估计值y₁,y₂…y_k，并将y₁,y₂…y_k代入q(y)＝||By-g||²，将使q(y)取得最小值的y_j作为全局最优解，从中读出移动端的位置向量p的全局最优解，其中1≤j≤k。

下面结合仿真实验，对本发明的技术效果作详细说明。

1.仿真条件：

蜂窝基站的数量为7，待定位的移动端数量为1；7个基站的位置坐标p_i随机产生，且服从范围为[50m,80m]*[50m,80m]上的二维均匀分布，移动端的位置坐标p服从范围为[-20m,-10m]*[-20m,-10m]上的二维均匀分布；基站与移动端的时钟相偏均服从范围[0.5s,1s]上的一维均匀分布；测量误差Δr_i服从均值为零且方差为σ²的高斯分布；本次仿真实验中，σ的取值为10^-3,10^-2,10^-1，1，针对每次σ的不同取值，进行10000次的独立仿真。

2.仿真内容与结果分析：

本次仿真利用均方误差MSE以及克拉米罗下界CRLB两个性能指标对本发明与现有技术中基于迭代的最小二乘定位方法的定位精度进行了一个对比仿真，其效果如图2所示。

参照图2，本发明定位性能明显优于基于迭代的最小二乘定位方法，定位均方误差低于基于迭代的最小二乘定位方法大约2-3dB,且更加接近CRLB界，说明本发明能够对移动端的位置进行有效估计。

Claims (5)

1.一种TDOA条件下最优联合时间同步与定位的定位方法，包括如下步骤：

(1)设定系统参数：

参与移动端定位且位置已知的蜂窝基站数量为m，移动端数量为1，基站与移动端均位于n维坐标系；移动端在外部时钟为0时发送信号给m个基站，选定坐标零点处的第1个基站作为参考基站，相对于外部参考时钟的时钟相偏为τ₁，接收到移动端信号的外部参考时钟为T₁，观测噪声为ΔT₁；其余各个蜂窝基站的位置向量为p_i，各个基站相对于外部参考时钟的时钟相偏为τ_i，各个基站接收到移动端信号的外部参考时钟为T_i，观测噪声为ΔT_i；待定位的移动端的位置向量为p，移动端相对于外部参考时钟的时钟相偏为τ，其中，m≥4，i＝2,3…m，p_i∈Rⁿ，p∈Rⁿ，τ₁∈R，τ_i∈R，τ∈R，T₁∈R，T_i∈R，R为实数域，位置向量p的单位为米，时钟相偏τ的单位为秒；

(2)建立联合时间同步与定位到达时间差TDOA的观测模型：

根据设定的系统参数，将联合时间同步与定位到达时间值TOA的观测模型中所有T_i分别与T₁作差，得到联合时间同步与定位到达时间差TDOA的观测模型：

<mrow> <msub> <mi>T</mi> <mi>i</mi> </msub> <mo>-</mo> <msub> <mi>T</mi> <mn>1</mn> </msub> <mo>=</mo> <mfrac> <mrow> <mo>|</mo> <mo>|</mo> <mi>p</mi> <mo>-</mo> <msub> <mi>p</mi> <mi>i</mi> </msub> <mo>|</mo> <mo>|</mo> <mo>-</mo> <mo>|</mo> <mo>|</mo> <mi>p</mi> <mo>|</mo> <mo>|</mo> </mrow> <mi>c</mi> </mfrac> <mo>+</mo> <msub> <mi>&amp;tau;</mi> <mn>1</mn> </msub> <mo>-</mo> <msub> <mi>&amp;tau;</mi> <mi>i</mi> </msub> <mo>+</mo> <msub> <mi>&amp;Delta;T</mi> <mi>i</mi> </msub> <mo>-</mo> <msub> <mi>&amp;Delta;T</mi> <mn>1</mn> </msub> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow>

其中，c为信号的传播速度，||p-p_i||为移动端与第i个基站的欧几里得距离，||p||为移动端与参考基站的欧几里得距离；

(3)通过联合时间同步与定位到达时间值TDOA的观测模型，建立移动端位置估计的有约束目标函数：

(3a)对公式(1)两端同乘c，得到关于距离差的观测模型：

r_i,1＝||p-p_i||-||p||+τ_1,i+Δr_i，1 (2)

其中，r_i,1＝cT_i-cT₁，τ_1,i＝cτ₁-cτ_i，Δr_i,1＝cΔT_i-cΔT₁，Δr_i,1为测量误差，其服从均值为零且方差为σ²的高斯分布，且σ²＞0；

(3b)将公式(2)转化为关于移动端位置向量p的线性等式：

对公式(2)两边分别求平方：

2r_i,1||p||+2p_i ^Tp-(p_i ^Tp_i-r_i,1 ²)＝2||p-p_i||Δr_i,1 (3)

定义一个与移动端位置向量p相关的中间变量β，β＝||p||，β≥0，并将其带入(3)式中，得到m-1组关于移动端位置向量p的线性等式：

<mrow> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <mn>2</mn> <msub> <mi>r</mi> <mrow> <mn>2</mn> <mo>,</mo> <mn>1</mn> </mrow> </msub> <mi>&amp;beta;</mi> <mo>+</mo> <mn>2</mn> <msup> <msub> <mi>p</mi> <mn>2</mn> </msub> <mi>T</mi> </msup> <mi>p</mi> <mo>-</mo> <mrow> <mo>(</mo> <msup> <msub> <mi>p</mi> <mn>2</mn> </msub> <mi>T</mi> </msup> <msub> <mi>p</mi> <mn>2</mn> </msub> <mo>-</mo> <msup> <msub> <mi>r</mi> <mrow> <mn>2</mn> <mo>,</mo> <mn>1</mn> </mrow> </msub> <mn>2</mn> </msup> <mo>)</mo> </mrow> <mo>=</mo> <mn>2</mn> <mo>|</mo> <mo>|</mo> <mi>p</mi> <mo>-</mo> <msub> <mi>p</mi> <mn>2</mn> </msub> <mo>|</mo> <mo>|</mo> <msub> <mi>&amp;Delta;r</mi> <mrow> <mn>2</mn> <mo>,</mo> <mn>1</mn> </mrow> </msub> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mn>2</mn> <msub> <mi>r</mi> <mrow> <mn>3</mn> <mo>,</mo> <mn>1</mn> </mrow> </msub> <mi>&amp;beta;</mi> <mo>+</mo> <mn>2</mn> <msup> <msub> <mi>p</mi> <mn>3</mn> </msub> <mi>T</mi> </msup> <mi>p</mi> <mo>-</mo> <mrow> <mo>(</mo> <msup> <msub> <mi>p</mi> <mn>3</mn> </msub> <mi>T</mi> </msup> <msub> <mi>p</mi> <mn>3</mn> </msub> <mo>-</mo> <msup> <msub> <mi>r</mi> <mrow> <mn>3</mn> <mo>,</mo> <mn>1</mn> </mrow> </msub> <mn>2</mn> </msup> <mo>)</mo> </mrow> <mo>=</mo> <mn>2</mn> <mo>|</mo> <mo>|</mo> <mi>p</mi> <mo>-</mo> <msub> <mi>p</mi> <mn>3</mn> </msub> <mo>|</mo> <mo>|</mo> <msub> <mi>&amp;Delta;r</mi> <mrow> <mn>3</mn> <mo>,</mo> <mn>1</mn> </mrow> </msub> </mrow> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <mrow> <mn>2</mn> <msub> <mi>r</mi> <mrow> <mi>m</mi> <mo>,</mo> <mn>1</mn> </mrow> </msub> <mi>&amp;beta;</mi> <mo>+</mo> <mn>2</mn> <msup> <msub> <mi>p</mi> <mi>m</mi> </msub> <mi>T</mi> </msup> <mi>p</mi> <mo>-</mo> <mrow> <mo>(</mo> <msup> <msub> <mi>p</mi> <mi>m</mi> </msub> <mi>T</mi> </msup> <msub> <mi>p</mi> <mi>m</mi> </msub> <mo>-</mo> <msup> <msub> <mi>r</mi> <mrow> <mi>m</mi> <mo>,</mo> <mn>1</mn> </mrow> </msub> <mn>2</mn> </msup> <mo>)</mo> </mrow> <mo>=</mo> <mn>2</mn> <mo>|</mo> <mo>|</mo> <mi>p</mi> <mo>-</mo> <msub> <mi>p</mi> <mi>m</mi> </msub> <mo>|</mo> <mo>|</mo> <msub> <mi>&amp;Delta;r</mi> <mrow> <mi>m</mi> <mo>,</mo> <mn>1</mn> </mrow> </msub> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>4</mn> <mo>)</mo> </mrow> </mrow>

(3c)将m-1组关于移动端位置向量p的线性等式转化为矩阵形式：

By-g＝ε (5)

其中，ε为m-1行的列向量，且A为m-1行n+2列的矩阵，且g为m-1行的列向量，且待估向量y＝[p β]^T，y∈Rⁿ⁺¹；

(3d)根据最小二乘算法，将(5)式转化为无约束目标函数：

<mrow> <munder> <mrow> <mi>m</mi> <mi>i</mi> <mi>n</mi> </mrow> <mrow> <mi>y</mi> <mo>&amp;Element;</mo> <msup> <mi>R</mi> <mrow> <mi>n</mi> <mo>+</mo> <mn>1</mn> </mrow> </msup> </mrow> </munder> <mo>{</mo> <mo>|</mo> <mo>|</mo> <mi>B</mi> <mi>y</mi> <mo>-</mo> <mi>g</mi> <mo>|</mo> <msup> <mo>|</mo> <mn>2</mn> </msup> <mo>}</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>6</mn> <mo>)</mo> </mrow> </mrow>

(3e)为公式(6)添加两个约束条件得到移动端位置估计的双约束目标函数：

<mrow> <munder> <mrow> <mi>m</mi> <mi>i</mi> <mi>n</mi> </mrow> <mrow> <mi>y</mi> <mo>&amp;Element;</mo> <msup> <mi>R</mi> <mrow> <mi>n</mi> <mo>+</mo> <mn>1</mn> </mrow> </msup> </mrow> </munder> <mo>{</mo> <mo>|</mo> <mo>|</mo> <mi>B</mi> <mi>y</mi> <mo>-</mo> <mi>g</mi> <mo>|</mo> <msup> <mo>|</mo> <mn>2</mn> </msup> <mo>:</mo> <msup> <mi>y</mi> <mi>T</mi> </msup> <mi>C</mi> <mi>y</mi> <mo>=</mo> <mn>0</mn> <mo>,</mo> <msub> <mi>y</mi> <mrow> <mi>n</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>&amp;GreaterEqual;</mo> <mn>0</mn> <mo>}</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>7</mn> <mo>)</mo> </mrow> </mrow>

其中，I_n为n行n列的单位矩阵，0_n×1为n行1列的零矩阵，0_1×n为1行n列的行向量，y_n+1为y向量的第n+1个元素；

(4)求解移动端定位估计的双约束目标函数，得到待估向量y的全局最优解：

(4a)获取待估向量y的解算条件式：

设c(y)，q(y)均为y的函数，且c(y)＝y^TCy，q(y)＝||By-g||²，并将(7)式转化为：

<mrow> <munder> <mrow> <mi>m</mi> <mi>i</mi> <mi>n</mi> </mrow> <mrow> <mi>y</mi> <mo>&amp;Element;</mo> <msup> <mi>R</mi> <mrow> <mi>n</mi> <mo>+</mo> <mn>2</mn> </mrow> </msup> </mrow> </munder> <mo>{</mo> <mi>q</mi> <mrow> <mo>(</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>:</mo> <mi>c</mi> <mrow> <mo>(</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>=</mo> <mn>0</mn> <mo>}</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>8</mn> <mo>)</mo> </mrow> </mrow>

令λ满足以下三个条件，λ∈R：

<mrow> <mo>&amp;dtri;</mo> <mi>q</mi> <mrow> <mo>(</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>+</mo> <mi>&amp;lambda;</mi> <mo>&amp;dtri;</mo> <mi>c</mi> <mrow> <mo>(</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>=</mo> <mn>0</mn> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>9</mn> <mo>)</mo> </mrow> </mrow>

c(y)＝0 (10)

<mrow> <msup> <mo>&amp;dtri;</mo> <mn>2</mn> </msup> <mi>q</mi> <mrow> <mo>(</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>+</mo> <mi>&amp;lambda;</mi> <msup> <mo>&amp;dtri;</mo> <mn>2</mn> </msup> <mi>c</mi> <mrow> <mo>(</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>&gt;</mo> <mn>0</mn> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>11</mn> <mo>)</mo> </mrow> </mrow>

其中符号表示求一阶导数，符号表示求两阶导数；

将c(y)和q(y)带入公式(9)、(10)和(11)中，得到待估向量y的解算条件式：

(B^TB+λC)y＝B^Tg (12)

y^TCy＝0 (13)

其中，表示B^TB+λC为正定矩阵；

(4b)根据待估向量y的解算条件式，将移动端位置估计的双约束目标函数的求解问题转换为多项式零点的求解问题：

根据式(12)，计算待估向量y：y＝(B^TB+λC)^-1B^Tg (15)

设是关于λ的函数，令：

将(15)式带入公式(16)得：

将公式(16)代入(13)，得方程实现将移动端位置估计的双约束目标函数的求解问题转化为多项式零点的求解问题；

(4c)利用式(14)计算的定义域I₁，并证明在定义域I₁内严格单调递减；

(4d)根据函数在定义域I₁上严格单调递减，利用二分法在定义域I₁内对λ进行有限次的迭代，得到的零点值；

(4e)将的零点值代入式(15)，得到待估向量y的估计值,并从中读出y_n+1；

(4f)判断y_n+1≥0是否成立，若是，则y就为全局最优解，从中读出移动端的位置向量p的全局最优解；否则执行步骤(4g)；

(4g)将待估向量y的解算条件式(14)更换为：B^TB+λC最多有一个负特征值，并在该条件下计算定义域I₂；

(4h)在定义域I₂内求出方程的所有实数根λ₁,λ₂…λ_k，其中1＜k≤n；

(4i)将λ₁,λ₂…λ_k代入y＝(B^TB+λC)^-1B^Tg，解出待估向量y的多个估计值y₁,y₂…y_k，并将y₁,y₂…y_k代入q(y)＝||By-g||²，将使q(y)取得最小值的y_j作为全局最优解，从中读出移动端的位置向量p的全局最优解，其中1≤j≤k。

2.根据权利要求1所述的一种TDOA条件下最优联合时间同步与定位的定位方法，其特征在于：步骤(2)中所述的将联合时间同步与定位到达时间值TOA的观测模型中所有T_i分别与T₁作差，其各个基站接收到移动端信号的外部参考时钟的表达式为：

<mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <msub> <mi>T</mi> <mn>1</mn> </msub> <mo>=</mo> <mi>&amp;tau;</mi> <mo>-</mo> <msub> <mi>&amp;tau;</mi> <mn>1</mn> </msub> <mo>+</mo> <mfrac> <mrow> <mo>|</mo> <mo>|</mo> <mi>p</mi> <mo>|</mo> <mo>|</mo> </mrow> <mi>c</mi> </mfrac> <mo>+</mo> <msub> <mi>&amp;Delta;T</mi> <mn>1</mn> </msub> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>T</mi> <mn>2</mn> </msub> <mo>=</mo> <mi>&amp;tau;</mi> <mo>-</mo> <msub> <mi>&amp;tau;</mi> <mn>2</mn> </msub> <mo>+</mo> <mfrac> <mrow> <mo>|</mo> <mo>|</mo> <mi>p</mi> <mo>-</mo> <msub> <mi>p</mi> <mn>2</mn> </msub> <mo>|</mo> <mo>|</mo> </mrow> <mi>c</mi> </mfrac> <mo>+</mo> <msub> <mi>&amp;Delta;T</mi> <mn>2</mn> </msub> </mrow> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>T</mi> <mi>m</mi> </msub> <mo>=</mo> <mi>&amp;tau;</mi> <mo>-</mo> <msub> <mi>&amp;tau;</mi> <mi>m</mi> </msub> <mo>+</mo> <mfrac> <mrow> <mo>|</mo> <mo>|</mo> <mi>p</mi> <mo>-</mo> <msub> <mi>p</mi> <mi>m</mi> </msub> <mo>|</mo> <mo>|</mo> </mrow> <mi>c</mi> </mfrac> <mo>+</mo> <msub> <mi>&amp;Delta;T</mi> <mi>m</mi> </msub> </mrow> </mtd> </mtr> </mtable> </mfenced>

3.根据权利要求1所述的一种TDOA条件下最优联合时间同步与定位的定位方法，其特征在于：步骤(4c)中所述的计算的定义域I₁，并证明在定义域I₁内严格单调递减，实现步骤为：

计算的定义域I₁：

(4c1)令矩阵得到矩阵B^TB+λC与矩阵合同，其中I为n+1行n+1列的单位矩阵；

(4c2)将矩阵表示为(C,B^TB)；

(4c3)令得到矩阵与矩阵C合同；

(4c4)利用线性代数中矩阵与矩阵C合同定义，矩阵与矩阵C拥有相同的正负惯性，即矩阵与矩阵C拥有相同的正特征值个数，负特征值个数和零特征值个数，根据得到矩阵有一个负特征值，n个正特征值，即：

λ₁(C,B^TB)＞λ₂(C,B^TB)＞…＞λ_n(C,B^TB)＞0＞λ_n+1(C,B^TB)

其中，λ_j(C,B^TB)表示矩阵的第j个特征值，j＝1,2…n+2，λ_j∈R，,且λ_j降序排列；

(4c5)再次利用线性代数中矩阵B^TB+λC与矩阵I+λ(C，B^TB)合同定义，即矩阵B^TB+λC与矩阵I+λ(C，B^TB)拥有相同的正特征值个数，负特征值个数和零特征值个数，根据公式(14)，得到矩阵I+λ(C,B^TB)的特征值均为正数，即：

1+λλ_j(C,B^TB)＞0

其中，λ_j(C,B^TB)表示矩阵的第j个特征值，j＝1,2…n+2，λ_j∈R，,且λ_j降序排列；

(4c6)根据λ₁(C,B^TB)＞λ₂(C,B^TB)＞…＞λ_n(C,B^TB)＞0＞λ_n+1(C,B^TB)，j＝1,2…n+1，解不等式1+λλ_j(C,B^TB)＞0的步骤如下：

当j∈[1,n]，

<mrow> <msub> <mi>&amp;lambda;</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mi>C</mi> <mo>,</mo> <msup> <mi>B</mi> <mi>T</mi> </msup> <mi>B</mi> <mo>)</mo> </mrow> <mo>&gt;</mo> <mn>0</mn> <mo>&amp;DoubleRightArrow;</mo> <mi>&amp;lambda;</mi> <mo>&gt;</mo> <mo>-</mo> <mfrac> <mn>1</mn> <mrow> <msub> <mi>&amp;lambda;</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mi>C</mi> <mo>,</mo> <msup> <mi>B</mi> <mi>T</mi> </msup> <mi>B</mi> <mo>)</mo> </mrow> </mrow> </mfrac> <mo>&amp;DoubleRightArrow;</mo> <mi>&amp;lambda;</mi> <mo>&gt;</mo> <mo>-</mo> <mfrac> <mn>1</mn> <mrow> <msub> <mi>&amp;lambda;</mi> <mn>1</mn> </msub> <mrow> <mo>(</mo> <mi>C</mi> <mo>,</mo> <msup> <mi>B</mi> <mi>T</mi> </msup> <mi>B</mi> <mo>)</mo> </mrow> </mrow> </mfrac> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>17</mn> <mo>)</mo> </mrow> </mrow>

当j＝n+1，

<mrow> <msub> <mi>&amp;lambda;</mi> <mrow> <mi>n</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mrow> <mo>(</mo> <mi>C</mi> <mo>,</mo> <msup> <mi>B</mi> <mi>T</mi> </msup> <mi>B</mi> <mo>)</mo> </mrow> <mo>&lt;</mo> <mn>0</mn> <mo>&amp;DoubleRightArrow;</mo> <mi>&amp;lambda;</mi> <mo>&lt;</mo> <mo>-</mo> <mfrac> <mn>1</mn> <mrow> <msub> <mi>&amp;lambda;</mi> <mrow> <mi>n</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mrow> <mo>(</mo> <mi>C</mi> <mo>,</mo> <msup> <mi>B</mi> <mi>T</mi> </msup> <mi>B</mi> <mo>)</mo> </mrow> </mrow> </mfrac> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>18</mn> <mo>)</mo> </mrow> </mrow>

取不等式(17)和(18)的交集得λ的定义域I₁：

<mrow> <msub> <mi>I</mi> <mn>1</mn> </msub> <mo>=</mo> <mrow> <mo>(</mo> <mo>-</mo> <mfrac> <mn>1</mn> <mrow> <msub> <mi>&amp;lambda;</mi> <mn>1</mn> </msub> <mrow> <mo>(</mo> <mi>C</mi> <mo>,</mo> <msup> <mi>B</mi> <mi>T</mi> </msup> <mi>B</mi> <mo>)</mo> </mrow> </mrow> </mfrac> <mo>,</mo> <mo>-</mo> <mfrac> <mn>1</mn> <mrow> <msub> <mi>&amp;lambda;</mi> <mrow> <mi>n</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mrow> <mo>(</mo> <mi>C</mi> <mo>,</mo> <msup> <mi>B</mi> <mi>T</mi> </msup> <mi>B</mi> <mo>)</mo> </mrow> </mrow> </mfrac> <mo>)</mo> </mrow> </mrow>

证明在定义域I₁内严格单调递减：

(4c7)对公式(12)两端同时对λ求一阶导得：

<mrow> <mo>(</mo> <msup> <mi>B</mi> <mi>T</mi> </msup> <mi>B</mi> <mo>+</mo> <mi>&amp;lambda;</mi> <mi>C</mi> <mo>)</mo> <mo>(</mo> <mfrac> <mrow> <mi>d</mi> <mi>y</mi> </mrow> <mrow> <mi>d</mi> <mi>&amp;lambda;</mi> </mrow> </mfrac> <mo>)</mo> <mo>=</mo> <mo>-</mo> <mi>C</mi> <mi>y</mi> </mrow>

(4c8)根据公式(9)和(12)，得：

<mrow> <mo>(</mo> <msup> <mi>B</mi> <mi>T</mi> </msup> <mi>B</mi> <mo>+</mo> <mi>&amp;lambda;</mi> <mi>C</mi> <mo>)</mo> <mo>(</mo> <mfrac> <mrow> <mi>d</mi> <mi>y</mi> </mrow> <mrow> <mi>d</mi> <mi>&amp;lambda;</mi> </mrow> </mfrac> <mo>)</mo> <mo>=</mo> <mo>-</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> <mfrac> <mrow> <mi>d</mi> <mi>c</mi> <mrow> <mo>(</mo> <mi>y</mi> <mo>)</mo> </mrow> </mrow> <mrow> <mi>d</mi> <mi>y</mi> </mrow> </mfrac> </mrow>

对上式进行移项得：

(4c9)根据公式(16)，将的两端同时对λ求一阶导：

(4c10)根据矩阵B^TB+λC为正定矩阵，利用正定矩阵的定义得即在定义域I₁内严格单调递减。

4.根据权利要求1所述的一种TDOA条件下最优联合时间同步与定位的定位方法，其特征在于：步骤(4g)中所述的将待估向量y的解算条件式(14)更换为B^TB+λC最多有一个负特征值，并在该条件下计算定义域I₂，实现步骤为：

(4g1)计算定义域I₂的前四步与计算定义域I₁的前四个步骤完全一致；

(4g2)根据矩阵B^TB+λC最多只有一个负特征值，j＝1,2…n+1，且λ₁(C,B^TB)＞λ₂(C,B^TB)＞…＞λ_n(C,B^TB)＞0＞λ_n+1(C,B^TB)得：

当λ≥0，j∈[1,n]时，1+λλ_j(C,B^TB)＞0恒成立，要使B^TB+λC最多有一个负特征值，需要满足不等式：

1+λλ_n+1(C,B^TB)＞0 (19)

解不等式(19)得：

<mrow> <mi>&amp;lambda;</mi> <mo>&gt;</mo> <mo>-</mo> <mfrac> <mn>1</mn> <mrow> <msub> <mi>&amp;lambda;</mi> <mrow> <mi>n</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mrow> <mo>(</mo> <mi>C</mi> <mo>,</mo> <msup> <mi>B</mi> <mi>T</mi> </msup> <mi>B</mi> <mo>)</mo> </mrow> </mrow> </mfrac> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>20</mn> <mo>)</mo> </mrow> </mrow>

当λ＜0，j＝n+1时，1+λλ_n+1(C,B^TB)＞0恒成立，当j∈[1,n]时，1+λλ₁(C,B^TB)＜1+λλ₂(C,B^TB)＜…1+λλ_n(C,B^TB)，要使B^TB+λC最多有一个负特征值，需满足1+λλ_j(C,B^TB)中有n-1个正数和一个负数，即需要满足不等式：

1+λλ₁(C,B^TB)＜0＜1+λλ₂(C,B^TB) (21)

解不等式(21)得：

<mrow> <mo>-</mo> <mfrac> <mn>1</mn> <mrow> <msub> <mi>&amp;lambda;</mi> <mn>2</mn> </msub> <mrow> <mo>(</mo> <mi>C</mi> <mo>,</mo> <msup> <mi>B</mi> <mi>T</mi> </msup> <mi>B</mi> <mo>)</mo> </mrow> </mrow> </mfrac> <mo>&lt;</mo> <mi>&amp;lambda;</mi> <mo>&lt;</mo> <mo>-</mo> <mfrac> <mn>1</mn> <mrow> <msub> <mi>&amp;lambda;</mi> <mn>1</mn> </msub> <mrow> <mo>(</mo> <mi>C</mi> <mo>,</mo> <msup> <mi>B</mi> <mi>T</mi> </msup> <mi>B</mi> <mo>)</mo> </mrow> </mrow> </mfrac> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>22</mn> <mo>)</mo> </mrow> </mrow>

取不等式(20)和(22)的交集得λ的定义域I₂：

<mrow> <msub> <mi>I</mi> <mn>2</mn> </msub> <mo>=</mo> <mrow> <mo>(</mo> <mo>-</mo> <mfrac> <mn>1</mn> <mrow> <msub> <mi>&amp;lambda;</mi> <mn>2</mn> </msub> <mrow> <mo>(</mo> <mi>C</mi> <mo>,</mo> <msup> <mi>B</mi> <mi>T</mi> </msup> <mi>B</mi> <mo>)</mo> </mrow> </mrow> </mfrac> <mo>,</mo> <mo>-</mo> <mfrac> <mn>1</mn> <mrow> <msub> <mi>&amp;lambda;</mi> <mn>1</mn> </msub> <mrow> <mo>(</mo> <mi>C</mi> <mo>,</mo> <msup> <mi>B</mi> <mi>T</mi> </msup> <mi>B</mi> <mo>)</mo> </mrow> </mrow> </mfrac> <mo>)</mo> </mrow> <mo>&amp;cup;</mo> <mrow> <mo>(</mo> <mo>-</mo> <mfrac> <mn>1</mn> <mrow> <msub> <mi>&amp;lambda;</mi> <mrow> <mi>n</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mrow> <mo>(</mo> <mi>C</mi> <mo>,</mo> <msup> <mi>B</mi> <mi>T</mi> </msup> <mi>B</mi> <mo>)</mo> </mrow> </mrow> </mfrac> <mo>,</mo> <mo>+</mo> <mi>&amp;infin;</mi> <mo>)</mo> </mrow> <mo>.</mo> </mrow>

5.根据权利要求1所述的一种TDOA条件下最优联合时间同步与定位的定位方法，其特征在于：步骤(4h)中所述的在定义域I₂内求出方程的所有实数根λ₁,λ₂…λ_k，其中1＜k≤n，实现步骤为：

(4h1)利用矩阵的相似对角化规则求解出一个正交矩阵P，该矩阵能够将矩阵B^TB和矩阵C同时对角化：

P^T(B^TB)P＝diag(γ₁,γ₂…γ_n+1)

P^TCP＝diag(δ₁,δ₂…δ_n+1)

其中P为n+1行n+1列的正交矩阵，diag(…)表示n+1行n+1列的对角矩阵，γ₁,γ₂…γ_n+1，δ₁,δ₂…δ_n+1为对角矩阵的对角线上的元素；

(4h2)将矩阵方程式转换为代数方程式：

<mrow> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>z</mi> <mo>=</mo> <mn>1</mn> </mrow> <mrow> <mi>n</mi> <mo>+</mo> <mn>1</mn> </mrow> </munderover> <mfrac> <mrow> <msup> <msub> <mi>f</mi> <mi>z</mi> </msub> <mn>2</mn> </msup> <msub> <mi>&amp;delta;</mi> <mi>z</mi> </msub> </mrow> <msup> <mrow> <mo>(</mo> <msub> <mi>&amp;gamma;</mi> <mi>z</mi> </msub> <mo>+</mo> <msub> <mi>&amp;lambda;&amp;delta;</mi> <mi>z</mi> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mfrac> <mo>=</mo> <mn>0</mn> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>23</mn> <mo>)</mo> </mrow> </mrow>

其中f为n+1行的列向量，f＝P^TB^Tg，f_z为列向量f的第z个元素，1≤z≤n+1；

(4h3)根据公式(23)，利用高阶多项式方程求根法求出方程的所有实数根λ₁,λ₂…λ_k，其中1＜k≤n。