静校正学习笔记-1

提示：文章写完后，目录可以自动生成，如何生成可参考右边的帮助文档

---

零、序

与许多特殊的地球物理课题不同，静校正是数据采集、资料处理、和资料解释的基础。

---

一、绪论

1.定义

1. 静校正（Static correction）：用来校正地形起伏或近地表层 ----注重反射法
2. 静校时移（Static shift）：用来对数据进行近地表电阻率不均匀性校正或地形校正
3. 静自然电位（Static SP）：对厚的净含水沙层与页岩之间的电压的一种度量
4. 静校正定义：对地震资料所做的校正，用于补偿由高程、风化层厚度以及风化层速度产生的影响，把资料校到一个指定的基准面上。其目的通常是获得在一个平面上进行采集，且没有风化层或低速介质存在时的反射波到达时间。
5. 静校正与数据采集（基准面或野外静校正）、数据处理（剩余静校正）以及剖面解释等有密切关系，因此它不属于这三个分支中的任何一支。但它会影响每一个步骤
6. 注：负的静校正量减少反射时间

对于某一点A来说，其总的静校正量$T_A$应定义为
$$T_A=\pm t_{Aw}\pm t_{Ae}\text{\hspace{1em}(1 - 1)}$$
其中，$t_{Aw}$表示A点处穿过风化层的旅行时，$t_{Ae}$表示从风化层底到基准面的旅行时，$\pm$表示旅行时的正负。

二、剩余静校正

1. 引言

1. 基准面静校正一般是根据相当简单的近地表模型进行的，如果缺乏详细的近地表资料，所建立的近地表模型的精度就会降低，从而导致基准面静校正精度降低。

2. 剩余静校正的目的是为了校正近地表模型中一些小的误差，应用剩余静校正后的最终叠加剖面应优于只做过基准面静校正的剖面。

3. 剩余静校正不能取代基准面静校正，一般而言，做好剩余静校正的前提是基准面静校正已经做的比较好。

4. 平滑数据的静校正方法假设大多数同相轴共有的不规则变化是由近地表的变化引起的，因此通过地震道时移的静校正应该可以减小这种不规则性的影响，大多数自动剩余静校正都采用统计方法。

5. 以下为应用了基准面静校正和剩余静校正的共中心点叠加剖面，（a）为基准面剩余静校正，(b) 为剩余静校正。
   

6. 在剩余静校正时，一般在这前后都需要进行速度分析，如果二次速度分析有较大的变化，则使用二次速度分析再进行剩余静校正，基本流程：
   野外带 $⟶$ 预处理 $⟶$ 速度分析 $⟶$ 初步叠加 $⟶$ 剩余静校正 $⟶$ 速度分析 ($\Lsh$) $⟶$ 最终叠加

2. 基本方法与公式

a. 基本方程

1. 基本旅行时方程：旅行时（炮检距为$x$） = 法向入射时间 + 动校正时差 + 炮点静校正量 + 接收点静校正量 + 误差项，这也是求解剩余静校正量的重要公式。一般表示为：
   $$T_{ijk}=G_k+S_i+R_j+M_k{X}_{ij}^2+N\text{\hspace{1em}(2 - 2 - 1)}$$
   式中，$T$ 表示经过动校正后的总反射时间，$i$ 表示炮点位置，$j$ 表示接收点位置，$k$ 表示CMP点位置 $[k=0.5(i+j)]$，$G$ 代表构造项或地质项（从基准面到反射面的双程旅行时），$S$ 表示地表一致性炮点静校正量（包括剩余静校正量），$R$ 表示地表一致性接收点静校正量（包括剩余静校正量），$M$ 表示剩余动校正时差系数，$X$ 表示炮检距，$N$ 表示噪声。
