专利摘要

本发明公开了一种基于广义Pareto分布杂波统计建模的雷达目标检测方法，目的是提高雷达目标检测性能。技术方案是，将对数最大似然估计与粒子群优化算法相结合，通过对数最大似然估计准则来设计代价函数，并利用粒子群优化算法来最小化代价函数以获取最优的参数估计结果；利用最优的参数估计结果确定检测门限，完成雷达目标检测。采用本发明可确保待估计参数在最大似然意义下解的最优性，可提高参数估计的精度和鲁棒性，并提高检测门限的准确性，确保雷达目标检测结果的可靠性。

权利要求

1.一种基于广义Pareto分布杂波统计建模的雷达目标检测方法，其特征在于包括以下步骤：

第一步，定义并初始化粒子群，方法是：

1.1定义粒子群G：

粒子群为由一群粒子即待估计参数组成的集合，粒子群定义为

G＝{p_i＝(σ_i,k_i),v_i＝(δσ_i,δk_i)；i＝1,2,…,I} (1)

式中，I为G中的粒子个数，p_i为G中第i个粒子的位置特征，σ_i为G中第i个粒子的尺度参数，k_i表示G中第i个粒子的形状参数；v_i为G中第i个粒子的速度特征，δσ_i为G中第i个粒子的尺度参数σ_i变化量，δk_i为G中第i个粒子的形状参数k_i变化量；

1.2初始化粒子位置特征：

将从杂波区观测得到的N个杂波数据样本记为z₁,...,z_n,...,z_N，1≤n≤N，n为整数，则N个杂波数据样本的均值和方差分别为 和 在没有先验信息的条件下，用 和 分别表示参数σ和k的粗估计值，k表示形状参数，σ表示尺度参数，σ＞0，将第i个粒子0时刻位置特征初始化为

式中，I为正整数；上标0表示参数的初始状态即0时刻， 为G中第i个粒子0时刻的位置特征， 为G中第i个粒子0时刻的尺度参数， 为G中第i个粒子0时刻的形状参数； 表示均匀分布函数，区间下限为a，上限为b； 表示从均匀分布函数 中随机取数，即在区间[a,b]之间随机取一个数。 表示 是从均匀分布函数 中随机取数，即在区间 之间随机取一个数， 为下限， 为上限； 表示 是从均匀分布函数 中随机取数， 为下限， 为上限；

1.3初始化粒子速度特征：

为G中第i个粒子0时刻的速度特征， 为G中第i个粒子0时刻的尺度参数变化量， 为G中第i个粒子0时刻的形状参数变化量， 表示 是从均匀分布函数 中取的数， 表示 是从均匀分布函数 中取的数；

第二步，令迭代次数变量t＝0；

第三步，求取第t次迭代时粒子的代价函数，方法是：

3.1将广义Pareto分布模型的复杂指数分布转化为公式(8)所示的对数形式的似然函数，k＜0，即

式中，ln(·)表示以自然常数e为底的对数函数，L(z|k,σ)表示包含未知分布参数k和σ的杂波信号幅度z的似然函数；

参数k和σ的最大似然估计通过求解公式(8)的偏导数得到，此时有

3.2根据公式(9)和(10)构建公式(11)所示的代价函数：

式中，|·|为取绝对值符号，通过该代价函数的最小化使T(σ,k)不断向0逼近，这一过程等效于使公式(9)和公式(10)的解达到最优；

3.3通过公式(11)求取G中I个粒子对应的适应性值，I个粒子经第t次迭代得到的I个粒子第t组适应性值，表示为

第四步，计算第t次迭代中粒子群的个体极值pbest^t与全局极值gbest^t，方法是：

4.1根据 找到第t次迭代中与最小粒子适应性值对应的粒子，即

式中， 表示先找到 中最小的值，并找到这个最小值对应的粒子，表示为参数对(σ′,k′)；

4.2根据 找出0～t次所有迭代过程中使得适应值最小的粒子，作为前t次迭代全局极值gbest^t，即

(σ^*,k^*)为前t次最小适应值对应的粒子；

第五步，令t＝t+1；

第六步，根据获取的个体极值pbest^t与全局极值gbest^t，更新G中I个粒子的位置特征与速度特征；

