CN102736090B - 位置计算方法和位置计算装置 - Google Patents

Abstract

本发明提供了一种新的位置计算方法和位置计算装置。该位置计算方法包括：接收作为定位信号中的一种的GPS卫星信号。并且，进行位置计算运算，以计算GPS接收器的位置。该位置计算运算基于至少将GPS接收器的位置的变化作为随机变量的概率分布模型即状态变化正态分布模型的运算。

Description

位置计算方法和位置计算装置

技术领域

本发明涉及位置计算方法和位置计算装置。

背景技术

作为利用定位信号的定位系统，公知有GPS(全球定位系统)，被用在内置于移动电话或汽车导航装置等的位置计算装置中。GPS基于多个GPS卫星的位置或从各GPS卫星到位置计算装置的伪距等信息进行位置计算运算以得到位置计算装置的位置坐标和时钟误差。

在利用卫星定位系统的位置计算运算中，众所周知的是利用所谓的最小二乘法的位置计算运算，即，使利用关于多个定位卫星的接收信号测定的伪距中可能包含的误差(以下称为“伪距误差”)的平方和最小化(例如专利文献1)。

现有技术文献

专利文献

专利文献1：JP特开第2009-97897号

在现有的位置计算运算中，通常是在每个进行位置计算运算的时刻(以下称为“运算时刻”)利用最小二乘法等一次次地单独进行位置计算运算。因此，例如，如果在测定的伪距中瞬间混入了很大的误差，缺点是此时的位置计算运算的精度急剧下降，而优点是对之后的位置计算运算没有影响。虽然是瞬间的，但作为运算结果得到的是与真实位置偏离很多的位置，即产生所谓的位置跳跃，有时会造成很大的问题。

发明内容

本发明鉴于上述问题而做出，其目的是提供一种用于提高位置计算的准确性的新的位置计算方法。

解决上述问题的第一方式是一种位置计算方法，包括：接收定位信号；以及，利用所接收的定位信号进行位置计算运算，该位置计算运算是基于至少将过去的运算结果的变化作为随机变量的给定的概率分布模型的运算。

另外，作为其他方式，也可以构造位置计算装置，该位置计算装置具有：接收部，接收定位信号；以及位置计算部，利用接收部接收的定位信号进行位置计算运算，该位置计算运算是基于至少将过去的运算结果的变化作为随机变量的给定的概率分布模型的运算。

根据该第一方式等，利用所接收的定位信号进行位置计算运算，该位置计算运算是基于至少将过去的运算结果的变化作为随机变量的给定的概率分布模型的运算。通过这样，可以不使各运算结果独立，而使运算结果具有时间上的相关性而进行位置计算。因此，即使测定量瞬间混入很大的误差，也可以防止位置计算精度降低。

另外，作为第二方式，也可以构造以下位置计算方法：在第一方式的位置计算方法中，还包括：基于上述过去的运算结果，对上述变化的标准进行更新；以及，利用上述标准计算上述变化。

根据该第二方式，基于过去的运算结果对运算结果的变化的标准进行更新，从而可以使这一次的位置计算运算结果反映过去的位置计算运算的结果。利用更新后的标准来计算运算结果的变化，从而提高位置计算的准确性。

另外，作为第三方式，也可以构造以下位置计算方法：在第一或第二方式的位置计算方法中，上述位置计算运算是基于至少将上述过去的运算结果作为随机变量的正态分布模型的运算，该位置计算方法还包括基于上述过去的运算结果对上述正态分布模型进行更新。

根据该第三方式，进行基于至少将过去的运算结果作为随机变量的正态分布模型的位置计算运算。即，假定运算结果的变化分布按照正态分布。而且，不是固定地设定正态分布模型，而是基于过去的位置计算运算的结果对正态分布模型进行更新，从而可以基于与时间变化对应的适当的正态分布模型进行位置计算。

另外，作为第四方式，也可以构造以下位置计算方法：在第三方式的位置计算方法中，还包括基于所接收的信号测量伪距，上述位置计算运算是基于将过去的运算结果的变化作为随机变量的正态分布模型和将伪距中包含的误差作为随机变量的概率分布模型的运算。

根据该第四方式，基于所接收的定位信号测量伪距。另外，进行基于将过去的运算结果的变化作为随机变量的正态分布模型和将伪距中包含的误差作为随机变量的概率分布模型的位置计算运算。除了将运算结果的变化作为随机变量的正态分布模型以外，还使用将伪距中包含的误差作为随机变量的概率分布模型，从而可以进行适当的位置计算运算。

更具体的是，作为第五方式，可构造以下位置计算方法：将第四方式的位置计算方法中的将上述伪距中包含的误差作为随机变数的概率分布模型形成为正态分布模型或正态混合分布模型。

在通常的定位环境下，可以假定伪距中包含的误差分布遵从正态分布。但是，如果像多径环境那样，假定伪距中包含的误差分布是正态分布有时未必适当。因此，如第五方式那样，将伪距中包含的误差作为随机变数的概率分布模型作为正态分布模型或正态混合分布模型进行位置计算运算，从而可以实现与定位环境相符的适当的位置计算。

另外，作为第六方式，也可以构造以下位置计算方法：在第三至第五方式中的任一个的位置计算方法中，上述位置计算运算是计算位置信息、钟差、速度信息以及时钟漂移中的至少任意一个的运算，上述正态分布模型是将过去的运算结果的各元素的变化作为随机变量的模型。

根据该第六方式，除了位置信息以外，还进行计算钟差、速度信息和时钟漂移中的任意个的位置计算运算。此时，使用将运算结果的各要素的变化作为随机变量的正态分布模型，从而可以使各要素具有时间相关性地进行位置计算。

另外，作为第七方式，也可以构造以下位置计算方法：在第一至第六方式中的任一个的位置计算方法中，还包括：通过利用所接收的信号进行预定的收敛运算来计算计算参数中的预定参数的值，该计算参数能够通过使用上述位置计算运算来计算，该位置计算运算是利用过去通过使用收敛运算而算出的该预定参数的值进行的运算。

根据该第七方式，使用接收信号进行预定的收敛运算，计算能够通过位置计算运算来计算的计算参数中的预定参数的值。并且，使用过去通过收敛运算算出的预定参数的值来进行位置计算运算。通过这样，就特定的计算参数而言，不使其具有时间相关性就可以进行预定的收敛运算、计算其值。并且，将该计算结果反映在基于给定的随机分布模型而确定的位置计算运算中，从而可以有效地进行位置计算。

附图说明

图1是状态变化正态分布模型的说明图。

图2是伪距误差正态分布模型的说明图。

图3是多径说明图。

图4是示出自相关的一个例子的图。

图5是示出多径环境中的自相关的一个例子的图。

图6是示出多径环境中的自相关的一个例子的图。

图7是第一伪距误差正态混合分布模型的说明图。

图8是第二伪距误差正态混合分布模型的说明图。

图9是示出位置计算结果的一个例子的图。

图10是示出位置计算结果的一个例子的图。

图11是示出位置计算结果的一个例子的图。

图12是示出移动电话的功能结构的一个例子的框图。

图13是示出基带处理电路部的电路结构的一个例子的图。

图14是示出基带处理的流程的流程图。

图15是示出第一位置计算处理的流程的流程图。

图16是示出第二位置计算处理的流程的流程图。

图17是示出第三位置计算处理的流程的流程图。

图18是示出第四位置计算处理的流程的流程图。

图19是示出多普勒定位处理的流程的流程图。

具体实施方式

以下参照附图就本发明适当的实施方式的一个例子进行说明。本实施方式是应用于一种卫星定位系统即GPS(全球定位系统)的实施方式。另外，可以使用本发明的实施方式当然不限于以下说明的实施方式。

1.原理

就本实施方式中的位置计算方法(位置计算原理)进行说明。在利用GPS的卫星定位系统中，GPS卫星(其为一种定位卫星)将包含年历和星历等卫星轨道数据的导航数据置于GPS卫星信号(其为一种定位卫星信号)中进行发送。

GPS卫星信号是通过C/A(Coarse and Acquisition)码(其为一种扩频码)、通过被称为扩频方式的CDMA(Code Division MultipleAccess)方式进行调制的1.57542[GHz]的通信信号。C/A码是使码长为1023个码片作为1PN帧的重复周期为1ms的伪随机噪声码，是每个GPS卫星固有的码。

GPS利用多个GPS卫星的位置和从每个GPS卫星到GPS接收器(位置计算装置)的伪距等与GPS卫星信号的接收信号相关的测定量，进行计算GPS接收器的位置(位置坐标)和时钟误差(钟差)的位置计算运算。伪距是进行卫星搜索、捕获GPS卫星信号、利用作为结果得到的测量信息进行测量而获得的。

卫星搜索包括对接收GPS卫星信号而得到的信号进行的相位方向的相关运算(以下称为“相位搜索”)和频率方向的相关运算(以下称为“频率搜索”)。通过进行相位搜索，可以获取GPS接收器在接收GPS卫星信号时的C/A码的相位(以下称为“码相位”)作为测量信息。另外，通过进行频率搜索，可以获取GPS接收器接收的GPS卫星信号的接收频率和多普勒频率作为测量信息。

在观念上，在GPS卫星与GPS接收器之间可以排列多个C/A码。GPS卫星与GPS接收器之间的距离不限于是C/A码的整倍数的长度，也可以产生小数部分。与该伪距的小数部分对应的是码相位。伪距的整数部分通过GPS接收器和GPS卫星的大致位置得到。

在本实施方式中，将GPS接收器(位置计算装置)的位置矢量“P”和钟差(时钟误差)“b”作为未知数，将以这些为分量的GPS接收器的状态矢量“X”定义为下式(1)。

[式1]

$X=\left[\begin{array}{c}P\\ b\end{array}\right]=\left[\begin{array}{c}x\\ y\\ z\\ b\end{array}\right]=\left[\begin{array}{c}{X}_1\\ X_2\\ X_3\\ X_4\end{array}\right]...\left(1\right)$

在式(1)中，“x、y、z”是进行定位的坐标系中的GPS接收器(位置计算装置)的位置矢量“P”的各轴的位置分量。为了后续方便，将位置分量“x、y、z”和钟差“b”在作为状态矢量“X”的分量这个意义上标记为“X_i＝(x，y，z，b)^T＝(X₁，X₂，X₃，X₄)^T”。下标的“i＝1、2、3、4”表示状态矢量“X”的分量编号。上标的“T”表示矩阵的转置。另外，在本实施方式的数学式中省略了表示矢量的箭头标记，以简要地记载。

设GPS接收器的真正的位置矢量为“P＝P₀+δP”，真正的钟差为“b＝b₀+δb”。“δP”和“δb”是用于各个位置矢量和钟差的初始值的未知的变化量(校正量)，分别称为“位置变化量矢量”和“钟差变化量”。另外，汇总“位置变化量矢量δP”和“钟差变化量δb”，将状态变化量矢量“δX”定义为下式(2)。

[式2]

$\mathrm{\δX}=\left[\begin{array}{c}\mathrm{\δP}=P-P_0\\ \mathrm{\δb}=b-b_0\end{array}\right]=\left[\begin{array}{c}\mathrm{\δx}=x-x_0\\ \mathrm{\δy}=y-y_0\\ \mathrm{\δz}=z-z_0\\ \mathrm{\δb}=b-b_0\end{array}\right]=\left[\begin{array}{c}\mathrm{\δ}{X}_1\\ {\mathrm{\δX}}_2\\ {\mathrm{\δX}}_3\\ {\mathrm{\δX}}_4\end{array}\right]...\left(2\right)$