2. 但由于有横向倾角等因素的影响，后又对该旅行时方程做了改进，即
   $$T_{ijkh}=G_{kh}+S_i+R_j+M_{kh}{X}_{ij}^2+D_{kh}{Y}_{ij}+N\text{\hspace{1em}(2 - 2 - 2)}$$
   式中，$D_{kh}$ 表示CMP号为$k$、层位为$h$的横向倾角系数，$Y_{ij}$表示从CMP点到实际剖面所在测线的垂直偏离距离。如果在剩余静校正之前应用了基准面校正，则（2 - 1）和（2 - 2）中的炮点/接收点的静校正量可改为剩余静校正量$\Delta S_i$。
3. fd

b. 互相关技术

互相关技术可以用来估算两道间的反射时差，式 (2 - 3) 给出了两个输入道 $G(t)$ 和 $H(t)$ 之间的归一化互相关函数 ${\varphi }_{GH}(\tau )$。
$${\varphi }_{GH}(\tau )=\frac{{\sum }_{t=t_1}^{t_2}G(t)H(t+\tau )dt}{{\left[{\sum }_{t=t_1}^{t_2}{G}^2(t)dt{\sum }_{t=t_1}^{t_2}{H}^2(t)dt\right]}^{1/2}}\text{\hspace{1em}(2 - 2 - 3)}$$
式中，$\tau$ 表示滞后时间，$t_1,t_2$ 分别表示数据的起始时间和结束时间。根据该公式，拾取最大峰值对应的时间可以获得两个输入地震道之间的时差。当两个地震道相同时，最大值为1，对应时间被称为互相关系数。
注：互相关时窗内应包含主要的反射波，并且最好避开较强的相干噪声、随机噪声以及多次波。
互相关算法也有一些难点，如信噪比较低的资料拾取精度不高，同时拾取波峰时可能会出现多个波峰，甚至变为波谷，这些都会影响互相关算法的工作。

剩余静校正算法一般利用互相关系数对时移进行加权，保证互相关系数越大对最终结果的影响越大，一般有以下加权方式：
$$W={\varphi }_{max}\text{\hspace{1em}(2 - 2 - 4.1)}$$
或
$$W=\frac{{\varphi }_{max}}{1-{\varphi }_{max}^2}\text{\hspace{1em}(2 - 2 - 4.2)}$$
当按照式（2 - 2 - 4.2）进行加权时，相关系数为0.8的加权系数是0.4的五倍，相关系数为0.9的加权系数是0.6的五倍。

3. 共地面点法

暂无

4. 共中心点方法

第3节介绍了共地面点方法，这种方法是在共炮点和共接收点道平面内进行的。但目前最常用的一种剩余静校正方法仍然是使用共中心点方法。

共中心点法是假设同一个CMP道集内各道之间的时差不包含构造项。其缺陷在于仅限于剩余静校正量小于半个主周期的情况。

a. 旅行时方程

这里的旅行时方程仍然采用式（2 - 2 - 2）。假设 $N_s$ 表示一条测线上炮点的个数，$N_t$ 表示每炮的接收道数，$R_m$ 表示沿测线移动一炮后地面接收点新增加的个数，$R_g$ 是排列内缺失的接收点数，$N_r$ 是一条测线上接收点数，$N_g$ 是一条测线上受构造项影响的总CMP点数，$N_m$ 是一条测线上收剩余动校正量影响的总CMP点数，$N_d$ 是一条测线上受横向倾角分量影响的总CMP点数，$N_u$ 是未知量或参数的个数。

假设观测值的个数$N_o$等于侧线上所有的记录道数，则有
$$N_o=N_s{N}_t\text{\hspace{1em}(2 - 4 - 1.1)}$$
$$N_r=N_t+R_g+R_m(N_s-1)\text{\hspace{1em}(2 - 4 - 1.2)}$$
$$N_g=N_t+R_g+2R_m(N_s-1)\text{\hspace{1em}(2 - 4 - 1.3)}$$
假设构造分量、动校正量和横向倾角分量的CMP个数相同，即$N_m=N_g=N_d$，则有
$$N_{u1}=N_s+N_r+N_g+N_m\text{\hspace{1em}(2 - 4 - 1.4)}$$
$$N_{u1}=3N_t+N_s(5R_m+1)-5R_m+3R_g\text{\hspace{1em}(2 - 4 - 1.5)}$$
综上所述，不难发现：用于定义静校正问题的线性联立方程组，看上去是超定的，但实际上是欠定的。