第七步，判断t是否等于最大迭代次数t_max，若满足，将全局极值gbest^t，也即(σ^*,k^*)作为参数对(σ,k)的最终估计结果，执行第八步；否则，转第三步；

第八步，利用(σ^*,k^*)进行雷达目标检测，方法是：

8.1利用(σ^*,k^*)重构与广义Pareto分布模型相对应的概率密度函数 即

式中， 表示由(σ^*,k^*)重构得到的关于自变量z的近似函数；

8.2给定目标检测虚警率P_f，根据公式(15)求检测门限th

8.3目标检测，方法是：

8.3.1通过雷达实时观测，获得J个观测数据，J≥1，令j＝1；

8.3.2判断雷达实际获取的第j个观测数据y_j是否有目标，方法是：

若y_j≥th，输出“第j个观测数据有目标”的结论；若y_j＜th，输出“第j个观测数据无目标”的结论；

8.3.3判定j是否小于J，若满足，令j＝j+1，转8.3.2；否则，表示雷达实时观测数据处理结束，完成目标检测。

2.如权利要求1所述的基于广义Pareto分布杂波统计建模的雷达目标检测方法，其特征在于所述I取200～500。

3.如权利要求1所述的基于广义Pareto分布杂波统计建模的雷达目标检测方法，其特征在于3.3步所述求取G中I个粒子对应的适应性值的方法是：依次将第t次迭代中G中I个粒子的位置特征代入公式(11)，得到I个粒子第t次迭代的代价函数值，作为I个粒子第t次迭代的适应性值；对于G中第i个粒子，具体做法是将其第t次迭代中的位置特征 即 参数对代入公式(11)，计算 并令第i个粒子第t次迭代的粒子适应性值

4.如权利要求1所述的基于广义Pareto分布杂波统计建模的雷达目标检测方法，其特征在于第六步所述更新G中I个粒子的位置特征与速度特征的方法是：

6.1令i＝1；

6.2第i个粒子的位置特征和速度特征按公式(17)和公式(16)进行更新：

其中，w^t-1＝0.9-0.5·(t-1)/t_max，为t-1次时的惯性权因子，t_max为最大迭代次数，为正整数；rand表示[0,1]之间的均匀分布随机数；c₁与c₂为学习因子，为正整数；

6.3令i＝i+1；

6.4判定i≤I是否成立，若成立，转6.2；否则，表示已更新完G中I个粒子的位置特征和速度特征，结束。

5.如权利要求4所述的基于广义Pareto分布杂波统计建模的雷达目标检测方法，其特征在于所述t_max取100～200；c₁与c₂均取为2。

6.如权利要求1所述的基于广义Pareto分布杂波统计建模的雷达目标检测方法，其特征在于所述P_f取10^-4～10^-2。

7.如权利要求1所述的基于广义Pareto分布杂波统计建模的雷达目标检测方法，其特征在于所述对公式(15)求检测门限th通过牛顿-科茨即Newton-Cotes数值积分公式来解算。

说明书

技术领域

本发明属雷达目标检测领域，涉及一种利用粒子群优化算法对广义Pareto分布进行精确统计建模的雷达目标检测方法。

背景技术

杂波数据统计建模(即利用合适的分布或函数对杂波数据直方图进行精确拟合)是雷达目标检测领域面临的一个重要问题，选取合适的杂波数据拟合分布模型(或称拟合函数)并准确获取模型参数是目标检测算法设计(如检测门限、虚警概率、检测概率计算等)的重要基础。若模型参数估计存在较大的误差，就无法得到准确的检测门限，从而导致检测概率下降和虚警概率上升。

常用的杂波数据拟合分布模型有Rayleigh分布模型、Weibull分布模型，以及K分布模型、Gamma分布模型等。这些分布模型有其自身的适应范围，如Rayleigh分布模型常用于描述低分辨均匀杂波数据，K分布模型则常用于描述低入射余角条件下的高分辨海杂波数据。