在现已使用的GPS的定位计算中，例如利用下式(3)中的观测方程式来计算状态变化量矢量“δX”。

[式3]

δρ＝GδX+E…(3)

在本说明书中，将在GPS卫星中GPS接收器捕获成功的卫星(以下称为“捕获卫星”)的编号标记为“k”，原则上用下标来标记各个量。另外，在本实施方式中，就GPS接收器对总共K个GPS卫星捕获成功为例进行说明。即就“k＝1、2、……、K”进行说明。

式(3)的左边的“δρ”是伪距校正量矢量，表示为下式(4)。

[式4]

$\mathrm{\δ\ρ}=\left[\begin{array}{c}{\mathrm{\δ\ρ}}_1\\ {\mathrm{\δ\ρ}}_2\\ .\\ .\\ .\\ {\mathrm{\δ\ρ}}_K\end{array}\right]...\left(4\right)$

即，伪距校正量矢量“δρ”是以各捕获卫星的伪距校正量矢量“δρ_k”为分量的矢量。

另外，式(3)的右边的矩阵“G”是确定捕获卫星的卫星配置的几何矩阵，表示为下式(5)。

[式5]

$G=\left[\begin{array}{cc}{\left({-l}_1\right)}^T& 1\\ {\left({-l}_2\right)}^T& 1\\ .& .\\ .& .\\ .& .\\ {\left({-l}_K\right)}^T& 1\end{array}\right]...\left(5\right)$

其中，“l_k”是从GPS接收器向着第k个捕获卫星的视线方向的矢量(以下称为“视线矢量”)。

另外，式(3)的右边的“E”是伪距误差矢量，表示为下式(6)。

[式6]

$E=\left[\begin{array}{c}{\mathrm{\ϵ}}_1\\ {\mathrm{\ϵ}}_2\\ .\\ .\\ .\\ {\mathrm{\ϵ}}_K\end{array}\right]...\left(6\right)$

其中，“ε_k”是从第k个捕获卫星到GPS接收器的伪距中包含的误差(以下称为“伪距误差”)。伪距误差矢量“E”是以各捕获卫星的伪距误差“ε_k”为分量的矢量。

式(3)中的未知数一共四个，即位置变化量矢量“δP＝(δx，δy，δz)”和钟差变化量“δb”。因此，在K≥4的情况下可以解出式(3)。如果K＞4，则式(3)是超定方程。这种情况下，一般使用最小二乘法，即求K个捕获卫星的伪距误差“ε”的平方和为最小的解。具体例如通过下式(7)求出式(3)的最小二乘解。

[式7]

δX＝(G^TG)^-1G^Tδρ…(7)

在利用式(3)的观测方程式的位置计算运算中，在每个运算时刻接收GPS卫星信号、测量伪距，解析求出式(7)的最小二乘解。这种情况下，在每个运算时刻独立地计算状态变化量矢量“δX”。

但是，GPS接收器(位置计算装置)通常连续地移动。例如，如果携带着GPS接收器的使用者正常行走或正常行驶，则每个运算时刻的GPS接收器的移动距离应该不会有很大的差。即，只要使用者不进行猛烈的动作，就可以将GPS接收器的位置变化量矢量“δP”的变化几乎当作是固定的。将GPS接收器设置在汽车等移动体上的情况下也是一样。因此，本发明的发明者想到定义将位置计算运算的运算结果的变化作为随机变量的给定的概率分布模型。

例如，使用者进行正常行走，如果按照“一秒钟的间隔”进行位置计算运算，则每秒的GPS接收器的移动距离为“几米左右”，即使考虑到个体差异，该误差幅度也是“几十公分左右”。因此，假定步行的情况下，在位置计算运算的运算结果中对位置变化量矢量“δP”应该可以假设期望值为“几米左右”，标准差为“几十公分左右”。另外，如果考虑将GPS接收器设置在汽车上，则位置变化量矢量“δP”应该可以假设期望值为“几十米左右”，标准差为“几米左右”。

也同样考虑钟差“b”，如果按短时间间隔考虑，则很难想象钟差“b”有多大的变化。因此，对于钟差变化量“δb”也可以假定概率分布模型。

因此，本发明的发明者对于以位置矢量“P”和钟差“b”为分量的状态矢量“X”，假定了将其变化即状态变化量矢量“δX”为随机变量的正态分布型的概率分布模型(以下称为“状态变化正态分布模型”)。

图1是状态变化正态分布模型(状态变化GM(高斯模型))的说明图。图1是表示着眼于状态变化量矢量“δX”的一个分量“δX_i”的模型。横轴是状态变化量的分量“δX_i”，纵轴是其概率密度“p(δX_i)”。将状态变化量矢量“δX_i”的期望值“μ”称为“状态变化量期望值”，将状态变化量矢量“δX_i”的标准差“σ”称为“状态变化量标准差”。

状态变化量矢量“δX_i”作为随机变量，假定其分布按照正态分布N(μ，σ²)。但是，在本实施方式中，状态变化量矢量“δX”包括位置变化量矢量“P”的各分量“(δx，δy，δz)”和钟差变化量“δb”。因此，扩展一维正态分布模型，作为下式(8)中提供的多维状态变化正态分布模型“p_δX ^(m)(δX)”使用。

[式8]

$p_{\mathrm{\δX}}^{\left(m\right)}\left(\mathrm{\δX}\right)=\frac{1}{{\left(\sqrt{2\mathrm{\π}}\right)}^m\sqrt{\left|S_{\mathrm{\δX}}\right|}}\exp (-\frac{1}{2}{(\mathrm{\δX}-{\mathrm{\μ}}_{\mathrm{\δX}})}^T{S}_{\mathrm{\δX}}^{-1}(\mathrm{\δX}-{\mathrm{\μ}}_{\mathrm{\δX}}))...\left(8\right)$

在式(8)中，上标的括弧“m”表示作为随机变量处理的计算参数的数量。计算参数是通过位置计算运算可计算的参数，在本实施方式中，“(δx，δy，δz，δb)”这四个参数相当于此。如果将四个计算参数都作为随机变量，则“m＝4”。也可以将四个计算参数中的一部分计算参数作为随机变量，具体在后面进行说明。这种情况下，作为随机变量的计算参数的数量为“m”。

式(8)的多维状态变化正态分布模型“p_δX ^(m)(δX)”是提供非独立的m个随机变量的同时发生概率的概率密度函数。“μ_δX”是状态变化量矢量“δX”的各分量的期望值即状态变化量期望值矢量。另外，“S_δX”是状态变化量矢量“δX”的分量之间的协方差即状态变化量协方差矩阵。下面将利用数学式就这些进行具体说明。

在本实施方式中，使用上述的多维状态变化正态分布模型“p_δX ^(m)(δX)”，利用最大似然估计法，求出状态变化量矢量“δX”的最大似然估计值。这是接收定位信号并利用定位信号的接收信号的运算，相当于基于将位置计算运算的运算结果的变化作为随机变量的给定的概率分布模型而确定的位置计算运算。如果求出状态变化量矢量“δX”的最大似然估计值，通过将状态矢量“X”的初始值“X₀”加上该最大似然估计值，则可以估计真正的状态矢量“X”。

并且，在本实施方式的位置计算运算中，基于GPS卫星信号的接收信号来测量伪距。从而，进行基于上述状态变化正态分布模型以及将伪距中包含的误差作为随机变量的概率分布模型而确定的位置计算运算。即，伪距误差“ε”也作为随机变量进行处理，假设给定的概率分布模型进行运算。在本实施方式中，将伪距误差“ε”作为随机变量的概率分布模型成为正态分布模型或正态混合分布模型。

图2是伪距误差正态分布模型(伪距误差GM)的说明图。在图2中，横轴是伪距误差“ε”，纵轴是其概率密度“p_ε(ε)”。将伪距误差“ε”的期望值“μ_ε”称为“伪距误差期望值”，将伪距误差“ε”的标准差“σ_ε”称为“伪距误差标准差”。

假定伪距误差“ε”按照正态分布N(μ_ε，(σ_ε)²)，求出观测方程式的解的方法是最小二乘法。广泛使用的是“μ_ε＝0，σ_ε＝1”的伪距误差正态分布模型，此时的最小二乘解与式(7)一致。

上述的伪距误差正态分布模型“p_ε(ε)”在GPS接收器的天空打开的环境(例如天空开放环境)下有用。这是因为在天空开放环境下可以从GPS卫星接收高质量的GPS卫星信号，因此，可以期待伪距误差“ε”接近零。但是，在GPS接收器的天空未打开的环境下，伪距误差“ε”的动作发生变化，其典型的例子是多径环境。

图3是多径环境的说明图。虚线表示直接波信号，点划线表示间接波信号。由于间接波信号的存在，GPS接收器测量的伪距产生误差。伪距虽然是通过GPS接收器利用进行相位搜索而获得的码相位来计算的，但在多径环境下，在该码相位上叠加误差。该码相位的误差相对真正的码相位有正误差和负误差。

图4是码相位检测原理的说明图。在图4中，横轴为码相位，纵轴为相关值，示出了C/A码的自相关值的简要例子。在以下的说明中，提到的相关值是指相关值的大小(绝对值)。

理想的C/A码的自相关值是以峰值为顶点的左右对称的大致三角形的形状。这种情况下，与相关值的峰值(以下称为“相关峰值”)对应的相位是所接收的GPS卫星信号的C/A码的相位。为了检测相关峰值，例如，考虑相对某个码相位超前一定量的码相位(以下称为“超前相位”)和滞后一定量的码相位(以下称为“滞后相位”)。并且，将超前相位的相关值与滞后相位的相关值相等的相关值检测为相关峰值。

在图4中，由于相关值是左右对称的大致三角形状，因此超前相位的相关值与滞后相位的相关值的中心码相位的相关值成为相关峰值，检测对应的码相位作为“峰值相位”。将该检测的峰值相位称为“检测峰值相位”。虽然图4的相关值的形状是理想形状，但在多径环境下相关值的形状是变化的。

图5和图6是示出多径环境下的相关值形状的一个例子的图。图5是间接波信号以与直接波信号相同的相位到达时的相关值曲线的一个例子，图6是间接波信号以与直接波信号相反的相位到达时的相关值曲线的一个例子。在这些图中，示出了分别与直接波信号、间接波信号、合成该直接波信号和间接波信号而得的多径信号对应的相关值曲线。横轴是码相位，纵轴是相关值。

对于间接波信号的相关值虽然与对于直接波信号的相关值一样形成为大致三角形形状，但间接波信号的相关峰值的大小小于直接波信号的相关峰值。这是由于从GPS卫星发送的GPS卫星信号在建筑物或地面上反射或穿过障碍物，从而发出时的信号强度在接收时变弱的缘故。

另外，间接波信号的峰值相位比直接波信号的峰值相位滞后。这是由于从GPS卫星发送的GPS卫星信号在建筑物或地面上反射或绕射障碍物，从而从GPS卫星到GPS接收器的传输距离变得更长的缘故。由于对于多径信号的相关值是直接波信号的相关值与间接波信号的相关值之和，因此三角形状扭曲、不形成以峰值为中心的左右对称。