b. 旅行时方程的分解

经过旅行时分解，观测时间和计算的哥个分量之间总是存在一个剩余时差。根据旅行时方程，可将时间误差定义为：
$${ϵ}_{ijkh}=T_{ijkh}-(G_{kh}+S_i+R_j+M_{kh}{X}_{ij}^2+D_{kh}{Y}_{ij})\text{\hspace{1em}(2 - 4 - 2)}$$
误差能量为：
$$E=\sum _{(i,j,k,h)}{\left({ϵ}_{ijkh}\right)}^2=\sum _{(i,j,k,h)}{\left(T_{ijkh}-G_{kh}-S_i-R_j-M_{kh}{X}_{ij}^2-D_{kh}{Y}_{ij}\right)}^2\text{\hspace{1em}(2 - 4 - 3)}$$

1. 加权法：略
2. 广义线性反演法：将旅行时方程表示为矩阵形式进行分解
   $$\mathbf{A}\mathbf{P}=t\text{\hspace{1em}(2 - 4 - 4)}$$
   其中 $\mathbf{A}$ 表示各种观测系统参数，$\mathbf{P}$ 则表示未知量矩阵，通过加入小因子以及最小二乘解，I代表单位矩阵，可得到其解为:
   $$\mathbf{P}=({\mathbf{A}}^{\mathrm{T}}\mathbf{A}+\lambda \mathbf{I}{)}^{-1}{\mathbf{A}}^{\mathrm{T}}t\text{\hspace{1em}(2 - 4 - 5)}$$
3. 独立解和不耦合现象
4. 高斯-赛德尔迭代法
5. 零均值法

c. 观测反射时间或时移的估计

在估算反射时间之前，需要对数据进行预处理

1. 数据预处理：在作互相关前要进行数据预处理，包括噪声衰减、反褶积、动校正等。
2. 选择分析的数据：在很多情况下，对给定的CMP数据，采用多个视窗进行剩余静校正是合理的。
3. 互相关法：…
4. 模型道：…
5. 互相关叠加法：
6. 叠加能量法：一组较精确地剩余静校正量能使CMP叠加道的振幅增大。该方法的目的是在分析视窗内使叠加能量最大。若有两个地震道 $F(t)$ 和 $G(t)$，其中一道的时移为 $\Delta t$，两道之和的能量 $E(\Delta t)$为：
   $$E(\Delta t)=\sum _t[F(t-\Delta t)+G(t){]}^2\text{\hspace{1em}(2 - 4 - 6)}$$
   或者可以写成:
   $$E(\Delta t)=\sum _t[F^2(t-\Delta t)+G^2(t)]+2\sum _t[F(t-\Delta t)\times G(t)]\text{\hspace{1em}(2 - 4 - 7)}$$
7. 反攻倒算

5. 蒙特卡罗和全局优化算法

a. 蒙特卡洛

前面提到的一些算法，基本上都是用来微调基准面静校正量并且假设其值较小，这些方法都属于线性方法，往往会收敛到局部解，如果假设不成立，就需要一些非线性的方法。

想解决该问题，有一个比较笨的方法，那就是为每一个点赋予一个可能的静校正量，然后一个个的试，总有一组数据是最优的。但这样做耗时耗力，例如：若有100个点，每个点的静校正量在10以内，则每个点有21（-10 - 10）种可能，总共有 ${21}^{100}$ 种可能，使用穷举法一般不现实。

有人曾经提出了一种蒙特卡罗方法，用蒙特卡罗半随机地选取可能的剩余静校正量，从而避免陷入局部叠加能量极大值。其定义为：随机地选择可能的值兵进行反复计算的一种数学方法，其结果是统计意义下的解。

b. 模拟退火法

如果将一种熔融的物质在经过临界温度时漫漫地仔细地冷却就会形成大的晶体，如果冷却速度太快，就会形成非结晶的玻璃体或小结晶的混合物。模拟退火法就是基于这样一个事实被提出来的，该方法主要是寻找全局最小。其流程图如下：