广义Pareto分布模型是一种重要的统计分布模型，是根据意大利经济学家Vilfredo Pareto的名字命名的。该模型最初主要应用于经济学、物理学、水文学与地震学等领域，现已逐步应用于高分辨雷达杂波数据的统计建模，如文献1：G.V.Weinberg，“Assessing Pareto fit to high-resolution high-grazing-angle sea clutter”，Electronics Letters，2011，47(8)：516-517(G.V.Weinberg于2011在《电子快报》第47卷第8期发表的论文，中文题目为“帕累托分布对高入射余角海杂波的拟合能力评估”)中的研究结果表明，广义Pareto分布模型对大入射余角条件下的高分辨、大拖尾海杂波数据具有良好的模型匹配和统计分布拟合能力。由此可见，广义Pareto分布模型参数估计是解决大入射余角特殊应用条件下、海杂波背景中目标检测问题涉及的一项重要技术。

广义Pareto分布模型可用下式表示为

式中，f(·)表示概率密度函数，z为自变量(代表杂波信号幅度，且其值大于零)，k表示形状参数，σ(σ＞0)表示尺度参数，exp(·)表示以自然常数e为底的指数函数，f(z|k,σ)表示包含未知分布参数(k和σ)的杂波信号幅度z的概率密度函数。

显然，当k＝0时，广义Pareto分布就退化为指数分布模型，参数σ的求解比较简单，不在本发明的讨论范围。

在k＜0的情况下，广义Pareto分布模型的参数(即k、σ)通常利用杂波信号幅度z(也称杂波数据样本)的r阶中心矩E(z^r)(E(·)表示取数学期望，r≥1且为整数)进行估计，如目前常见的矩估计法(Method of Moments，MoM)、概率加权矩方法(Probability Weighted Moments Method，PWMM)、似然矩方法(Likelihood Moment Method，LMM)等。文献2：P.de Zea Bermudez，Samuel Kotz，“Parameter estimation of the generalized Pareto distribution Parts I&II”，Stat.Plann.Inference 140，2010：1353-1388(P.deZea Bermudez等人于2010年在《统计规划与推理期刊》第140卷发表的论文，中文题目为“广义帕累托分布的模型参数估计”)中对若干种广义Pareto分布模型的参数估计方法进行了总结、分析与对比，并指出最大似然估计(Maximum Likelihood Estimation，MLE)是最有效的估计方法。

事实上，当雷达杂波呈现大拖尾分布时，其实测数据统计原点矩与相应阶数的杂波模型理论原点矩之间会出现一定的偏差，且该偏差随矩阶数(r)的增大和杂波数据样本数量的减小而增大，这将对基于统计矩的估计方法带来不利的影响。最大似然估计具有很高的估计精度，但是当雷达杂波分布模型(如广义Pareto分布模型，也即概率密度函数表达式)比较复杂时，该方法很难得到参数估计的解析表达式；在这种情况下，虽然可以通过如文献3：Tjalling J.Ypma，“Historical development of the Newton-Raphson method”，SIAM Review，1995，37(4)：531-551(Tjalling J.Ypma于1995年在《美国工业与应用数学学会评论》期刊第37卷第4期发表的论文，中文题目为牛顿-拉夫森方法的历史发展)中的Newton-Raphson数值算法进行求解，但这种数值方法容易陷入局部极值而无法得到最优参数估计。

粒子群优化(Particle Swarm Optimization，PSO)是一种重要的智能优化算法，具有计算简单、收敛速度快、稳定性好等特点，且具有全局寻优能力。通过PSO算法，可望解决复杂模型的参数高精度估计问题，从而提高雷达目标检测性能，但目前尚未见到有公开文献涉及如何运用PSO算法进行Pareto杂波模型参数高精度估计和高性能目标检测。

发明内容

本发明要解决的技术问题是：提供一种基于广义Pareto分布杂波统计建模的雷达目标检测方法，提高雷达目标检测性能。

技术方案是，将对数最大似然(Logarithmic Maximum Likelihood，LML)估计与粒子群优化(PSO)算法结合起来，通过对数最大似然估计准则来设计代价函数，并利用粒子群优化算法来最小化代价函数以获取最优的参数估计结果；利用最优的参数估计结果确定检测门限，完成雷达目标检测。

具体步骤如下：

第一步，定义并初始化粒子群。方法是：

1.1定义粒子群G：

粒子群为由一群粒子组成的集合，此处的粒子即为待估计参数，粒子群可定义为

G＝{p_i＝(σ_i,k_i),v_i＝(δσ_i,δk_i)；i＝1,2,…,I} (2)