如果间接波信号以与直接波信号相同的相位到达GPS接收器，直接波信号与间接波信号则彼此加强。因此，合成波信号的相关值是对于直接波信号的相关值的大小与对于间接波信号的相关值的大小的总和。这种情况下，相关值的形状例如是图5所示的形状。在图5中，检测峰值相位是比真正的峰值相位滞后的相位。

与此相反，如果间接波信号与直接波信号相比形成例如半个周期以上且不足一个周期的范围等的相反相位，则直接波信号与间接波信号相互抵消。因此，合成波信号的相关值是对于直接波信号的相关值的大小减去对于间接波信号的相关值的大小而得的值。这种情况下，相关值的形状例如是图6所示的形状。在图6中，检测峰值相位是比真正的峰值相位超前的相位。

在此，将检测峰值相位与真正的峰值相位的相位差定义为“码相位误差”。并且，为了方便起见，将检测峰值相位比真正的峰值相位滞后时的码相位误差的符号定义为“正”，将检测峰值相位比真正的峰值相位超前时的码相位误差的符号定义为“负”。这种情况下，在图5中，码相位误差为“正”，在图6中，码相位误差为“负”。

如上所述，由于使用码相位计算伪距，因此码相位的误差作为误差与伪距叠加。即，如图5所示，如果码相位包含正误差，伪距误差“ε”就形成正误差。而如图6所示，如果码相位包含负误差，伪距误差“ε”就形成负误差。

这样，在多径环境等中，根据码相位误差的正负符号，伪距误差“ε”可以成为正值，也可以成为负值。伪距误差“ε”除了码相位误差，当然也包含其它误差分量。因此，伪距误差“ε”可以根据这些误差分量的大小采用各种值。

因此，本发明的发明者判断，例如在多径环境下，不假定伪距误差“ε”的分布是正态分布，而假定其为混合了多个正态分布的正态混合分布更为确切。作为假定伪距误差“ε”是正态混合分布的概率分布模型，在本实施方式中举例说明了两种模型。

图7是作为伪距误差GMM的一个例子的第一伪距误差正态混合分布模型(第一伪距误差GMM(高斯混合模型))的说明图。第一伪距误差GMM“f(ε)”定义为合成了第一正态分布模型“p₁(ε)”和第二正态分布模型“p₂(ε)”的两种正态分布模型的函数。在图7中，横轴是伪距误差“ε”(单位为米)，纵轴是其概率密度“f(ε)”。另外，用虚线表示第一正态分布模型“p₁(ε)”，用点划线表示第二正态分布模型“p₂(ε)”，用粗实线表示第一伪距误差GMM“f(ε)”。

第一正态分布模型“p₁(ε)”是假设GPS接收器接收直接波信号时的伪距误差“ε”的概率分布模型。如果没有间接波信号的影响，可以期望伪距误差“ε”为“0米”左右的值。因此，根据将伪距误差期望值设为“0米”(μ_ε＝0)、将伪距误差标准差设为比较小的值“10米”(σ_ε＝10)的正态分布函数“N(0，10²)”，将第一正态分布模型“p₁(ε)”进行模型化。

第二正态分布模型“p₂(ε)”是假设GPS接收器接收直接波信号和间接波信号时的伪距误差“ε”的模型函数。如上所述，根据码相位误差的正负符号，伪距误差“ε”可以采用正负的任意值。因此，将第二正态分布模型“p₂(ε)”的伪距误差期望值设为正负的中心即“0米”(μ_ε＝0)。另外，根据码相位误差的大小，伪距误差“ε”可具有非常大的误差。因此，较大地估计伪距误差标准差，设为“100米”(σ_ε＝100)。即，通过正态分布函数“N(0，100²)”将第二正态分布模型“p₂(ε)”模型化。

以如下方式决定第一正态分布模型“p₁(ε)”和第二正态分布模型“p₂(ε)”的混合比例(合成比例)。通常认为GPS接收器受多径影响的比例并不是那么高。因此，假定GPS接收器按一成的概率受到多径的影响，设第一正态分布模型“p₁(ε)”的比例为“9/10”(九成)，设第二正态分布模型“p₂(ε)”的比例为“1/10”(一成)。因此，作为混合(合成)9/10的第一正态分布模型“p₁(ε)”的函数和1/10的第二正态分布模型“p₂(ε)”的函数而得的函数，第一伪距误差GMM“f(ε)”被模型化。

图8是作为伪距误差GMM的其它例子的第二伪距误差正态混合分布模型(第二伪距误差GMM)的说明图。第二伪距误差GMM“f(ε)”与图7的第一伪距误差GMM一样，是混合(合成)了第一正态分布模型“p₁(ε)”和第二正态分布模型“p₂(ε)”两种正态分布模型的概率密度函数。

第一正态分布模型“p₁(ε)”以作为假定GPS接收器接收直接波信号时的伪距误差“ε”的模型函数的正态分布函数“N(0，10²)”进行了模型化。另外，第二正态分布模型“p₂(ε)”以作为假定在多径环境下GPS接收器接收间接波信号时的伪距误差“ε”的模型函数的正态分布函数“N(100，20²)”进行了模型化。

与图7的第一伪距误差GMM“f(ε)”的不同点在于，第二正态分布模型“p₂(ε)”未定义为假定接收直接波信号和间接波信号时的模型函数，而是定义为假定只接收间接波信号时的模型函数。这种情况例如在GPS接收器的周围有高层建筑物等障碍物的环境下有可能发生。由于GPS接收器的天空没有完全打开，因此从GPS卫星发送的信号在障碍物反射，作为间接波信号到达GPS接收器。但是是被障碍物遮挡、直接波信号不直接到达GPS接收器的环境。城市峡谷环境是典型的例子。

本发明的发明者的实验结果揭示了在城市峡谷环境中，伪距误差“ε”有包含“+100～+150米”左右的正误差的趋势。因此，在第二伪距误差GMM“f(ε)”中，将第二正态分布模型函数“p₂(ε)”的伪距误差期望值设为“100米”(μ_ε＝100)，将伪距误差标准差设为“20米”(σ_ε＝20)。即，作为“p₂(ε)＝N(100，20²)”进行了模型化。

例如假定GPS接收器按两成的概率受到多径的影响，来决定第一和第二正态分布模型的合成比例。即，设第一正态分布模型“p₁(ε)”的比例为“8/10”(八成)，设第二正态分布模型“p₂(ε)”的比例为“2/10”(二成)。因此，作为混合(合成)8/10的第一正态分布模型“p₁(ε)”的函数和2/10的第二正态分布模型“p₂(ε)”的函数而得的函数，第二伪距误差GMM“f(ε)”被模型化。

如上所述，作为将伪距误差“ε”作为随机变量的概率分布模型，定义了正态分布模型和正态混合分布模型两种模型。在进行位置计算时，根据定位环境来区分使用伪距误差正态分布模型和伪距误差正态混合分布模型，进行位置计算运算。以下利用数学式就使用每个概率分布模型时的具体位置计算方法进行详细说明。

1-1.如果使用伪距误差正态分布模型(伪距误差GM)

(1)第一位置计算方法

第一位置计算方法是使用伪距误差正态分布模型、将状态变化量矢量“δX”的全部计算参数作为随机变量来计算位置的方法。在该方法中，状态变化量矢量“δX”的四个计算参数“(δx，δy，δz，δb)”全部作为随机变量。这种情况下，式(8)的多维状态变化正态分布模型“p_δX ^(m)(δX)”成为四维状态变化正态分布模型“p_δX ⁽⁴⁾(δX)”。

最大似然估计法是统计学的方法，即，基于被称为似然值的指标值，利用给定的观测数据推算作为推算对象的参数的最大似然值。似然值是就某观测数据，用概率表示所使用的概率模型的适合程度的度量。

在本实施方式中，用式(6)的伪距误差矢量“E”的K个卫星的伪距误差“ε_k”(k＝1，2，……，K)的同时发生概率和状态变化量矢量“δX”的四个计算参数“δX_i＝(δx，δy，δz，δb)”(i＝1、2、3、4)的同时发生概率的乘积表示似然值L。即，用下式(9)表示似然值L。

[式9]

$L=\left[\underset{k=1}{\overset{K}{\mathrm{\Π}}}{p}_{\mathrm{\ϵ},k}\left({\mathrm{\ϵ}}_k\right)\right]\left[\underset{i=1}{\overset{4}{\mathrm{\Π}}}{p}_{\mathrm{\δX},i}\left({\mathrm{\δX}}_i\right)\right]$

$=\left[\underset{k=1}{\overset{K}{\mathrm{\Π}}}{p}_{\mathrm{\ϵ},k}\left({\mathrm{\ϵ}}_k\right)\right]\left[p_{\mathrm{\δX}}^{\left(4\right)}\left(\mathrm{\δX}\right)\right]...\left(9\right)$

在式(9)中，“p_ε，k(ε)”表示第k个捕获卫星的伪距误差正态分布模型。另外，在图2中定义的伪距误差正态分布模型“p_ε(ε)”可以作为对全部捕获卫星使参数值通用的通用模型，也可以作为对每个捕获卫星使参数值个性化的不同的模型。

另外，在式(9)中，“p_δX，i(δX_i)”表示第i个计算参数的状态变化正态分布模型。在从第一行向第二行变化时，改写成给出四个计算参数的同时发生概率的四维状态变化正态分布模式“p_δX ⁽⁴⁾(δX)”。

求出使式(9)的似然值L最大化的状态变化量矢量“δX”作为最大似然估计值。由于原封不动地使用似然值L进行运算很难，因此在本实施方式中，利用取其对数的对数似然值logL进行运算。如下式(10)所示地运算对数似然值logL。

[式10]

$\log L=\log \left[\underset{k=1}{\overset{K}{\mathrm{\Π}}}{p}_{\mathrm{\ϵ},k}\left({\mathrm{\ϵ}}_k\right)\right]\left[p_{\mathrm{\δX}}^{\left(4\right)}\left(\mathrm{\δX}\right)\right]$

$=\mathrm{const}-\frac{1}{2}[{(E-{\mathrm{\μ}}_E)}^T{W}_E(E-{\mathrm{\μ}}_E)+{(\mathrm{\δX}-{\mathrm{\μ}}_{\mathrm{\δX}})}^T{S}_{\mathrm{\δX}}^{-1}(\mathrm{\δX}-{\mathrm{\μ}}_{\mathrm{\δX}})]$

$=\mathrm{const}-\frac{1}{2}[{(E-{\mathrm{\μ}}_E)}^T{W}_E(E-{\mathrm{\μ}}_E)+{(\mathrm{\δX}-{\mathrm{\μ}}_{\mathrm{\δX}})}^T{W}_{\mathrm{\δX}}(\mathrm{\δX}-{\mathrm{\μ}}_{\mathrm{\δX}})]...\left(10\right)$

在式(10)中，“const”是常数项。另外，“μ_E”是以K个卫星的伪距误差“ε”的期望值“μ_ε，k＝(μ_ε，1，μ_ε，2，……，μ_ε，K)”为分量的伪距误差期望值矢量，通过下式(11)得到。

[式11]

${\mathrm{\μ}}_E=\left[\begin{array}{c}{\mathrm{\μ}}_{\mathrm{\ϵ},1}\\ {\mathrm{\μ}}_{\mathrm{\ϵ},2}\\ .\\ .\\ .\\ {\mathrm{\μ}}_{\mathrm{\ϵ},K}\end{array}\right]...\left(11\right)$