在剩余静校正过程中，我们需要实现能量最大，也可看做负能量最小，因此可采用模拟退火法。一般来说要对整条测线进行剩余静校正分析，但大多数情况下仅仅是对怀疑存在大的剩余静校正量的部分进行分析，对发现有问题的测线段，可以针对性的进行分析。对于数据好的地段，剩余静校正量的搜寻范围可以缩小。

值得注意的是，为了节约计算机时间，常常采用大的采样间隔对剩余静校正量进行采样，这可能需要在模拟退火之后再使用互相关的方法进行微调。其伪代码如下：

/*
* J(y)：在状态y时的评价函数值
* Y(i)：表示当前状态
* Y(i+1)：表示新的状态
* r： 用于控制降温的快慢
* T： 系统的温度，系统初始应该要处于一个高温的状态
* T_min ：温度的下限，若温度T达到T_min，则停止搜索
*/
while( T > T_min )
{
　　dE = J( Y(i+1) ) - J( Y(i) ) ; 

　　if ( dE >=0 ) //表达移动后得到更优解，则总是接受移动
       Y(i+1) = Y(i) ; //接受从Y(i)到Y(i+1)的移动
　　else
　　{
       // 函数exp( dE/T )的取值范围是(0,1) ，dE/T越大，则exp( dE/T )也
      if ( exp( -dE/T ) > random( 0 , 1 ) )
      Y(i+1) = Y(i) ; //接受从Y(i)到Y(i+1)的移动
　　}
　　T = r * T ; //降温退火 ，0<r<1 。r越大，降温越慢；r越小，降温越快
　　/*
　　* 若r过大，则搜索到全局最优解的可能会较高，但搜索的过程也就较长。若r过小，则搜索的过程会很快，但最终可能会达到一个局部最优值
　　*/
　　i ++ ;
}

1. 两步法

• 首先根据用户定义的取值范围产生一个随机数，然后把这个随机数赋予某一炮点作为其临时的或推荐的剩余静校正量。
• 使用上文提出的能量方法计算与这个炮点有关的一些CMP道集的叠加能量。
• 把当前的剩余静校正量所对应的叠加能量与推荐的剩余静校正量所对应的叠加能量进行比较，并且分析它们的能量差 $\Delta E$。
• 如果 $\Delta E\le 0$，表明推荐的剩余静校正量可以接受，负责就按照概率 $P(-\Delta E)$ 接受这个推荐值。
  $$P(-\Delta E)=e^{\frac{-\Delta E}{KT}}\text{\hspace{1em}(2 - 5 - 1)}$$
  式中，T是温度，k是常数，为了方便起见，k常常为1。
• 如果接受了推荐的值，利用新的剩余静校正量可以对CMP叠加进行更新。
• 然后选择新的地面点重复上述过程，直到测线上所有的炮点和接收点都用该方法检测过。

可以看出这种方法分两步走，首先选取推荐的值和进行叠加能量计算，然后决定是否采纳推荐值。

• 再用稍小一点的温度进行下一轮迭代。
• 常用的温度函数有：
  $$T_k=\frac{C}{\log (1+k)}\text{\hspace{1em}(2 - 5 - 2)}$$
  其中k代表迭代次数，C是不依赖与k的常量，其值通常很大。
  $$T_k={\alpha }^k{T}_0\text{\hspace{1em}(2 - 5 - 3)}$$
  其中，T0为初始温度，$\alpha$ 为一个常数。
• 当经过数次迭代之后，可接受的推荐值的数量趋向于0，就可以停止迭代。

同时，还有另一种广义的模拟退火法，其主要目的是：当逐渐逼近全局最小值时，确保推荐值在错误方向被接受的概率趋向于零。这需要引进一个因子，这个因子与叠加能量 $\varphi$ 和它的全局最小值 ${\varphi }_{\mathrm{m}\mathrm{i}\mathrm{n}}$ 有关，此时的(2 - 5 - 1)可描述为：
$$p(\Delta \varphi )=e^{\frac{-\Delta \varphi (\varphi -{\varphi }_{\mathrm{m}\mathrm{i}\mathrm{n}}{)}^g}{kT}}\text{\hspace{1em}(2 - 5 - 4)}$$
式中g为任意负数。