式中，I为G中的粒子个数，p_i为G中第i个粒子的位置特征，σ_i为G中第i个粒子的尺度参数，k_i表示G中第i个粒子的形状参数；v_i为G中第i个粒子的速度特征，δσ_i为G中第i个粒子的尺度参数σ_i变化量，δk_i为G中第i个粒子的形状参数k_i变化量。

1.2初始化粒子位置特征：

将从杂波区观测得到的N个杂波数据样本记为z₁,...,z_n,...,z_N(1≤n≤N，n为整数)，则N个杂波数据样本的均值和方差分别为 和 在没有先验信息的条件下，可用 和 分别表示参数σ和k的粗估计值，由此将第i个粒子0时刻位置特征初始化为

式中，I为正整数，I的大小由经验值给定，一般取200～500为宜；上标0表示参数的初始状态(即0时刻)， 为G中第i个粒子0时刻的位置特征， 为G中第i个粒子0时刻的尺度参数， 为G中第i个粒子0时刻的形状参数； 表示均匀分布函数，区间下限为a，上限为b； 表示从均匀分布函数 中随机取数，即在区间[a,b]之间随机取一个数。 表示 是从均匀分布函数 中随机取数(即在区间 之间随机取一个数)， 为下限， 为上限； 表示 是从均匀分布函数 中随机取数， 为下限， 为上限。

1.3初始化粒子速度特征：

为G中第i个粒子0时刻的速度特征， 为G中第i个粒子0时刻的尺度参数变化量， 为G中第i个粒子0时刻的形状参数变化量。 表示 是从均匀分布函数 中取的数， 表示 是从均匀分布函数 中取的数。

第二步，令迭代次数变量t＝0。

第三步，求取第t次迭代时粒子的代价函数。方法是：

3.1将广义Pareto分布模型的复杂指数分布(k＜0时的情形)转化为公式(9)所示的对数形式的似然函数，即

式中，ln(·)表示以自然常数e为底的对数函数，L(z|k,σ)表示包含未知分布参数(k和σ)的杂波信号幅度z的似然函数。

由此可见，参数k和σ的最大似然估计可通过求解公式(9)的偏导数得到，此时有

很明显，通过公式(10)和公式(11)两式无法求得关于参数k和σ的显式表达式，也就是说，无法通过将公式(10)和公式(11)两式得到的偏导数结果置零求出参数k和σ的解，本发明转而通过粒子群优化算法来求数值解，需构建代价函数。

3.2根据公式(10)和(11)构建公式(12)所示的代价函数

式中，|·|为取绝对值符号，通过该代价函数的最小化可使T(σ,k)不断向0逼近，这一过程等效于使公式(10)和公式(11)的解达到最优。

3.3通过公式(12)求取G中I个粒子对应的适应性值，方法是

依次将第t次迭代中G中I个粒子的位置特征代入公式(12)，得到I个粒子第t次迭代的代价函数值，作为I个粒子第t次迭代的适应性值。对于G中第i个粒子，具体做法是将其第t次迭代中的位置特征 (即 参数对)代入公式(12)，计算 并令第i个粒子第t次迭代的粒子适应性值

I个粒子经第t次迭代得到的I个粒子第t组适应性值，表示为

第四步，计算第t次迭代中粒子群的个体极值pbest^t与全局极值gbest^t。方法是：

4.1根据 找到第t次迭代中与最小粒子适应性值对应的粒子，即