“W_E”是伪距误差逆协方差矩阵，作为由K个卫星的伪距误差“ε_k”的协方差组成的伪距误差协方差矩阵的逆矩阵进行运算。具体通过下式(12)得到。

[式12]

“μ_δX”是在式(8)中出现的状态变化量期望值矢量，通过以下的式(13)得到。

[式13]

${\mathrm{\μ}}_{\mathrm{\δX}}=\left[\begin{array}{c}{\mathrm{\μ}}_{\mathrm{\δx}}\\ {\mathrm{\μ}}_{\mathrm{\δy}}\\ {\mathrm{\μ}}_{\mathrm{\δz}}\\ {\mathrm{\μ}}_{\mathrm{\δb}}\end{array}\right]...\left(13\right)$

其中，“μ_δx，μ_δy，μ_δz，μ_δb”分别是“δx，δy，δz，δb”的期望值。

另外，“S_δX”是在式(8)中出现的状态变化量协方差矩阵，通过下式(14)得到。

[式14]

$S_{\mathrm{\δX}}=\left[\begin{array}{cccc}{\mathrm{\σ}}_{\mathrm{xx}}^2& {\mathrm{\σ}}_{\mathrm{xy}}^2& {\mathrm{\σ}}_{\mathrm{xz}}^2& {\mathrm{\σ}}_{\mathrm{xb}}^2\\ {\mathrm{\σ}}_{\mathrm{yx}}^2& {\mathrm{\σ}}_{\mathrm{yy}}^2& {\mathrm{\σ}}_{\mathrm{yz}}^2& {\mathrm{\σ}}_{\mathrm{yb}}^2\\ {\mathrm{\σ}}_{\mathrm{zx}}^2& {\mathrm{\σ}}_{\mathrm{zy}}^2& {\mathrm{\σ}}_{\mathrm{zz}}^2& {\mathrm{\σ}}_{\mathrm{zb}}^2\\ {\mathrm{\σ}}_{\mathrm{bx}}^2& {\mathrm{\σ}}_{\mathrm{by}}^2& {\mathrm{\σ}}_{\mathrm{bz}}^2& {\mathrm{\σ}}_{\mathrm{bb}}^2\end{array}\right]...\left(14\right)$

其中，4×4矩阵的各分量是“δx，δy，δz，δb”的四个分量之间的协方差。

另外，“W_δX”是状态变化量逆协方差矩阵，作为状态变化量协方差矩阵“S_δX”的逆矩阵进行运算。具体通过下式(15)得到。

[式15]

$W_{\mathrm{\δX}}=S_{\mathrm{\δX}}^{-1}...\left(15\right)$

为了求出使式(10)最大化的“δX”，解方程“d(logL)/d(δX)＝0”即可。省略中间运算过程，最终通过下式求出最大似然估计值“δX”。

[式16]

δX＝(G^TW_EG+W_δX)^-1(G^TW_E(δρ-μ_E)+W_δXμ_δX)…(16)

其中，“δρ”是式(4)的伪距变化量矢量，“G”是式(5)的几何矩阵。

分别就每一个运算时刻，根据式(16)求出状态变化量矢量“δX”的最大似然估计值。然后利用求出的状态变化量矢量“δX”的最大似然估计值计算状态矢量“X”。

另外，在本实施方式的位置计算运算中，计算上述最大似然估计值的同时，基于过去的位置计算运算的结果对状态变化正态分布模型进行更新处理。在此，将运算时刻“t”中的状态变化量矢量“δX”标记为“δX(t)”。并且，将运算时刻“t”中的状态变化量期望值矢量“μ_δX”和状态变化量协方差矩阵“S_δX”分别标记为“μ_δX(t)”和“S_δX(t)”。

这种情况下，根据下式(17)和(18)分别对状态变化量期望值矢量“μ_δX”和状态变化量协方差矩阵“S_δX”进行更新。

[式17]

μ_δX(t+1)＝δX(t)…(17)

[式18]

$S_{\mathrm{\δX}}(t+1)={(G^T{W}_{E}G+W_{\mathrm{\δX}})}^{-1}[G^T{W}_{E}G+W_{\mathrm{\δX}}{S}_{\mathrm{\δX}}\left(t\right)W_{\mathrm{\δX}}^T]{(G^T{W}_{E}G+W_{\mathrm{\δX}})}^{-1}...\left(18\right)$

式(17)表示用上一次运算时刻计算的状态变化量矢量“δX(t)”更新这次的运算时刻的状态变化量期望值矢量“μ_δX(t+1)”。这相当于基于过去的位置计算运算的结果更新运算结果(状态变化量矢量“δX”)的变化标准(状态变化量期望值矢量“μ_δX(t)”)。另外，式(18)表示用上一次运算时刻运算的状态变化量协方差矩阵“S_δX(t)”等各量更新这次的运算时刻的状态变化量协方差矩阵“S_δX(t+1)”。

状态变化量期望值矢量“μ_δX”是与状态变化正态分布模型的期望值相关的矢量。另外，状态变化量协方差矩阵“S_δX”是与状态变化正态分布模型的标准差(协方差)相关的矩阵。因此，可以将式(17)和(18)称为基于过去的位置计算运算的结果对状态变化正态分布模型进行更新的式子。

(2)第二位置计算方法

第二位置计算方法是伪距误差“ε”使用正态分布模型、将状态变化量矢量“δX”的一部分计算参数作为随机变量进行位置计算的方法。在该方法中，使构成状态变化量矢量“δX”的四个计算参数“δx、δy、δz、δb”中的一部分计算参数为随机变量，其余的计算参数为非随机变量。

根据所使用的系统，可以适当地选择作为随机变量和非随机变量的计算参数。这里，设构成使用者的位置变化量矢量“δP”的“(δx，δy，δz)”的三个计算参数为随机变量，设钟差变化量“δb”为非随机变量，以此为例进行说明。这种情况下，式(8)的多维状态变化正态分布模型“p_δX ^(m)(δX)”成为三维状态变化正态分布模型“p_δX ⁽³⁾(δX)”。

这种情况下，由于未假定钟差变化量“δb”是正态分布模型，因此，将相当于式(10)的对数似然值“logL”更改为下式(19)。

[式19]

$\log L=\log \left[\underset{k=1}{\overset{K}{\mathrm{\Π}}}{p}_{\mathrm{\ϵ},k}\left({\mathrm{\ϵ}}_k\right)\right]\left[p_{\mathrm{\δX}}^{\left(3\right)}\left(\mathrm{\δX}\right)\right]$

$=\mathrm{const}-\frac{1}{2}[{(E-{\mathrm{\μ}}_E)}^T{W}_E(E-{\mathrm{\μ}}_E)+{(\mathrm{\δX}-{\mathrm{\μ}}_{\mathrm{\δX}})}^T{S}_{\mathrm{\δX}}^{-1}(\mathrm{\δX}-{\mathrm{\μ}}_{\mathrm{\δX}})]$

$=\mathrm{const}-\frac{1}{2}[{(E-{\mathrm{\μ}}_E)}^T{W}_E(E-{\mathrm{\μ}}_E)+{(\mathrm{\δX}-{\mathrm{\μ}}_{\mathrm{\δX}})}^T{W}_{\mathrm{\δX}}(\mathrm{\δX}-{\mathrm{\μ}}_{\mathrm{\δX}})]...\left(19\right)$

另外，将在式(13)中定义的状态变化量期望值矢量“μ_δX”更改为下式(20)。

[式20]

${\mathrm{\μ}}_{\mathrm{\δX}}=\left[\begin{array}{c}{\mathrm{\μ}}_{\mathrm{\δx}}\\ {\mathrm{\μ}}_{\mathrm{\δy}}\\ {\mathrm{\μ}}_{\mathrm{\δz}}\\ 0\end{array}\right]...\left(20\right)$

其中，由于设钟差变化量“δb”为非随机变量，因此使与其对应的第四个分量为零。

另外，将在式(14)中定义的状态变化量协方差矩阵“S_δX”更改为下式(21)。

[式21]

$S_{\mathrm{\δX}}=\left[\begin{array}{cccc}{\mathrm{\σ}}_{\mathrm{xx}}^2& {\mathrm{\σ}}_{\mathrm{xy}}^2& {\mathrm{\σ}}_{\mathrm{xz}}^2& 0\\ {\mathrm{\σ}}_{\mathrm{yx}}^2& {\mathrm{\σ}}_{\mathrm{yy}}^2& {\mathrm{\σ}}_{\mathrm{yz}}^2& 0\\ {\mathrm{\σ}}_{\mathrm{zx}}^2& {\mathrm{\σ}}_{\mathrm{zy}}^2& {\mathrm{\σ}}_{\mathrm{zz}}^2& 0\\ 0& 0& 0& 0\end{array}\right]...\left(21\right)$

其中，由于设钟差变化量“δb”为非随机变量，因此使4×4矩阵的第四行和第四列的分量都为零。

这种情况下，也可以如式(16)所示地求出最大似然估计值“δX”。另外，分别根据式(17)和式(18)更新在式(20)中定义的状态变化量期望值矢量“μ_δX”和在式(21)中定义的状态变化量协方差矩阵“S_δX”。这种情况下，与钟差变化量“δb”对应的分量为零，不进行更新。

另外，对作为非随机变量的钟差变化量“δb”使用接收信号进行预定的收敛运算求解。具体来说，对式(3)的观测方程式使用伪距误差正态分布模型，例如利用牛顿法进行收敛运算，将求出钟差变化量“δb”的近似解。通过收敛运算计算的钟差变化量“δb”在下一次的运算时刻的位置计算运算中用于近似地求出伪距。

1-2.如果使用伪距误差正态混合分布模型(伪距误差GMM)

(3)第三位置计算方法

第三位置计算方法是伪距误差“ε”使用正态混合分布模型、将状态变化量矢量“δX”的全部计算参数作为随机变量进行位置计算的方法。概念与第一位置计算方法相同，但不同点在于伪距误差的概率分布模型不是正态分布模型，而使用正态混合分布模型。

这种情况下，利用伪距误差GMM“f(ε)”和四维状态变化正态分布模型“p_δX ⁽⁴⁾(δX)”将对数似然值“logL”定型为下式(22)。

[式22]

$\log L=\log \left[\underset{k=1}{\overset{K}{\mathrm{\Π}}}{f}_k\left({\mathrm{\ϵ}}_k\right)\right]\left[\underset{i=1}{\overset{4}{\mathrm{\Π}}}{p}_{\mathrm{\δX},i}\left({\mathrm{\δX}}_i\right)\right]$

$=\log \left[\underset{k=1}{\overset{K}{\mathrm{\Π}}}{f}_k\left({\mathrm{\ϵ}}_k\right)\right]\left[p_{\mathrm{\δX}}^{\left(4\right)}\left(\mathrm{\δX}\right)\right]...\left(22\right)$

在式(22)中，“f_k(ε_k)”表示第k个捕获卫星的伪距误差GMM，一般化为下式(23)。

[式23]

$f_k\left({\mathrm{\ϵ}}_k\right)=\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{\mathrm{\α}}_l{N}_l[{\mathrm{\μ}}_k^{\left(l\right)},{\left({\mathrm{\σ}}_k^{\left(l\right)}\right)}^2]$

$=\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{\mathrm{\α}}_l\left[\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_k^{\left(l\right)}}\exp (-\frac{{({\mathrm{\ϵ}}_k-{\mathrm{\μ}}_k^{\left(l\right)})}^2}{2{\left({\mathrm{\σ}}_k^{\left(l\right)}\right)}^2})\right]...\left(23\right)$