2. 一步法或热浴法
   在这种方法中，推荐值是根据互相关导出的概率分布选取的，并且这一推荐值总是被接受。余下的其他过程都一样。这里的概率函数变为：
   $$P(S_{i_{\tau }})=\frac{e^{\frac{E(S_{i_{\tau }})}{T}}}{{\sum }_{\tau =1}^n{e}^{\frac{E(S_{i_{\tau }})}{T}}}\text{\hspace{1em}(2 - 5 - 5)}$$
   其中，n是炮点Si处可能的剩余静校正量 ($\tau$) 的个数，$E(S_{i_{\tau }})$ 为炮点Si处作为炮点剩余静校正量($\tau$) 的函数的叠加能量。为了进一步简化计算，使用${\varphi }_{FG}(\tau )$ 来代替上式中的$E(S_{i_{\tau }})$，其定义为：
   $${\varphi }_{FG}(\tau )=\frac{{\sum }_{t}F(t)G(t+\tau )dt}{\left[{\sum }_t{F}^2(t)dt{\sum }_t{G}^2(t)dt\right]}\text{\hspace{1em}(2 - 5 - 6)}$$
   式中，$\tau$ 为延迟时间。
3. 发的

三、疑惑（忽略）

1. 剩余静校正到底位于哪个步骤
2. CMP道集何如形成，所谓24道接收，12次覆盖是什么意思
3. 剩余静校正使用的数据集是什么样，怎么使用

四、算法实现

1. 互相关算法

其基本计算方法为：
$${\varphi }_{GH}(\tau )=\frac{{\sum }_{t=t_1}^{t_2}G(t)H(t+\tau )dt}{{\left[{\sum }_{t=t_1}^{t_2}{G}^2(t)dt{\sum }_{t=t_1}^{t_2}{H}^2(t)dt\right]}^{1/2}}\text{\hspace{1em}(1)}$$
其中，$G(t)$ 和 $H(t)$代表两个地震道，$\varphi$ 表示最终的时移，$\tau$ 为滞后时间。

由前面的内容可知，剩余静校正量存在于每一个炮点/接收点中，而最终体现在每一道数据中，每道数据都有其自己的剩余静校正量。

因此，互相关算法就被用来解决剩余静校正问题。其基本思路是用两个地震道 $G(t)$ 和 $H(t)$ 进行互相关计算（互相关计算方法见式（1），最终可得到一条互相关曲线，如下图所示。图中Point表示实际计算出的值，Curve表示经过平滑后的互相关曲线，然后我们选取互相关值最大的地方，作为$G(t)$ 和 $H(t)$之间的时移。

显然，如果$G(t)$ 和 $H(t)$其中有一道是没有剩余静校正量的（也可以叫做模型道），那就可以计算出另一道的剩余静校正量，以此循环，最终可计算出道的剩余静校正量，即完成了校正过程。但这个方法也存在一定的问题，首先是计算复杂度，其次是模型道的获取是困难的，我们无法准确的判断出哪一道是没有剩余静校正量的道。

算法的python代码如下：

def cross_correlation(para_data):
    # 一维数组，用于保存最终的静校正量
    temp_final_static = []
    for episode in range(len(para_data)):
        para_model = para_data[0]
        para_second = para_data[episode]

        # 保存画图所需的数据
        temp_draw = []
        for i in range(2 * T + 1):
            # 保存自变量t
            temp_t = - T + i
            # 保存每一轮的互相关值
            temp_cross_up = 0
            temp_cross_down = 0
            # 开始求互相关的值
            for j in range(STATIC_LEN):
                temp_index = j + temp_t
                # 判断数组是否越界
                if temp_index < 0 or temp_index >= STATIC_LEN:
                    temp_up = 0
                else:
                    temp_up = para_model[j] * para_second[j + temp_t]
                temp_down = ((para_model[j] ** 2) * (para_second[j] ** 2)) ** 0.5
                # print(temp_down)
                # 计算分式上下的值并保存
                temp_cross_up += temp_up
                temp_cross_down += temp_down
            temp_cross = temp_cross_up / temp_cross_down
            temp_draw.append(temp_cross)

        # 根据互相关结果，求互相关最大值，并且求出静校正量
        temp_static_amount = np.argmax(temp_draw) - 5
        temp_final_static.append(temp_static_amount)

        # print(np.argmax(temp_draw) - 5)
        # 画图，查看互相关曲线
        # draw_cross_curve(temp_draw)
    return temp_final_static