式中， 表示先找到 中最小的值，并找到这个最小值对应的粒子(表示为参数对(σ′,k′))。(假设该最小值的序号为i，则找到的粒子为 即

4.2根据 找出0～t次所有迭代过程中使得适应值最小的粒子，作为前t次迭代全局极值gbest^t，即

式中，(σ^*,k^*)表示前t次迭代最小适应值对应的粒子。

第五步，令t＝t+1。

第六步，根据获取的个体极值pbest^t与全局极值gbest^t，更新G中I个粒子的位置特征与速度特征，方法是：

6.1令i＝1。

6.2第i个粒子的位置特征和速度特征按公式(15)和公式(16)进行更新：

其中，w^t-1＝0.9-0.5·(t-1)/t_max，为t-1次时的惯性权因子，t_max为最大迭代次数(为正整数，通常取100～200)；rand表示[0,1]之间的均匀分布随机数；c₁与c₂为学习因子(或加速常数)，为正实数，通常均取为2。

6.3令i＝i+1。

6.4判定i≤I是否成立，若成立，转6.2；否则，表示已更新完G中I个粒子的位置特征和速度特征，执行第七步。

第七步，判断t是否等于最大迭代次数t_max，若满足，则将全局极值gbest^t，也即(σ^*,k^*)作为参数对(σ,k)的最终估计结果，执行第八步；否则，转第三步。

第八步，利用(σ^*,k^*)进行雷达目标检测，方法是：

8.1利用(σ^*,k^*)重构与广义Pareto分布模型相对应的概率密度函数 即

式中， 表示由(σ^*,k^*)重构得到的关于自变量z的近似函数。

8.2给定目标检测虚警率P_f(P_f通常取10^-4～10^-2)，根据公式(18)求检测门限th

式中，积分算式可通过牛顿-科茨(Newton-Cotes)数值积分公式来解算，具体解法可参考文献3：叶其孝，沈永欢等，“实用数学手册(第2版)”，科学出版社，2006年，719～720。

8.3目标检测，方法是：

8.3.1通过雷达实时观测，获得J个观测数据，J≥1，令j＝1；

8.3.2判断雷达实际获取的第j个观测数据y_j是否有目标，方法是：

若y_j≥th，输出“第j个观测数据有目标”的结论；若y_j＜th，输出“第j个观测数据无目标”的结论；

8.3.3判定j是否小于J，若满足，令j＝j+1，转8.3.2；否则，表示雷达实时观测数据处理结束，完成目标检测。

本发明的有益效果：

1.本发明第三步通过对数最大似然函数构建Pareto分布模型参数的估计式(公式(9))，确保待估计参数在最大似然意义下解的最优性；

2.本发明第四步利用粒子群优化算法具有全局寻优的能力，克服传统数值方法(如Newton-Raphson等)在解决模型参数估计问题时存在的对初值敏感、易于收敛于局部极值的缺陷，可提高参数估计的精度和鲁棒性，确保雷达目标检测结果的可靠性。

3.本发明将PSO算法与广义Pareto分布模型参数的最大似然参数估计方法结合起来，提高了检测门限th的准确性，进而提高了雷达目标的检测性能。

附图说明

图1是本发明总体流程图。

具体实施方式

下面结合附图和实验对本发明的实施方式进行说明。

图1给出了本发明总体流程图，下面首先以杂波样本数N＝100为例来说明具体的实验过程。在此基础上，通过不同的实验条件设置(如改变杂波样本数N和模型形状参数k等)，以考察本发明在不同实验条件下的检测性能。

第一步，定义并初始化粒子群。从杂波区得到N＝100个杂波数据样本，设定粒子群大小I＝500，粒子群定义为G＝{p_i＝(σ_i,k_i),v_i＝(δσ_i,δk_i)；i＝1,2,…,I}粒子群的初始位置分量 分别在 和 区间随机取样，初始速度分量 统一在 之间随机取样，其中i＝1,2,…,I；最大迭代次数为t_max＝200。

第二步，令t＝0。

第三步，计算第t(t≥0)次迭代中的粒子适应性值。根据公式(12)中的代价函数T(σ,k)求取G中I个粒子对应的适应性值，得到I个粒子第t次迭代中所有I个粒子的适应性值，其中第i个粒子的适应性值 I个粒子经第t次迭代得到的I个粒子第t组适应性值，表示为

第四步，分别利用公式(13)和公式(14)寻找t次迭代过程中粒子群的个体极值pbest^t与全局极值gbest^t，gbest^t＝(σ^*,k^*)；

第五步，令t＝t+1。

第六步，按照公式(15)和公式(16)更新I个粒子的速度和位置分量。

第七步，判断t是否达到最大迭代次数t_max。若满足，将(σ^*,k^*)作为参数估计结果输出，并转第八步；否则返回第三步。

第八步，雷达目标检测。首先根据公式(17)构建概率密度函数 在此基础上，利用公式(18)确定检测门限th，最后通过雷达实时观测，获得J个观测数据(实验中J＝500)，对J个观测数据按步骤8.3.2～8.3.3进行有无目标的判断从而完成目标检测。

为说明本发明在不同形状参数、不同杂波样本数量条件下对广义Pareto分布模型参数的估计性能，实施例中选取了三个典型参数对(σ,k)(对应三组不同的实验)进行仿真验证。参数估计精度用平均误差(Mean Error，ME)与均方根误差(Root of Mean Square Error，RMSE)两个指标来评估，其表达式分别为：

式中， 表示第m次仿真中得到的参数υ(υ代表尺度参数或形状参数)的估计，M表示Monte Carlo仿真次数。

考虑到尺度参数σ对估计结果的影响不大，而形状参数k对广义Pareto分布影响较大，且k值越小，广义Pareto分布拖尾越大。本发明重点分析评估所述方法对形状参数的估计效果，为此设定σ＝1，k＝-0.3,-0.2,-0.1(即分别将k的值取为-0.3，-0.2，-0.1，进行三组不同的实验)，M＝500，对本发明基于粒子群优化的广义Pareto杂波模型目标检测方法(LML-PSO)进行测试，并比较其与传统方法(如MoM、PWMM、LMM等)之间的性能差异。实验在通用计算机平台上进行，利用Matlab软件实现，所得参数(σ,k)在不同杂波数据样本数(N＝25,50,100,500)、不同估计方法(MoM、PWMM、LMM和LML-PSO)条件下的仿真结果如表1～表3所示。其中，表1～表3分别对应k＝-0.3、k＝-0.2和k＝-0.1条件下的平均误差(ME)和均方根误差(RMSE)。

从表1～表3可看出，本发明LML-PSO的平均误差与均方根误差均明显小于另三种方法。

表1

表2

表3

为了进一步说明形状参数k估计结果对目标检测的影响，以k＝-0.3为例，在虚警概率P_f＝0.01的应用条件下，若不存在参数估计误差，通过公式(18)可解算出真实的目标检测门限为th＝9.898。当存在不同程度的形状参数估计误差(即ME和RMSE均不为0)时，得到相应的目标检测门限如表4所示，表中误差大小范围为-0.12～0.12(误差间隔为0.002)，即对应真实参数±40％(-0.3×40％～0.3×40％)的误差范围。当估计误差为-0.12时，目标检测门限为th＝14.023；当估计误差为0.12时，目标检测门限为th＝7.149。

表4

从表中可以看出：

1)总体上，各种方法的估计误差(或误差绝对值)随杂波数据样本数量(N)的增加而降低，说明杂波数据样本数量越大，参数估计的精度就越高，估计性能也就越好；