在式(23)中，“l”是混合(合成)的正态分布模型的序号，记载为在混合L个(l＝1，2，……，L)正态分布模型的情况下的式子。“α_l”是第l个正态分布模型“N_l”的比重(合成比例)。“μ_k ^(l)”是与第k个捕获卫星相关的第l个正态分布模型的期望值，“σ_k ^(l)，是与第k个捕获卫星相关的第l个正态分布模型的标准差。

并且，上述伪距误差GMM“f(ε)”可以作为对全部捕获卫星使参数通用的通用模型，也可以作为对每个捕获卫星使参数值个性化的不同的模型。

为了求出使式(22)的对数似然值“logL”最大化的“δX”，解方程“d(logL)/d(δX)＝0”即可。但是，在正态混合分布(GMM)中，由于存在对数的求和等复杂的部分，因此很难解析地求解。因此，在本实施方式中，使用EM(Expectation Maximum，期望最大化)算法(其为对于优化问题的一种算法)求出最大似然估计值。

EM算法是基于最大似然估计法推断概率分布模型的参数的一种方法。在除了可观测数据(样本)以外，还存在不可观测的数据(隐藏参数)的情况下，在该隐藏参数取决于概率分布模型时使用的方法。EM算法是一种迭代法，重复两个步骤，即，计算与似然函数的条件概率相关的期望值的E(期望)步骤以及求出使期望值最大化的解的M(最大化)步骤，以求出最大似然估计值。

省略中间运算过程，如果使用EM算法，通过下式(24)求出状态变化量矢量“δX”的最大似然估计值。

[式24]

$\mathrm{\δX}={(G^T\left(\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{M}_l\left({\mathrm{\δX}}_0\right)\right)G+W_{\mathrm{\δX}})}^{-1}(G^T\left(\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{M}_l\left({\mathrm{\δX}}_0\right)(\mathrm{\δ\ρ}-{\mathrm{\μ}}_l)\right)+W_{\mathrm{\δX}}{\mathrm{\μ}}_{\mathrm{\δX}})...\left(24\right)$

其中，“δX₀”是各运算时刻的状态变化量矢量“δX”的初始值。

通过下式(25)得到式(24)中的“M_l(δX₀)”。

[式25]

另外，通过下式(26)得到式(25)中的“h_k ^(l)(δX₀)”。

[式26]

另外，通过下式(27)得到式(24)中的“μ_l”。

[式27]

${\mathrm{\μ}}_l=\left[\begin{array}{c}{\mathrm{\μ}}_1^{\left(l\right)}\\ {\mathrm{\μ}}_2^{\left(l\right)}\\ .\\ .\\ .\\ {\mathrm{\μ}}_K^{\left(l\right)}\end{array}\right]...\left(27\right)$

在每个运算时刻，根据式(24)至式(27)求出最大似然估计值“δX”。另外，与使用伪距误差正态分布模型时一样，在每个运算时刻更新在式(13)中定义的状态变化量期望值矢量“μ_δX”和在式(14)中定义的状态变化量协方差矩阵“S_δX”。

具体地，分别根据下式(28)和(29)更新状态变化量期望值矢量“μ_δX(t)”和状态变化量协方差矩阵“S_δX(t)”。

[式28]

μ_δX(t+1)＝δX(t)…(28)

[式29]

S_δX(t+1)＝H(δX(t))F(δX)^TH(δX(t))^T…(29)

通过下式(30)和(31)分别得到式(29)中的“H(δX)”和“F(δX)”。

[式30]

$H\left(\mathrm{\δX}\right)={(G^T\left(\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{M}_l\left(\mathrm{\δX}\right)\right)G+W_{\mathrm{\δX}})}^{-1}...\left(30\right)$

[式31]

$F\left(\mathrm{\δX}\right)=G^T\left(\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{M}_l\left(\mathrm{\δX}\right)R_l{M}_l{\left(\mathrm{\δX}\right)}^T\right)G+W_{\mathrm{\δX}}{S}_{\mathrm{\δX}}\left(t\right)W_{\mathrm{\δX}}^T...\left(31\right)$

并且，通过下式(32)得到式(31)中的“R_l”。

[式32]

(4)第四位置计算方法

第四位置计算方法是伪距误差“ε”使用正态混合分布模型、将状态变化量矢量“δX”的一部分计算参数作为随机变量进行位置计算的方法。在该方法中，与第二位置计算方法一样，设构成状态变化量矢量“δX”的四个计算参数“δx、δy、δz、δb”中的一部分计算参数为随机变量，其余的计算参数为非随机变量。与第二位置计算方法的不同点在于伪距误差的概率分布模型不是正态分布模型，而使用正态混合分布模型。

这种情况下，也可以使用与第三位置计算方法相同的数学式求出状态变化量矢量“δX”的最大似然估计值。与第三位置计算方法的不同点在于，对作为随机变量处理的计算参数只使用正态分布模型进行运算。即，与第二位置计算方法相同，设状态变化量矢量“δX”的分量中的与作为非随机变量处理的计算参数对应的分量为零进行运算。

例如，设钟差变化量“δb”为非随机变量的情况下，使与钟差变化量“δb”对应的分量为零进行运算。并且，对钟差变化量“δb”使用接收信号进行预定的收敛运算来求解。具体是对式(3)的观测方程式使用伪距误差正态混合分布模型(伪距误差GMM)，例如利用EM算法进行收敛运算，将求出钟差变化量“δb”的近似解。此时，通过收敛运算计算的钟差变化量“δb”在下一次的运算时刻的位置计算运算中用于近似地求出伪距。

2.实验结果

就利用本实施方式的位置计算方法进行位置计算运算的实验结果进行说明。分别使用现有方法与本实施方式的方法进行位置计算运算，进行实际计算位置的实验。使用者携带位置计算装置，一面沿着事先确定的绕行路径步行一面计算位置。

图9是表示位置计算结果的一个例子的实验结果。图表的横向表示东西方向，纵向表示南北方向。单位都是米“m”。这里所示的结果是状态变化的概率分布模型使用正态分布、伪距误差的概率分布模型使用正态混合分布(GMM)的结果。白色圆圈的标记表示通过现有的位置计算方法得到的轨迹，黑色圆圈的标记表示通过本实施方式的位置计算方法得到的轨迹。

从该结果可见，在使用现有的位置计算方法时，有些地方发生跳位，未能获得非常准确的位置。但在使用本实施方式的位置计算方法时，得到了沿着绕行路径的平滑轨迹。通过不是将位置计算的运算结果独立、而是使运算结果具有时间上的相关性来进行位置计算，计算位置的连续性从而有飞跃性地提高。

图10是在上述步行实验中测量位置运算精度的图表。图表的横向表示东西方向的位置误差，纵向表示南北方向的位置误差。单位都是米(m)。图表的中心相当于位置误差“0(米)”，标记表示越靠近图表的中心位置计算精度越高。

从图表可见，在使用现有的位置计算方法时，标记整体上差别很大，位置计算精度降低。相反，在使用本实施方式的位置计算方法时，标记整体上集中在中心部，与现有的方法相比，位置计算精度明显提高。

图11是表示位置计算精度的另外的实验结果的一个例子。这里所示的结果是状态变化的概率分布模型和伪距误差的概率分布模型都使用了正态分布模型的结果。图示方式与图10相同。从该结果可见，在伪距误差的概率分布模型使用正态分布模型的情况下，表示位置误差的标记集中在图表的中心部，位置计算精度明显提高。

通过上述实验结果，验证了通过使用本实施方式的位置计算方法，位置计算的准确性得以提高。

3.实施例

以下就利用上述的位置计算方法进行位置计算的位置计算装置的实施例进行说明。这里以作为具有位置计算装置的电子设备的一个例子的移动电话作为实施例进行说明。当然，可以使用本发明的实施例不限于以下说明的实施例。

3-1.移动电话的结构

图12是示出移动电话1的功能结构的一个例子的框图。移动电话1被构成为包括GPS天线5、GPS接收部10、主机处理部30、操作部40、显示部50、移动电话用天线60、移动电话用无线通信电路部70、存储部80以及时钟部90。

GPS天线5是接收包含从GPS卫星发送的GPS卫星信号的RF(Radio Frequency)信号的天线，向GPS接收部10输出接收信号。

GPS接收部10是基于从GPS天线5输出的信号计算移动电话1的位置的电路或装置，是相当于GPS接收器的功能块。在本实施例中，该GPS接收部10相当于位置计算装置。

GPS接收部10被构成为包括RF接收电路部11和基带处理电路部20。GPS接收部10相当于接收GPS卫星信号的接收部。并且，RF接收电路部11和基带处理电路部20可以作为分开的LSI(大规模集成电路)制造，也可以作为一个芯片来制造。

RF接收电路部11是RF信号的接收电路。电路结构例如可以构成为利用A/D转换器将从GPS天线5输出的RF信号转换成数字信号，并处理数字信号的接收电路的结构。另外，也可以构成为在保持为模拟信号的情况下对从GPS天线5输出的RF信号进行信号处理，最后再利用A/D转换器将数字信号输出至基带处理电路部20的结构。

在后一种情况下，可以如下所述地构成RF接收电路部11。即，通过分频或倍频预定的振荡信号，生成用于RF信号的乘法的振荡信号。然后，将生成的振荡信号与从GPS天线输出的RF信号相乘，从而将RF信号下转换到中间频率的信号(以下称为“IF(中频)信号”)，然后放大IF信号后通过A/D转换器转换成数字信号，向基带处理电路部20输出。

基带处理电路部20对从RF接收电路部11输出的接收信号进行载波去除和相关运算等以捕获GPS卫星信号。然后，利用从捕获到的GPS卫星信号中提取的时刻信息或卫星轨道信息等来计算移动电话1的位置和时钟误差。

主机处理部30是根据存储于存储部80的系统程序等各种程序统一控制移动电话1的各部分的处理器，被构成为具有CPU(中央处理单元)等处理器。主机处理部30根据从基带处理电路部20获取的位置坐标，在显示部50上显示指示当前位置的地图，或将该位置坐标用于各种应用程序处理。

操作部40是例如由触摸屏或按钮开关等构成的输入装置，将所按下的键或按钮的信号向主机处理部30输出。通过操作该操作部40，输入通话请求、邮件发送请求、位置计算请求等各种指示。

显示部50由LCD(液晶显示器)等构成，是基于从主机处理部30输入的显示信号进行各种显示的显示装置。在显示部50上显示位置显示页面和时间信息等。

移动电话用天线60是在移动电话1与其通信服务提供商设置的无线基站之间进行移动电话用无线信号的收发的天线。

移动电话用无线通信电路部70是由RF转换电路、基带处理电路等构成的移动电话的通信电路部，通过进行移动电话用无线信号的调制和解调等，实现通话或邮件的收发等。

存储部80被构成为具有ROM(只读存储器)或闪存ROM、RAM(随机存储器)等存储装置，存储主机处理部30用以控制移动电话1的系统程序、用以执行各种应用程序处理的各种程序或数据等。

时钟部90是移动电话1的内部时钟，被构成为具有晶体振荡器等，晶体振荡器由晶体振子和振荡电路构成。时钟部90的计时时刻被随时向基带处理电路部20和主机处理部30输出。时钟部90使用基带处理电路部20计算的时钟误差进行校正。

3-2.基带处理电路部的结构