使用互相关算法计算出的剩余静校正量如下图所示，为了分辨，选取了效果不好的时候作图，其中蓝色代表原来的，黄色代表校正出来的。

如果将该算法施加道模拟的数据中，其效果较为理想，上面的是校正前的，下面是矫正后的。

2. 模拟退火算法

模拟退火算法是一种最优化算法，受启发于晶体的熔融物质的冷却过程，其目标是寻找目标函数的最小值。主要思想为：将待求的参数编入某一个模型，然后每次根据当前的温度（主要参数）设置一个扰动，如果该扰动使得目标函数变小，则承认该扰动，否则就按照某一概率承认该扰动。这里需要注意一点，剩余静校正中的目标函数是最大化能量，所以使用模拟退火需要变为最小化负能量。

在本实验中，承认概率以及温度的设置如式（2）、（3）所示。
$$P(-\Delta E)=e^{\frac{-\Delta E}{KT}}\text{\hspace{1em}(2)}$$
式中，T是温度，k是常数，为了方便起见，k常常为1。

$$T_k=\frac{C}{\log (1+k)}\text{\hspace{1em}(2 - 5 - 2)}$$
其中k代表迭代次数，C是不依赖与k的常量，其值通常很大。

具体实现如下：

def simulate(para_k, para_shot_num, para_data):
    # 保存初始的叠加能量
    init_energy = energy([])
    # 保存能量的变化情况
    temp_energy_trend = []
    for episode in range(para_k):
        # 保存每一轮可接受的数量
        ever_accept_num = 0
        # 保存初始温度
        t = C / (math.log(1 + episode + 1))
        # 保存每一轮的静校正量
        ever_epi_static = np.zeros(STATIC_SIZE)
        for ever_shot in range(para_shot_num):
            # 产生随机数作为推荐剩余静校正量
            amount_static = np.random.randint(-3, 3)
            # 构造静校正量数组，用于计算叠加能量
            temp_static_array = np.zeros(STATIC_SIZE)
            temp_static_array[ever_shot] = amount_static

            # 计算叠加能量并判断
            temp_energy = - fitness_function(para_data, temp_static_array)
            delta_e = temp_energy - init_energy

            # 如果负叠加能量更小则接受，否则再根据温度判断是否接受
            # 一旦判定接受此推荐，则更新剩余静校正数组、初始能量、可接受数量
            if delta_e <= 0:
                ever_epi_static[ever_shot] = amount_static
                init_energy = temp_energy
                ever_accept_num += 1
            else:
                # 计算接受的概率
                temp_p = math.e ** (- delta_e / t)
                temp_num = np.random.random()
                # 如果处于可接收的范围内，则改变推荐的静校正量
                if temp_p > temp_num:
                    ever_epi_static[ever_shot] = amount_static
                    init_energy = temp_energy
                    ever_accept_num += 1

        temp_energy_trend.append(init_energy)
        
        # 设置终止条件
        if ever_accept_num <= 0:
            print('可接受推荐小于0，退出循环')
            return ever_epi_static, temp_energy_trend
    print('轮次结束，退出循环')
    return ever_epi_static, temp_energy_trend

可以看出模拟退火算法流程并不复杂，但和初始温度以及常数C的设置有较大的关系。其校正前后的对比如下图所示，可以看出它有一定的效果，但并不理想，推测可能跟我的参数设置有关：

---

---

---

五、机器学习

该部分为增加知识储备+复习