2)三组不同的实验中，在杂波数据样本数相同的条件下，本发明所述方法与其它三种估计方法相比，都具有最小的平均误差与均方根误差，说明其具有一致较高的参数估计精度；

3)通过对比可知，即便在杂波数据样本数目较小(如N＝50)的情况下，LML-PSO方法仍具有较高的估计精度，说明对小样本也有较好的适应性。

4)相比之下，LML-PSO方法具有最小的参数估计误差，说明由此计算得出的目标检测门限与检测门限的理论值最为接近，相应地检测性能也最佳。

综上，与传统方法相比，本发明所述方法对广义Pareto分布模型参数估计具有明显的精度性能优势，且能很好地适应模型参数的变化，可有效提高长拖尾分布杂波背景中的雷达目标检测性能。

虽然参照上述实施例详细描述了本发明，但是应该理解本发明并不限于所公开的实施例。对于本专业领域的技术人员来说，可以对其形式和细节进行各种改变。本发明涵盖了所附权利要求书的精神和范围内的各种变形。

一种基于广义Pareto分布杂波统计建模的雷达目标检测方法专利购买费用说明

Q：办理专利转让的流程及所需资料

A：专利权人变更需要办理著录项目变更手续，有代理机构的，变更手续应当由代理机构办理。

1：专利变更应当使用专利局统一制作的“著录项目变更申报书”提出。

2：按规定缴纳著录项目变更手续费。

3：同时提交相关证明文件原件。

4：专利权转移的，变更后的专利权人委托新专利代理机构的，应当提交变更后的全体专利申请人签字或者盖章的委托书。

Q:专利著录项目变更费用如何缴交

A：（1）直接到国家知识产权局受理大厅收费窗口缴纳,（2）通过代办处缴纳，（3）通过邮局或者银行汇款,更多缴纳方式

Q：专利转让变更，多久能出结果

A：著录项目变更请求书递交后，一般1-2个月左右就会收到通知，国家知识产权局会下达《转让手续合格通知书》。