图13是基带处理电路部20的电路结构以及存储部23的数据结构的说明图。基带处理电路部20具有处理部21和存储部23作为主要的功能结构。

处理部21是统一控制基带处理电路部20的各功能部的控制装置和运算装置，被构成为具有CPU或DSP(数字信号处理器)等处理器。处理部21具有卫星捕获部211和位置计算部213作为功能部。

卫星捕获部211是执行GPS卫星的捕获的功能部。具体地，对从RF接收信号部11输出的被数字化的接收信号进行载波去除和相关运算的数字信号处理，以捕获GPS卫星。并且对相关运算结果进行峰值判断，获取接收载波信号的多普勒频率或接收频率、接收C/A码的码相位等信息作为运算信息。

另外，卫星捕获部211基于相关运算结果对导航数据进行解码。检测接收载波信号的相位(载波相位)和接收C/A码的相位(码相位)，一旦成为获得了相关的状态，就可以根据相关值的时间变化对构成导航数据的各比特值进行解码。例如通过被称为锁相环的PLL(锁相环)处理进行该相位的同步。

位置计算部213利用与由卫星捕获部211捕获的各GPS卫星有关的测量信息或导航数据、时刻信息、卫星信息等各个量(諸量)，根据上述原理进行位置计算运算，计算移动电话1的位置(位置坐标)和时钟误差(钟差)。然后向主机处理部30输出计算的位置，同时利用计算的时钟误差来校正时钟部90。

位置计算部213相当于位置计算部，其利用GPS接收部10的接收信号进行位置计算运算，该位置计算运算至少基于将运算结果的变化作为随机变量的给定的概率分布模型而确定。

存储部23存储基带处理电路部20的系统程序、用于实现卫星捕获功能、位置计算功能等各种功能的各种程序、数据等。另外，具有临时存储各种处理的中间处理数据、处理结果等的工作区。

存储部23存储由处理部21读出并作为基带处理(参照图14)来执行的第一基带处理程序231作为程序。另外，基带处理程序231包括作为各种位置计算处理(参照图15至图18)来执行的位置计算程序2311作为子例程。存储与在原理中说明的第一至第四位置计算方法中的任一个对应的程序作为位置计算程序2311。

另外，作为主要数据，存储部23存储卫星轨道数据232、捕获到的各GPS卫星相关的每个卫星的数据233、状态变化正态分布模型234以及运算结果数据235。

卫星轨道数据232是存储了全部GPS卫星的大体的卫星轨道信息的年历、存储了各GPS卫星的详细的卫星轨道信息的星历等数据。该卫星轨道数据232除了通过对从GPS卫星接收的GPS卫星信号进行解码来获取以外，例如从移动电话1的无线基站或辅助服务器获取以作为辅助数据。

每个卫星的数据233包括测量信息233A、卫星信息233B以及定位计算用的各个量233C。测量信息233A是对该捕获卫星进行测量所得的GPS卫星信号的码相位、接收频率、多普勒频率等信息。卫星信息233B是该捕获卫星的位置、速度、移动方向等信息。定位计算用的各个量233C是根据上述原理进行定位计算时使用的各个量。

状态变化正态分布模型234是在图1中图示并说明的状态变化正态分布模型，模型参数值例如包括状态变化量期望值234A和状态变化量标准差234B。状态变化正态分布模型在每个运算时刻进行更新。

运算结果数据235是作为位置计算运算的结果得到的数据，包括移动电话1的位置矢量“P”和钟差“b”。

3-3.处理流程

图14是示出通过由处理部21读出存储在存储部23的基带处理程序231而在基带处理电路部20中执行的基带处理流程的流程图。

首先，处理部21根据作为子例程存储在存储部23中的位置计算程序2311开始位置计算处理(步骤A1)。位置计算处理是图15至图18的第一至第四位置计算处理中的任一个处理。

然后，处理部21判断是否是运算结果的输出定时(步骤A3)。可以将运算结果的输出定时设定为任意的定时。例如可以是与位置计算运算的时间间隔相同的时间间隔的定时，也可以是比位置计算运算的时间间隔长的时间间隔的定时。另外，也可以将使用者指示输出运算结果的定时作为输出定时。

如果判断是输出定时(步骤A3中的是)，处理部21就向主机处理部30输出最新的运算结果(步骤A5)。在步骤A3中如果判断出不是输出定时(步骤A3中的否)，或在步骤A5之后，处理部21判断是否结束基带处理(步骤A7)。

如果判断出不结束基带处理(步骤A7中的否)，则处理部21返回步骤A3。而如果判断出结束(步骤A7中的是)，则处理部21结束基带处理。

图15是示出与第一位置计算方法对应的第一位置计算处理的流程的流程图。首先，处理部21设定状态矢量“X”的初始值(步骤B1)。具体地，将通过服务器辅助从移动电话1的无线基站获取的位置(基站位置)或预先确定的位置(例如固定值)设定为初始位置。另外，对钟差“b”也设定预定的初始值(例如固定值)。

然后，处理部21进行捕获目标卫星选择处理(步骤B3)。具体地，利用存储在存储部23中的年历和星历等卫星轨道数据232，判断在时钟部90进行了计时的当前时刻位于初始位置的天空中的GPS卫星，选择作为捕获目标卫星。

然后，处理部21进行状态变化量估计处理(步骤B5至B17)。具体地，处理部21对各捕获目标卫星执行循环A(步骤B5至B13)。在循环A中，处理部21从该捕获目标卫星接收GPS卫星信号并测量伪距，将该数值作为伪距测量值“ρ_c”(步骤B7)。可以利用作为相位搜索的结果而得到的码相位来计算伪距测量值“ρ_c”。

另外，处理部21计算伪距近似值“ρ_a”(步骤B9)。具体地，利用从卫星轨道数据232求出的该捕获目标卫星的位置矢量“P_k”以及移动电话1的最新的位置矢量“P”和最新的钟差“b”，计算伪距近似值“ρ_a＝||P_k-P||+b”。其中，在该式中，将钟差“b”的单位换算成距离。

然后，处理部21计算伪距变化量“δρ”(步骤B11)。具体地，计算在步骤B7中测量的伪距测量值“ρ_c”和在步骤B9中计算的伪距近似值“ρ_a”的差，作为伪距变化量“δρ”。然后，处理部21转向对下一个捕获目标卫星的处理。对全部捕获目标卫星进行了步骤B7至B11的处理之后，处理部21结束循环A(步骤B13)。

然后，处理部21计算几何矩阵G(步骤B15)。具体地，利用从卫星轨道数据232求出的各捕获目标卫星的卫星位置矢量和移动电话1的最新的位置矢量，计算从移动电话1向着各捕获目标卫星的视线矢量“l”。然后求出式(5)的几何矩阵G。

然后，处理部21使用伪距误差正态分布模型和状态变化正态分布模型求出状态变化量“δX”的最大似然估计值(步骤B17)。具体地，根据式(16)计算状态变化量“δX”。由于假定伪距误差“ε”是正态分布，因此这里得到的最大似然估计值相当于最小二乘解。之后处理部21结束状态变化量估计处理。

然后，处理部21利用在步骤B17中作为最大似然估计值估计的状态变化量矢量“δX”来校正状态矢量“X”(步骤B19)。具体地，将作为最大似然估计值求出的状态变化量矢量“δX”与最新的状态矢量“X”相加，从而重新计算这次的运算时刻的状态矢量“X”。

然后，处理部21判断是否结束位置计算(步骤B23)，如果判断出未结束(步骤B23中否)，则对状态变化正态分布模型进行更新处理(步骤B25)。具体地，分别通过式(17)和式(18)更新状态变化量期望值矢量“μ_δX”和状态变化量协方差矩阵“S_δX”。

在进行了状态变化正态分布模型的更新处理之后，处理部21返回步骤B3，转移到下一个运算时刻的位置计算运算。另外，如果在步骤B23中判断出位置计算结束(步骤B23中的是)，则处理部21结束第一位置计算处理。

从步骤B3到B25的一系列处理相当于在一个运算时刻的位置计算运算。第一位置计算处理中的特征在于，在每个运算时刻的位置计算运算中省略了收敛运算。现有的位置计算方法是在利用最小二乘法求解时，例如利用牛顿法进行收敛运算求出近似解。但是，本实施方式中将运算结果的变化作为随机变量，并且在每个运算时刻更新概率分布模型，因此无需进行收敛运算。

图16是示出与第二位置计算方法对应的第二位置计算处理的流程的流程图。是将状态变化量矢量“δX”中的位置变化量矢量“δP”的分量“δx、δy、δz”设为随机变量、将钟差变化量“δb”设为非随机变量进行位置计算的流程图。与图15的第一位置计算处理相同的步骤使用相同的符号并省略说明。

处理部21在步骤B3中进行了捕获目标卫星选择处理之后，进行状态变化量估计处理(步骤C5)。具体地，处理部21进行与第一位置计算处理的状态变化量估计处理(步骤B5至B17)相同的处理。其中，由于将这里的钟差“δb”设为非随机变量，因此使状态变化量矢量“δX”的与钟差变化量“δb”对应的分量为零来计算。即，状态变化量估计处理成为估计位置变化量矢量“δP”的位置变化量估计处理。

另外，在状态变化量估计处理中，利用在上一次的位置计算运算中通过收敛运算计算的钟差“b”的值计算伪距近似值“ρ_a”(步骤B9)。具体地，利用从卫星轨道数据232求出的该捕获目标卫星的位置矢量“P_k”、在上一次的位置计算运算中计算的位置矢量“P”以及在上一次的位置计算运算中通过收敛运算计算的钟差“b”，来计算伪距近似值“ρ_a＝||P_k-P||+b”。这相当于利用过去通过收敛运算计算的预定的参数的值来进行位置计算运算。

然后，处理部21进行钟差变化量估计处理(步骤C7至C19)。钟差变化量估计处理是利用GPS卫星信号的接收信号进行预定的收敛运算、求出作为非随机变量的钟差变化量“δb”的近似解的处理。

具体地，处理部21对各捕获目标卫星执行循环B(步骤C7至C13)。在循环B中，处理部21计算伪距近似值“ρ_a”(步骤C9)。然后，处理部21计算伪距变化量“δρ”(步骤C11)。具体地，计算在状态变化量估计处理步骤B7中测量的伪距测量值“ρ_c”与在步骤C9中计算的伪距近似值“ρ_a”的差，作为伪距变化量“δρ”。

然后，处理部21转向对下一个捕获目标对象的处理。对所有的捕获目标卫星进行了步骤C9至C11之后，处理部21结束循环B(步骤C13)。

然后，处理部21使用伪距误差正态分布模型求出钟差变化量“δb”的最小二乘解(步骤C17)。然后处理部21就求出的钟差变化量“δb”判断预定的收敛条件是否成立(步骤C18)。具体地，例如，如果在这次迭代步骤中求出的钟差变化量“δb”与在上一次的迭代步骤中求出的钟差变化量“δb”的差的绝对值小于预定阈值，则判断收敛条件成立。

如果判断收敛条件不成立(步骤C18中的否)，则处理部21利用在步骤C17中求出的钟差变化量“δb”，校正并更新钟差“b”(步骤C19)。然后处理部21返回步骤C7。另外，如果判断收敛条件成立(步骤C18中的是)，则处理部21结束钟差变化量估计处理。