一、模型评估

1. 数据集划分

1. 留出法：对数据集进行互斥划分。常用$\frac{2}{3}$ ~ $\frac{4}{5}$作为训练，需要进行多次划分，取平均值
2. 交叉验证：将数据集划分为$k$个互斥子集，每次使用$k-1$个子集的并集作为训练集，k次后去平均值，通常$k$取值为10。也需要多次划分，其具体代码如下

def k_cross(clf, para_data):
    """
    10次10折交叉验证
    :param clf: 分类器
    :param para_data: 包含训练集和测试集的所有数据
    :return: None
    """
    final_acc = 0
    for i in range(10):
        np.random.shuffle(para_data)
        temp_data = np.vsplit(para_data, 10)
        temp_data = np.array(temp_data)
        temp_acc = 0
        for k in range(10):
            # 对数据进行拆分
            temp_x_test = temp_data[k][:, :-1]
            temp_data_ = np.delete(temp_data, k, axis=0)
            temp_y_test = np.squeeze(temp_data[k][:, -1:])
            temp_train = np.reshape(temp_data_, (-1, 5))
            temp_x_train = temp_train[:, :-1]
            temp_y_train = np.squeeze(temp_train[:, -1:])

            # 训练并计算测试准确率
            clf.fit(temp_x_train, temp_y_train)
            pred = clf.predict(temp_x_test)
            temp_acc += score(temp_y_test, pred)
        final_acc += temp_acc
        print('第{}次10折验证，准确率为{}'.format(i+1, temp_acc / 10))
    print('k折交叉验证最终测试准确率为：{}'.format(final_acc / 100))

3. 自助法：从$m$个样本的数据集中有放回抽样$m$次，抽样获得的数据用于训练，从未抽到的用于测试。容易引入估计偏差，不如前两者通用
4. 实现了留出法和交叉验证法之后，在计算测试准确率的时候，留出法总是要比交叉验证高一些，即误差更小，效果也就差一些

2. 性能度量

1. 错误率与精度：分类任务常用指标
2. 查准率P、查全率R、F1：关注所有正例有多少被判断正确类似问题
3. ROC与AUC：…
4. 代价敏感错误率：见前面博文

3. 比较检验

暂无

二、线性模型

$$f(\mathbf{x})={\mathbf{w}}^{\mathrm{T}}\mathbf{x}+b\text{\hspace{1em}(1)}$$

1.线性回归

优化目标
$$({\mathbf{w}}^{\ast },b^{\ast })=\underset{(\mathbf{w},b)}{\mathrm{argmin}}\sum _{i=1}^m(y_i-{\mathbf{w}}^{\mathrm{T}}{\mathbf{x}}_i-b{)}^2\text{\hspace{1em}(2)}$$

2. Logistics回归（分类问题）

$$y=\frac{1}{1+e^{-({\mathbf{w}}^{\mathrm{T}}\mathbf{x}+b)}}\text{\hspace{1em}(3)}$$

3. 线性判别分析LDA

1. 主要思想：设法将样例投影到一条直线上，是的同类样例尽可能近（在直线上的距离），异类样例尽可能远（同类异类距离远）
2. 优化目标
   $$\max \frac{{\mathbf{w}}^{\mathrm{T}}({\mu }_0-{\mu }_1)({\mu }_0-{\mu }_1{)}^{\mathrm{T}}\mathbf{w}}{{\mathbf{w}}^{\mathrm{T}}({\Sigma }_0+{\Sigma }_1)\mathbf{w}},\text{\hspace{1em}(4)}$$
   其中，${\mu }_i、{\Sigma }_i、\mathbf{w}$分别表示第$i$类实例的均值向量、协方差矩阵、直线。具体实现如下