然后，处理部21将在状态变化量估计处理中估计的位置变化量矢量“δP”和在钟差变化量估计处理中估计的时钟变化量“δb”作为状态变化量矢量“δX”，用其校正状态矢量“X”(步骤C21)。即，将状态变化量矢量“δX”与最新的状态矢量“X”相加，重新计算这次的运算时刻的状态矢量“X”。

然后，处理部21判断是否结束位置计算(步骤B23)，如果判断出还未结束(步骤B23中的否)，则对状态变化正态分布模型进行更新处理(步骤B25)。然后，处理部21返回步骤B3，转向在下一个运算时刻的位置计算运算。另外，如果在步骤B23中判断出位置计算结束(步骤B23中的是)，处理部21就结束第二位置计算处理。

在此情况下，从步骤B3到B25的一系列处理也相当于在一个运算时刻的位置计算运算。另外，步骤C7至C19的反复处理相当于使用接收信号的预定收敛运算。第二位置计算运算处理中的特征在于，对作为随机变量处理的位置变化量矢量“δP”省略了收敛运算，但对作为非随机变量处理的钟差变化量“δb”未省略收敛运算。这是因为，由于钟差变化量“δb”不作为随机变量，而作为完全的未知数处理，因此如果不进行收敛运算就不能得到最佳值。

图17是示出与第三位置计算方法对应的第三位置计算处理的流程的流程图。第三位置计算处理的基本流程与图15的第一位置计算处理相同。但是，第一位置计算处理的步骤B17置换为步骤D17，在这点上有很大的不同。

在步骤D17中，处理部21利用伪距误差正态混合分布模型和状态变化正态分布模型，求出状态变化量矢量“δX”的最大似然估计值。具体地，处理部21根据式(24)至(27)计算状态变化量矢量“δX”的最大似然估计值。该步骤是使用EM算法的计算步骤。

另外，在第三位置计算处理中，在步骤B25的状态变化正态分布模型更新处理中，处理部21根据式(28)更新状态变化量期望值矢量“μ_δX”。另外，处理部21根据式(29)至式(32)更新状态变化量协方差矩阵“S_δX”。

图18是示出与第四位置计算方法对应的第四位置计算处理的流程的流程图。第四位置计算处理的基本流程与图16的第二位置计算处理相同。但是，第二位置计算处理的步骤C5置换为步骤E5、步骤C17置换为步骤E17，在这点上有很大的不同。

在第四位置计算处理中，在步骤E5的状态变化量估计处理中，将钟差“δb”作为非随机变量处理，设状态变化量矢量“δX”的与钟差变化量“δb”对应的分量为零来计算。即，此时的状态变化量估计处理成为估计位置变化量矢量“δP”的位置变化量估计处理。

另外，与第二位置计算处理相同，在状态变化量估计处理中，使用在上一次的位置计算运算中通过收敛运算计算的钟差“b”的值计算伪距近似值“ρ_a”(步骤B9)。即，使用通过过去的收敛运算计算的预定参数的值进行位置计算运算。

另外，在第四位置计算处理中，在步骤C7至C19的钟差变化量估计处理中，使用伪距误差正态混合分布模型，例如利用EM算法进行收敛运算(步骤C18中的否→步骤C19→步骤C7)，求出钟差变化量“δb”的最大似然估计值(步骤E17)。

然后，处理部21就钟差变化量“δb”判断预定的收敛条件是否成立(步骤C18)，如果判断出不成立(步骤C18中的否)，就利用在步骤E17中求出的钟差变化量“δb”校正并更新钟差“δ”(步骤C19)，之后返回步骤C7。另外，如果判断出收敛条件成立(步骤C18中的是)，处理部21就转向步骤C21。

4.作用效果

根据本实施方式，接收作为定位用卫星信号中的一种的GPS卫星信号。并且进行位置计算运算，其是利用GPS卫星信号的接收信号的运算，是基于将GPS接收器的位置和钟差的变化作为随机变量的给定的概率分布模型而确定的运算。

具体地，设定将GPS接收器(位置计算装置)的位置和钟差的变化作为随机变量的状态变化正态分布模型。另外，将伪距中包含的误差的变化作为随机变量的概率分布模型设定为伪距误差正态分布模型或伪距误差正态混合分布模型。并且，利用这些概率模型，基于最大似然估计法进行运算，求出位置和钟差变化量的最大似然估计值。通过这样，能够不使每个运算时刻的运算结果独立、而使运算结果具有时间上的相关性地进行位置计算。这种情况下，即使在测量的伪距中混合了瞬间很大的误差，也可以防止位置计算精度下降。

另外，在本实施方式中，基于过去的位置计算运算的结果更新运算结果的变化标准，利用更新后的标准计算运算结果的变化。具体地，利用在上一次的运算时刻计算的位置和钟差的变化量来更新这次的运算时刻的位置和钟差的变化量的标准。这相当于更新状态变化正态分布模型的期望值。

另外，在本实施方式中，利用在上一次的运算时刻的位置和钟差的变化量的协方差更新在这次的运算时刻的位置和钟差的变化量的协方差。这相当于更新状态变化正态分布模型的标准差。即，不是固定地设定与运算结果的变化相关的正态分布模型，而是基于上一次的运算时刻的运算结果更新该正态分布模型。通过这样，可以敏捷地对应位置和钟差的时间变化，可以更准确地得到移动体的位置。

5.变形例

5-1.多普勒定位

在上述实施方式中，以利用伪距的位置计算方法为例进行了说明，但对于利用多普勒的位置计算方法实际上也同样可以使用本发明。将利用多普勒的位置计算方法称为“多普勒定位”。

关于多普勒，下式(33)的模型式成立。

[式33]

$h_k=(V_k-V)\·l_k+d+{\mathrm{\ϵ}}_k^h...\left(33\right)$

其中，将单位全部作为速度的量纲(m/s)来记载该式。

在式(33)中，“h_k”是就第k个卫星测量的频率的漂移(以下称为“频率漂移”)。即，作为GPS卫星信号的接收信号的频率(接收频率)与预定的载波频率(1.57542[GHz])的差来进行计算。可以通过进行频率方向的相关运算获得接收频率。

“V_k”是第k个捕获卫星的速度矢量，“V”是GPS接收器的速度矢量。因此，“V_R＝V_k-V”是GPS接收器与第k个捕获卫星的相对速度矢量。另外，“l_k”是从GPS接收器向着第k个捕获卫星的视线矢量。通过向视线矢量“l_k”的方向投影相对速度矢量“V_R＝V_k-V”，从而得到视线方向上的相对速度矢量(多普勒频移)。

另外，在式(33)中，“d”是GPS接收器的工作时钟漂移(时钟漂移)。另外，“ε_k ^h”是第k个卫星的频率漂移固有的观测误差。

考虑对式(33)使用牛顿法。定义以GPS接收器的位置矢量“P”、钟差“b”、速度矢量“V”以及时钟漂移“d”为分量的状态矢量“X＝(P，b，V，d)^T”。即，除了位置信息，也将钟差、速度信息以及时钟漂移作为状态矢量“X”的分量。

在此，定义以式(33)的(右边)-(左边)表示的下式(34)的函数“g(P，V)”。

[式34]

$g_k(P,V)=(V_k-V)\·l_k+d-h_k+{\mathrm{\ϵ}}_k^h...\left(34\right)$

函数“g”包括速度矢量“V”和视线矢量“l”。视线矢量“l”将位置矢量“P”作为变量。因此，函数“g”是将位置矢量“P”和速度矢量“V”作为变量的函数。

如果用函数“g(P，V)”的位置矢量“P”、速度矢量“V”、钟差“b”以及时钟漂移“d”求偏导，则分别得到下式(35)至(38)。

[式35]

$\frac{{\∂g}_k(P,V)}{\∂P}=-\frac{V_k-V}{\left|\right|P_k-P\left|\right|}+\frac{1}{{\left|\right|P_k-P\left|\right|}^3}\left(\begin{array}{c}(v_{k,x}-v_x){(x_k-x)}^2\\ (v_{k,y}-v_y){(y_k-y)}^2\\ (v_{k,z}-v_z){(z_k-z)}^2\end{array}\right)...\left(35\right)$

[式36]

$\frac{{\∂g}_k(P,V)}{\∂V}=-l_k...\left(36\right)$

[式37]

$\frac{{\∂g}_k(P,V)}{\∂b}=0...\left(37\right)$

[式38]

$\frac{{\∂g}_k(P,V)}{\∂d}=1...\left(38\right)$

在此，为了方便起见，将就已捕获的K个卫星总结了函数“g(P，V)”的函数“G(P，V)”称为观测函数，定义为下式(39)。

[式39]

$G(P,V)=\left[\begin{array}{c}{g}_1(P,V)\\ g_2(P,V)\\ .\\ .\\ .\\ g_K(P,V)\end{array}\right]...\left(39\right)$

这种情况下，下式(40)的观测方程式成立。

[式40]

$(-G(P_0,b_0,V_0,d_0))=\left(\frac{\∂G}{\∂(P,b,V,d)}{|}_{(P_0,b_0,V_0,d_0)}\right)\left[\begin{array}{c}\mathrm{\δP}\\ \mathrm{\δb}\\ \mathrm{\δV}\\ \mathrm{\δd}\end{array}\right]$

$\→-G=\mathrm{Y\δX}...\left(40\right)$

式(40)的左边是通过位置矢量的初始值“P₀”、钟差的初始值“b₀”、速度矢量的初始值“V₀”以及时钟漂移的初始值“d₀”确定的观测函数的值“G(P₀，b₀，V₀，d₀)”。

另外，在式(40)的右边，将矢量“(δP，δb，δV，δd)^T”定义为表示状态矢量“X”的变化的状态变化量矢量“δX”。另外，为了方便起见，将用观测函数“G”的各变量的偏导定义的矩阵定义为偏导矩阵“Y”。状态变化量矢量“δX”与由各变量的初始值决定的偏导矩阵“Y”的值相乘。

这样，从式(33)导出将状态变化量矢量“δX＝(δP，δb，δV，δd)^T”作为未知数的式(40)的观测方程式。利用该观测方程式与上述实施方式同样地求出状态变化量矢量“δX”的最大似然估计值。与上述实施方式的不同点在于，作为位置计算运算基础的观测方式从式(3)置换为式(40)。

这种情况下，定义将状态变化量矢量“δX＝(δP，δb，δV，δd)^T”中包含的计算参数“(δx，δy，δz，δb，δv_x，δv_y，δv_z，δd)”的全部或一部分作为随机变量的概率分布模型。除了位置信息，将钟差、速度信息以及时钟漂移的变化作为随机变量处理。

将8个计算参数全部作为随机变量的情况下，将式(8)的多维状态变化正态分布模型定义为8维状态变化正态分布模型“p_δX ⁽⁸⁾(δX)”。另外，如果将8个计算参数中的例如除了钟差“b”和时钟漂移“d”以外的6个计算参数作为随机变量，则将式(8)的多维状态变化正态分布模型定义为6维状态变化正态分布模型“p_δX ⁽⁶⁾(δX)”。

由于在式(40)中有8个未知数，因此为了求解，必需捕获至少8个卫星。因此，实用的是想办法减少未知数。例如，将速度矢量“V”设为固定值(例如零)，将位置矢量“P”、钟差“b”以及时钟漂移“d”设为未知数。这种情况下，式(40)的观测方程式可以换成下式(41)。

[式41]

$(-G(P_0,b_0,V_0,d_0))=\left(\frac{\∂G}{\∂(P,b,d)}{|}_{(P_0,b_0,V_0,d_0)}\right)\left[\begin{array}{c}\mathrm{\δP}\\ \mathrm{\δb}\\ \mathrm{\δd}\end{array}\right]...\left(41\right)$

在式(41)中，将表示状态矢量“X”的变化的状态变化量矢量定义为“δX＝(δP，δb，δd)^T”。通过将速度矢量“V＝(v_x，v_y，v_z)”设为固定值，可以将未知数从8个减少到5个。通过这样，如果能至少捕获5个卫星就可以求出式(41)的解。

图19是示出上述实施例中的处理部21进行的多普勒定位处理的流程的流程图。该多普勒定位处理是可在图14的基带处理的步骤A1中更换成位置计算处理来使用的处理。另外，该多普勒定位处理是基于式(40)的观测方程式的处理，是对应于将状态变化量矢量“δX”的所有计算参数作为随机变量的图15的第一位置计算处理的处理。

首先，处理部21设定状态矢量“X”的初始值(步骤F1)。即，对位置矢量“P”、钟差“b”、速度矢量“V”以及时钟漂移“d”设定预定的初始值。

然后，处理部21进行捕获目标卫星选择处理(步骤F3)。并且，处理部21进行状态变化量估计处理(步骤F5至F19)。具体地，处理部21就每个捕获目标卫星执行循环C(步骤F5至F13)。

在循环C中，处理部21从该捕获目标卫星接收GPS卫星信号并计算频率漂移“h”(步骤F7)。频率漂移“h”是作为通过频率搜索获得的接收信号的频率(接收频率)与预定的载波的频率(预定载波频率)的差而计算出的。该频率漂移“h”成为多普勒定位中的观测量。

然后，处理部21计算位置计算装置与该捕获目标卫星的相对速度矢量“V_R”(步骤F9)。利用位置计算装置的最新速度矢量“V”和该捕获目标卫星的最新速度矢量“V_k”计算相对速度矢量“V_R”。初次定位的情况下，将速度矢量的初始值“V₀”作为最新速度矢量，第二次以后定位时，将在上次定位时计算的速度矢量作为最新的速度矢量进行计算。

然后，处理部21计算从位置计算装置向着该捕获目标卫星的视线矢量“l”(步骤F11)。利用位置计算装置的最新的位置矢量“P”和该捕获目标卫星的最新的位置矢量“P_k”计算视线矢量“l”。初次定位的情况下，将位置矢量的初始值“P₀”作为最新的位置矢量，第二次以后定位时，将在上次定位时计算的位置矢量作为最新的位置矢量进行计算。

然后，处理部21转向对下一个捕获目标卫星的处理。就所有的捕获目标卫星进行了步骤F7至F11之后，处理部21结束循环C(步骤F13)。

在循环C之后，处理部21使用最新的位置矢量“P”、钟差“b”、速度矢量“V”以及时钟漂移“d”计算式(40)左边的观测函数“G”的值(步骤F15)。另外，处理部21使用最新的位置矢量“P”、钟差“b”、速度矢量“V”以及时钟漂移“d”计算式(40)右边的偏导矩阵“Y”的值(步骤F17)。

然后，处理部21根据式(40)的观测方程式，使用观测误差正态分布模型和状态变化正态分布模型求出状态变化量“δX”的最大似然估计解(步骤F19)。观测误差正态分布模型是将作为观测量的频率漂移“h”中包含的观测误差“ε^h”作为随机变量、将其分布假定为正态分布的模型。

进行了步骤F19之后，处理部21结束状态变化量估计处理。之后的步骤与第一位置计算处理相同。

另外，虽然此处省略了流程图的图示和说明，但对于与图16至图18的第二至第四位置计算处理相对应的多普勒定位处理，也可以同样地构造。

5-2.随机变量

在上述的实施方式中，定义了将位置计算运算的运算结果的变化作为随机变量的概率分布模型，但可以定义将位置计算运算的运算结果其本身作为随机变量的概率分布模型并进行位置计算运算。

重点考虑位置矢量“P”。在位置矢量的初始值“P₀”上加上位置变化量矢量“δP”得到位置矢量“P”。由于位置矢量的初始值“P₀”是固定值，因此将位置变化量矢量“δP＝P-P₀”作为随机变量处理和将位置矢量“P”其本身作为随机变量处理是等价的。因此，也可以用位置矢量“P”其本身替代位置变化量矢量“δP”作为随机变量进行位置计算运算。钟差“b”、速度矢量“V”以及时钟漂移“d”也同样。

这种情况下，例如定义将GPS接收器的状态矢量“X”其本身作为随机变量的状态正态分布模型。然后，利用该状态正态分布模型与上述的实施方式同样地进行位置计算运算。更新状态正态分布模型时，利用过去计算的位置矢量“P”或钟差“b”更新状态正态分布模型的期望值和标准差。

5-3.状态变化正态分布模型的更新

在上述的实施方式中，基于一个时刻前的位置计算运算的结果对状态变化正态分布模型进行了更新，但这只不过是一个例子。例如也可以基于两个时刻前或三个时刻前的位置计算运算的结果来更新状态变化正态分布模型。

另外，也可以对过去预定期间(例如5个时刻或10个时刻)的位置计算运算的结果进行平均处理，基于该结果更新状态变化正态分布模型。平均处理可以使用简单的算术平均或几何平均，也可以使用加权平均。如果使用加权平均，对离现在的运算时刻越近的时刻的位置计算运算的结果设置越高的权重以进行加权平均是有效的。

另外，也可以使用可从外部获取的参考信息来更新状态变化正态分布模型。在汽车等移动体上设置本发明的位置计算装置的情况下，在能够从车速检测系统(例如车速脉冲)或磁传感器等方位传感器获取移动体的速度或移动方向(即速度矢量)作为参考信息的情况下，也可以构成为利用该参考信息更新状态变化正态分布模型。

具体地，在最新的移动体的位置矢量“P”中加进速度矢量“V”，计算现在的运算时刻的移动体的位置矢量“P”。由于可以这样计算位置矢量“P”，因此可以计算每个运算时刻的位置矢量的变化量“δP”。因此，可以利用速度矢量“V”计算位置变化量矢量“δP”，并利用其更新状态变化量期望值矢量“μ_δX”。

另外，同样也可以更新状态变化正态分布模型的标准差。如果知道移动体的行进方向，由于对于该行进方向以外的方向几乎没有速度分量，因此可以认为对于该方向位置和速度的误差接近零。因此，以使状态变化量协方差矩阵“S_δX”中的与移动体的行进方向以外的方向对应的协方差的分量为较小的值(例如零)的方式来更新状态变化量协方差矩阵“S_δX”是有效的。

5-4.非随机变量

在上述的实施方式中，以将可通过位置计算运算计算的计算参数中的位置矢量“P”的各分量作为随机变量、将钟差“b”作为非随机变量的情况为例进行了说明。但是这只不过是一个例子，可以适当地选择作为随机变量和非随机变量的计算参数。

例如，也可以将三维的位置分量中的一部分分量的变化作为随机变量。具体地，也可以将位置变化量矢量“δP”中的高度方向的位置变化量“δz”作为非随机变量，利用接收信号对“δz”进行预定的收敛运算来求出近似解。另外，当然也可以将高度方向的位置变化量“δz”和钟差变化量“δb”这两种计算参数作为非随机变量，对“δz”以及“δb”进行收敛运算求出近似解。

5-5.最大似然估计算法

在上述的实施方式中，作为用于最大似然估计的算法，以使用EM算法为例进行了说明，当然也可以使用其它的算法。也可以使用例如称为登山法的梯度法。梯度法是对于优化问题的一种算法，是利用预定的评价函数的梯度进行参数优化的方法。

5-6.电子设备

在上述的实施例中，作为具有GPS接收器的电子设备，以移动电话为例进行了说明，但可以使用本发明的电子设备不仅限于此。例如，汽车导航系统、便携式导航设备、个人电脑、PDA(个人数字助理)、手表等其它电子设备同样也可以使用本发明。

5-7.卫星定位系统

另外，在上述的实施方式中，作为卫星定位系统，以GPS为例进行了说明，当然也可以是WAAS(广域增强系统)、QZSS(准天顶卫星系统)、GLONASS(全球导航卫星系统)、GALILEO等其它卫星定位系统。

5-8.处理的主体

在上述的实施例中，就设置于GPS接收部内部的处理部进行位置计算处理进行了说明。即，在上述的实施例中，就GPS接收部(基带处理电路部)具有用作位置计算装置的功能进行了说明。但作为电子设备的处理器的主机处理部也可以作为位置计算的主体。

这种情况下，例如，基带电路部的处理部捕获GPS卫星信号，计算并获取测量信息并输出至主机处理部。并且，主机处理部利用从处理部输入的测量信息，基于上述原理进行位置计算运算。这种情况下，电子设备(主机处理部)具有用作位置计算装置的功能。

符号说明

1移动电话，5GPS天线，10GPS接收部，11RF接收电路部，20基带处理电路部，21处理部，23存储部，30主机处理部，40操作部，50显示部，60移动电话用天线，70移动电话用无线通信电路部，80存储部，90时钟部。

Claims (6)

1.一种位置计算方法，包括：

接收定位信号；以及

利用所接收的定位信号进行位置计算运算，所述位置计算运算是基于至少将过去的运算结果的变化作为随机变量的给定的概率分布模型的运算，

所述位置计算运算是基于至少将所述过去的运算结果的变化作为随机变量的正态分布模型的运算，并且

所述位置计算方法还包括：基于所述过去的运算结果，对所述正态分布模型进行更新，

基于所接收的定位信号测量伪距，

其中，所述位置计算运算是基于将所述过去的运算结果的变化作为随机变量的正态分布模型和将所述伪距中包含的误差作为随机变量的概率分布模型的运算。

2.根据权利要求1所述的位置计算方法，还包括：

基于所述过去的运算结果，对所述变化的标准进行更新；以及利用所述标准计算所述变化。

3.根据权利要求1所述的位置计算方法，

其中，将所述伪距中包含的误差作为随机变量的概率分布模型是正态分布模型或正态混合分布模型。

4.根据权利要求1所述的位置计算方法，

其中，所述位置计算运算是计算位置信息、钟差、速度信息以及时钟漂移中的至少任意一个的运算，

其中，所述正态分布模型是将所述过去的运算结果的各元素的变化作为随机变量的模型。

5.根据权利要求1所述的位置计算方法，还包括：通过利用所接收的信号进行预定的收敛运算来计算计算参数中的预定参数的值，所述计算参数能够通过使用所述位置计算运算来计算，

其中，所述位置计算运算是利用过去通过使用所述收敛运算而算出的所述预定参数的值进行的运算。

6.一种位置计算装置，包括：

接收部，接收定位信号；以及

位置计算部，利用所述接收部接收的定位信号进行位置计算运算，所述位置计算运算是基于至少将过去的运算结果的变化作为随机变量的给定的概率分布模型的运算，

所述位置计算运算是基于至少将所述过去的运算结果的变化作为随机变量的正态分布模型的运算，并且

所述位置计算装置基于所述过去的运算结果，对所述正态分布模型进行更新，

所述位置计算装置还基于所接收的定位信号测量伪距，

其中，所述位置计算运算是基于将所述过去的运算结果的变化作为随机变量的正态分布模型和将所述伪距中包含的误差作为随机变量的概率分布模型的运算。