def lda(para_class_1, para_class_2):
    """
    线性判别分析实现二分类
    :param para_class_1: 第一类数据
    :param para_class_2: 第二类数据
    :return: None
    """
    # 计算均值，寻找中心点
    mu_1 = np.mean(para_class_1, axis=0)
    mu_2 = np.mean(para_class_2, axis=0)

    # 合并所有样本
    temp_data = np.vstack((para_class_1, para_class_2))
    mu_all = np.mean(temp_data, axis=0)

    # 计算类内离散度矩阵S_w
    n_1, n_2 = para_class_1.shape[0], para_class_2.shape[0]
    s_1, s_2 = 0, 0
    for i in range(n_1):
        s_1 += (para_class_1[i, :] - mu_1).T * (para_class_1[i, :] - mu_1)
    for i in range(n_2):
        s_2 += (para_class_2[i, :] - mu_2).T * (para_class_2[i, :] - mu_2)
    s_w = (n_1 * s_1 + n_2 * s_2) / (n_1 + n_2)

    # 计算类间离散度矩阵S_b
    s_b = (n_1 * (mu_all - mu_1).T * (mu_all - mu_1) + n_2 * (mu_all - mu_2).T * (mu_all - mu_2)) / (n_1 + n_2)

    # 求最大特征值对应的特征向量
    eig_value, eig_vector = np.linalg.eig(np.mat(s_w).I * s_b)
    index_vector = np.argsort(-eig_value)
    n_largest = index_vector[:1]
    w = eig_vector[:, n_largest]
    b = np.dot((mu_1 + mu_2) / 2, w)

    return w, b

4.多分类学习

1. 将二分类直接扩展为多分类
2. 拆解法：每两个类别训练一个分类器最后投票、一对多（一个类别作为正，其他作为负）训练分类器后根据置信度

5. 类别不平衡问题

1. 欠采样
2. 过采样
3. 阈值移动
   略

三、决策树

1.划分选择

构建决策树的关键在于如何选择最优划分属性。

1. 信息熵：使用信息熵，度量样本集合纯度，该值越小样本的纯度越高，若当前样本$D$中第$k$类的比例为$p_k$，则此时的信息熵为：
   $$Ent(D)=-\sum _{k=1}^{|Y|}{p}_k{\log }_2{p}_k\text{\hspace{1em}(5)}$$
2. 由于信息增益率倾向于可取值数目较多的属性，故又提出了增益率，该方法倾向于数目较少的属性
3. 基尼指数：CART决策树使用基尼指数选择划分属性，该值越小，纯度越高，计算方式如下：
   $$Gini(D)=1-\sum _{k=1}^{|Y|}{p}_k^2\text{\hspace{1em}(6)}$$

2.剪枝

决策树有时会将节点划分的过于仔细，导致过拟合，故常使用不同的剪枝策略降低该风险，主要分为预剪枝与后剪枝。

1. 预剪枝：在每一次划分时，判断该划分是否会提升决策树效果，若提升了才执行该划分。但也可能造成欠拟合
2. 后剪枝：构建完决策树之后，自底向上对每个节点考察，判断剪枝该节点是否提升效果，若是则剪枝。但时间开销稍大

3.连续与缺失

1. 连续值处理：常用方法为二分法，取是的信息增益率最大的地方作为划分点。若某个节点使用了>0.25，则子节点仍可使用>0.4
2. 缺失值处理：略

多变量决策树

上述决策树划分时，划分出的边界平行于坐标轴，如果数据较为复杂，则需要n多次划分，开销较大。

1. 若想降低开销，则需要划分一条“曲线”或者“斜线段”，这时不再是为每个节点寻找一个最优属性划分，而是将一个线性分类器作为节点，例如$0.25\ast a_1+-0.36\ast a_2\le 0.58$。

二、使用步骤

1.引入库

代码如下（示例）：

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import warnings
warnings.filterwarnings('ignore')
import  ssl
ssl._create_default_https_context = ssl._create_unverified_context

2.读入数据

代码如下（示例）：

data = pd.read_csv(
    'https://labfile.oss.aliyuncs.com/courses/1283/adult.data.csv')
print(data.head())

该处使用的url网络请求的数据。

---

总结

提示：这里对文章进行总结：
例如：以上就是今天要讲的内容，本文仅仅简单介绍了pandas的使用，而pandas提供了大量能使我们快速便捷地处理数据的函数和方法。